求教Fortran求解简单的线性方程，提示：fail
:@Q唉，用Fortran的高斯消元法解一个10未知量的线性方程，很简单的，就是
ax+by+cz+.....+=0，可是运行以后，提示是fail，失败？
为什么呢，是无解还是很多解呢？
请大家指教！！
下面是具体的方程
-2021.106282+0.8605360981e-4*a51+0.2733567103e-3*a52+0.6871398328e-3*a53+0.1611353844e-2*a54+0.1146461461e-4*a61+0.3650562961e-4*a62+0.9195739610e-4*a63+0.1061626344e-5*a71+0.3387872307e-5*a72+0.8550527491e-5*a73 = 0,
-857.4808972+0.3779405407e-4*a51+0.1187759143e-3*a52+0.2956048332e-3*a53+0.6871398328e-3*a54+0.5015867178e-5*a61+0.1581385748e-4*a62+0.3946772124e-4*a63+0.4628346079e-6*a71+0.1463494083e-5*a72+0.3661992028e-5*a73 = 0,
-339.5389867+0.1554591708e-4*a51+0.4826821982e-4*a52+0.1187759143e-3*a53+0.2733567103e-3*a54+0.2053775622e-5*a61+0.6403018372e-5*a62+0.1581385820e-4*a63+0.1887148741e-6*a71+0.5905826516e-6*a72+0.1463494151e-5*a73 = 0,
-114.7882830+0.5015867157e-5*a51+0.1581385820e-4*a52+0.3946772124e-4*a53+0.9195739610e-4*a54+0.6666193011e-6*a61+0.2107720891e-5*a62+0.5273699370e-5*a63+0.6158938563e-7*a71+0.1952474903e-6*a72+0.4896656633e-6*a73 = 0,
-106.4791792+0.5071755690e-5*a51+0.1554591708e-4*a52+0.3779405407e-4*a53+0.8605360981e-4*a54+0.6666247127e-6*a61+0.2053775528e-5*a62+0.5015867157e-5*a63+0.6096684187e-7*a71+0.1887148610e-6*a72+0.4628345948e-6*a73 = 0,
-45.35079793+0.2053775528e-5*a51+0.6403018372e-5*a52+0.1581385748e-4*a53+0.3650562961e-4*a54+0.2718343318e-6*a61+0.8506225952e-6*a62+0.2107720891e-5*a63+0.2502024843e-7*a71+0.7855948901e-7*a72+0.1952474905e-6*a73 = 0,
-14.18452482+0.6666247127e-6*a51+0.2053775622e-5*a52+0.5015867178e-5*a53+0.1146461461e-4*a54+0.8783104738e-7*a61+0.2718343318e-6*a62+0.6666193011e-6*a63+0.8050144469e-8*a71+0.2502024792e-7*a72+0.6158938431e-7*a73 = 0,
-10.67722920+0.4628345948e-6*a51+0.1463494151e-5*a52+0.3661992028e-5*a53+0.8550527491e-5*a54+0.6158938431e-7*a61+0.1952474905e-6*a62+0.4896656633e-6*a63+0.5696757617e-8*a71+0.1810238466e-7*a72+0.4549493106e-7*a73 = 0,
-4.209751084+0.1887148610e-6*a51+0.5905826516e-6*a52+0.1463494083e-5*a53+0.3387872307e-5*a54+0.2502024792e-7*a61+0.7855948901e-7*a62+0.1952474903e-6*a63+0.2306433509e-8*a71+0.7263906386e-8*a72+0.1810238466e-7*a73 = 0,
-1.313551051+0.6096684187e-7*a51+0.1887148741e-6*a52+0.4628346079e-6*a53+0.1061626344e-5*a54+0.8050144469e-8*a61+0.2502024843e-7*a62+0.6158938563e-7*a63+0.7392847122e-9*a71+0.2306433509e-8*a72+0.5696757617e-8*a73 = 0 如果高斯消去程序本身没有问题，则需要看看系数矩阵的条件数。过高的条件数可以导致高斯法不能运算或给出错误解（误差过大）。对于条件数比较的系数矩阵，可以采用双精度数据格式，情况会好一点。
最简单的高斯消去法程序如下。此程序经过检验没有错误，如果不能计算说明方程本身有问题。这个程序只适用于小规模并且系数矩阵规范的情况，不是很实用。推荐使用Lapack 3 或者 slatec 中的程序模块。这些模块都附带用来检验有效位数的子程序。Lapack 3 或者 slatec 都是免费的，在网上可以找到。
!==============================================================================
function sl(a,b)
! linear solver using gaussian elimination for full matrix
! solve a*x=b,a: n-by-n matrix, b: n-by-1 array
! note: this function is only efficient on small matrix!
! ------------------------------------------------------------------------------
implicit none
! ------------------------------------------------------------------------------
! declare dummy variables
double precision, dimension(:,:), intent(in) :: a
double precision, dimension(:), intent(in) :: b
double precision, dimension(size(b,1)):: sl
! ------------------------------------------------------------------------------
! declare local variables
integer:: i,j,k,n
double precision:: sum,xx
double precision, dimension(size(a,1),size(a,2)):: aa
double precision, dimension(size(b,1)):: bb
! ------------------------------------------------------------------------------
sl=0.0d0
n=size(b,1)
aa=a
bb=b
! convert to upper triangular form
do k = 1,n-1
if (dabs(aa(k,k)).gt.1.d-6) then
do i = k+1, n
xx = aa(i,k)/aa(k,k)
do j = k+1, n
aa(i,j) = aa(i,j) - aa(k,j)*xx
end do
bb(i) = bb(i) - bb(k)*xx
end do
else
write (*,'("zero pivot found in line ",i5)') k
stop
end if
end do
! ------------------------------------------------------------------------------
! back substitution
do i = n,1,-1
sum = bb(i)
if (i.lt.n) then
do j= i+1,n
sum = sum - aa(i,j)*bb(j)
end do
end if
bb(i) = sum/aa(i,i)
end do
sl=bb
! ------------------------------------------------------------------------------
end function sl
! ------------------------------------------------------------------------------
! note: to use the function, add the following interface to the caller
! interface
! function sl(a,b)
! double precision, dimension(:,:), intent(in) :: a
! double precision, dimension(:), intent(in) :: b
! double precision, dimension(size(b,1)):: sl
! end function sl
! end interface
! ==============================================================================
页:
[1]