论坛 › Fortran语法讨论 › 解方程组代码问题求助

806940059 发表于 2014-5-6 19:47:19

解方程组代码问题求助

用高斯消元法解线性方程组，为九元一次方程式组，遇到结果为NaN的情况，代码如下：
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 -111/31/3 -1/3 -1/31/926
10 -10 -2/31/32/3    0   -2/928
11 -1 -11/31/3 -1/3   1/31/920
1 -1001/3-2/3   0    2/3-2/930
1000 -2/3-2/30       0    4/935
11001/3 -2/30    -2/3   -2/927
1 -11 -11/31/31/3 -1/31/9   23
1010 -2/31/3 -2/3    0   -2/924
11111/31/31/31/31/9   20

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

806940059 发表于 2014-5-6 20:02:58

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

vvt 发表于 2014-5-6 21:20:05

我这里解出来是
线性方程组的解为：
   25.88889
-2.000000
-1.166667
0.7500000
-4.666667
-7.166667
   1.250000
-0.7500000
   2.750000
如果你在文件中输入：
1 -1 -111/31/3 -1/3 -1/31/926
10 -10 -2/31/32/3    0   -2/928
11 -1 -11/31/3 -1/3   1/31/920
1 -1001/3-2/3   0    2/3-2/930
1000 -2/3-2/30       0    4/935
11001/3 -2/30    -2/3   -2/927
1 -11 -11/31/31/3 -1/31/9   23
1010 -2/31/3 -2/3    0   -2/924
11111/31/31/31/31/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 发表于 2014-5-7 16:10:58

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

如果你在文件中输入：


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

查看完整版本: 解方程组代码问题求助