消主元解方程的程序问题
subroutine ch(iz)parameter nub0=2400
common/a4/rm(nub0,nub0),x(nub0),u(nub0)
write(*,*)rm
do 20 ii=1,iz-1
c=1.0/rm(ii,ii)
do 15 ij=ii,iz
15 rm(ii,ij)=rm(ii,ij)*c
u(ii)=u(ii)*c
do 20 jj=ii+1,iz
do 30 kk=ii+1,iz
rm(jj,kk)=rm(jj,kk)-rm(ii,kk)*rm(jj,ii)
30 continue
u(jj)=u(jj)-u(ii)*rm(jj,ii)
20 continue
x(iz)=u(iz)/rm(iz,iz)
do 10 ii=iz-1,1,-1
as1=0.0
do 40 jj=ii+1,iz
as1=as1+x(jj)*rm(ii,jj)
40 continue
x(ii)=u(ii)-as1
10 continue
do 55 ii=1,iz-1
write(*,*)u(ii)
55 continue
return
end
这是一个消主元解微分方程的程序 主程序里rm已经设定为系数矩阵了，为什么我读u（ii）能读出但是读x（ii）就结果就是NAN呢
你的代码风格太古老了，应该更新一下了。
我对代码做了现代化风格的改造，主要错误是应该先对 u 操作，然后对 rm 操作。
给你的建议：
1.要学习 Fortran95 或以上的新语法
2.放弃 Compaq Fortran，放弃固定格式，放弃 Common！！
Program Main
Implicit None
Integer , Parameter :: nub0 = 5
Real :: rm , x , u , rma(nub0,nub0)
Common /a4/rm(nub0, nub0), x(nub0), u(nub0)
call RANDOM_SEED()
call RANDOM_NUMBER(rm)
call RANDOM_NUMBER(u)
rma = rm
write(*,'("rm0:",/,5(5es14.6,/))') rma!//原始系数矩阵
write(*,'("U0:",/,5es14.6)') u!//原始右端向量
call ch( nub0 )
write(*,'("X:",/,5es14.6)') x!//解
write(*,'("rm_up:",/,5(5es14.6,/))') rm!//上三角化
write(*,'("U_Check:",/,5es14.6)') matmul( rma , reshape(x,) )!//验证
End Program Main
Subroutine ch(z)
Implicit None
Integer , Parameter :: nub0 = 5
Real :: rm , x , u
integer :: i , z , j
Common /a4/rm(nub0, nub0), x(nub0), u(nub0)
Do i = 1, z
u(i)=u(i) / rm(i,i)
rm(i,i:z) = rm(i,i:z) / rm(i,i)
Do j = i + 1, z
u(j)=u(j)-u(i)*rm(j,i)
rm(j,i:z) = rm(j,i:z) - rm(i,i:z) * rm(j,i)
End Do
End Do
x(z) = u(z)/rm(z, z)
Do i = z - 1, 1, -1
x(i) = u(i) - sum( x(i+1:z) * rm(i,i+1:z) )
End Do
End Subroutine ch
页:
[1]