Fortran Coder

查看: 218|回复: 1

[线性代数] 消主元解方程的程序问题

[复制链接]

1

帖子

1

主题

0

精华

新人

F 币
15 元
贡献
4 点
发表于 2018-6-30 10:25:13 | 显示全部楼层 |阅读模式
[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呢

回复

使用道具 举报

1314

帖子

12

主题

5

精华

论坛跑堂

Fcode跑堂伙计

F 币
156 元
贡献
219 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2018-6-30 11:51:01 | 显示全部楼层
你的代码风格太古老了,应该更新一下了。

我对代码做了现代化风格的改造,主要错误是应该先对 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


您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2018-11-16 22:32

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表