解方程组代码问题求助

806940059

用高斯消元法解线性方程组，为九元一次方程式组，遇到结果为NaN的情况，代码如下：

[Fortran] 纯文本查看 复制代码

Program linear_equation
Implicit none
integer::i,j,k,imax,t
real::max,n
!用矩阵(实型动态数组)将线性方程组表示出来
real,dimension(:,:),allocatable::a,m
real,dimension(:),allocatable::x
!线性方程组的元数t
t=9
!给动态数组分配内存
allocate(a(t+1,t+2),m(t,t+1),x(t))
!依次按行读入线性方程组的增广矩阵a
open(1,file='zengguangjuzhen1.txt')
read(1,*)((a(i,j),j=1,t+1),i=1,t)
!列出增广矩阵
!print*,'所读取增广矩阵为：'
!do j=1,t+1
!print*,a(i,j)
!end do
!下面选取每列的最大列主元素
do k=1,t-1
max=abs(a(k,k))
imax=k
Do i=k+1,t
if (abs(a(i,k))>max) then
max=abs(a(i,k))
imax=i
end if
end do
!将最大列元素所在的行与第K行进行交换
Do j=k,t+1
m(k,j)=a(k,j)
a(k,j)=a(imax,j)
a(imax,j)=m(k,j)
end do
!对方程组按X1,X2,…,Xt的顺序进行依次消元，将增广矩阵a化为上三角矩阵
Do i=k+1,t
m(i,k)=a(i,k)/a(k,k)
do j=k+1,t+1
a(i,j)=a(i,j)-m(i,k)*a(k,j)
end do
end do
end do
!guass主元消去法的第一步消去至此结束
!先计算出xt
x(t)=a(t,t+1)/a(t,t)
!下面进行回代计算
Do k=t-1,1,-1
n=0
do j=t,k+1,-1
n=n+a(k,j)*x(j)
end do
x(k)=(a(k,t+1)-n)/a(k,k)
end do
print*,'线性方程组的解为：'
do i=1,t
print*,x(i)
end do
end program linear_equation

所读取的增广矩阵为：
1 -1 -1  1  1/3  1/3 -1/3 -1/3  1/9  26
1  0 -1  0 -2/3  1/3  2/3    0   -2/9  28
1  1 -1 -1  1/3  1/3 -1/3   1/3  1/9  20
1 -1  0  0  1/3  -2/3   0    2/3  -2/9  30
1  0  0  0 -2/3  -2/3  0       0    4/9  35
1  1  0  0  1/3 -2/3  0    -2/3   -2/9  27
1 -1  1 -1  1/3  1/3  1/3 -1/3  1/9   23
1  0  1  0 -2/3  1/3 -2/3    0   -2/9  24
1  1  1  1  1/3  1/3  1/3  1/3  1/9   20

结果为：
线性方程组的解为：
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
NaN
麻烦给看一下错在哪了，有劳了！谢谢！

806940059

补充：正确结果应该是25.8889
   -2.0000
   -1.1667
    0.7500
   -4.6667
   -7.1667
    1.2500
   -0.7500
    2.7500

vvt

我这里解出来是

线性方程组的解为：
   25.88889
  -2.000000
  -1.166667
  0.7500000
  -4.666667
  -7.166667
   1.250000
-0.7500000
   2.750000

如果你在文件中输入：
1 -1 -1  1  1/3  1/3 -1/3 -1/3  1/9  26
1  0 -1  0 -2/3  1/3  2/3    0   -2/9  28
1  1 -1 -1  1/3  1/3 -1/3   1/3  1/9  20
1 -1  0  0  1/3  -2/3   0    2/3  -2/9  30
1  0  0  0 -2/3  -2/3  0       0    4/9  35
1  1  0  0  1/3 -2/3  0    -2/3   -2/9  27
1 -1  1 -1  1/3  1/3  1/3 -1/3  1/9   23
1  0  1  0 -2/3  1/3 -2/3    0   -2/9  24
1  1  1  1  1/3  1/3  1/3  1/3  1/9   20

这是无法被识别的，你需要输入

1.000         -1.000         -1.000         1.000         0.333         0.333         -0.333         -0.333         0.111         26.000
1.000         0.000         -1.000         0.000         -0.667         0.333         0.667         0.000         -0.222         28.000
1.000         1.000         -1.000         -1.000         0.333         0.333         -0.333         0.333         0.111         20.000
1.000         -1.000         0.000         0.000         0.333         -0.667         0.000         0.667         -0.222         30.000
1.000         0.000         0.000         0.000         -0.667         -0.667         0.000         0.000         0.444         35.000
1.000         1.000         0.000         0.000         0.333         -0.667         0.000         -0.667         -0.222         27.000
1.000         -1.000         1.000         -1.000         0.333         0.333         0.333         -0.333         0.111         23.000
1.000         0.000         1.000         0.000         -0.667         0.333         -0.667         0.000         -0.222         24.000
1.000         1.000         1.000         1.000         0.333         0.333         0.333         0.333         0.111         20.000

806940059

vvt 发表于 2014-5-6 21:20
我这里解出来是

如果你在文件中输入：

谢谢！确实是这样！十分感谢！