消主元解方程的程序问题

Fortran会长大

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

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呢

fcode

你的代码风格太古老了，应该更新一下了。

我对代码做了现代化风格的改造，主要错误是应该先对 u 操作，然后对 rm 操作。

给你的建议：
1.要学习 Fortran95 或以上的新语法
2.放弃 Compaq Fortran，放弃固定格式，放弃 Common！！

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

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,[nub0,1]) )!//验证
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