Fortran Coder

查看: 4447|回复: 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呢

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2033

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1641 元
贡献
709 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

沙发
发表于 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


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

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-12-26 05:05

Powered by Tencent X3.4

© 2013-2024 Tencent

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