解方程组代码问题求助

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

用高斯消元法解线性方程组，为九元一次方程式组，遇到结果为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 · 发表于 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 -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 · 发表于 2014-5-7 16:10:58

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

如果你在文件中输入：

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

[求助] 解方程组代码问题求助

规矩勋章

元老勋章

新人勋章

水王勋章

热心勋章