| 查看: 1689 | 回复: 10 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【求助】fortran出错 编译通过 运行溢出
|
|||
|
subroutine adams(t,h,m,y,n,z,f,d,eps,b,e,s,g,kp1,kp2) external f,rkt dimension y(n),d(n),z(n,m),b(4,n),e(n),s(n),g(n) real*8 kp1,kp2 a=t do 10 i=1,n 10 z(i,1)=y(i) call f(t,y,n,d,kp1,kp2) do 20 i=1,n 20 b(1,i)=d(i) do 50 i=2,4 if (i.le.m)then t=a+(i-1)*h call rkt(t,h,y,n,f,eps,d,e,s,g,kp1,kp2) do 30 k=1,n 30 z(k,i)=y(k) call f(t,y,n,d,kp1,kp2) do 40 k=1,n 40 b(i,k)=d(k) end if 50 continue do 100 i=5,m do 60 j=1,n q=55*b(4,j)-59*b(3,j)+37*b(2,j)-9*b(1,j) y(j)=z(j,i-1)+h*q/24.0 b(1,j)=b(2,j) b(2,j)=b(3,j) b(3,j)=b(4,j) 60 continue t=a+(i-1)*h call f(t,y,n,d,kp1,kp2) do 70 k=1,n 70 b(4,k)=d(k) do 80 j=1,n q=9*b(4,j)+19*b(3,j)-5*b(2,j)+b(1,j) z(j,i)=z(j,i-1)+h*q/24.0 y(j)=z(j,i) 80 continue do 90 k=1,n 90 b(4,k)=d(k) 100 continue t=a return end subroutine rkt(t,h,y,n,f,eps,d,b,c,g,kp1,kp2) dimension y(n),d(n),a(4),b(n),c(n),g(n) real*8 kp1,kp2 h1=h m=1 p=1+eps x=t do 10 i=1,n 10 c(i)=y(i) 20 if (p.ge.eps)then a(1)=h1/2.0 a(2)=a(1) a(3)=h1 a(4)=h1 do 30 i=1,n g(i)=y(i) y(i)=c(i) 30 continue dt=h/m t=x do 80 j=1,m call f(t,y,n,d,kp1,kp2) do 40 i=1,n 40 b(i)=y(i) do 60 k=1,3 do 50 i=1,n y(i)=y(i)+a(k)*d(i) b(i)=b(i)+a(k+1)*d(i)/3.0 50 continue tt=t+a(k) call f(tt,y,n,d,kp1,kp2) 60 continue do 70 i=1,n 70 y(i)=b(i)+h1*d(i)/6.0 t=t+dt 80 continue p=0.0 do 90 i=1,n q=abs(y(i)-g(i)) if(q.gt.p) p=q 90 continue h1=h1/2.0 m=m+m goto 20 endif t=x return end program madams c ----------------------------------------------------------------------- c this rpogram is used to solve a systwem of differential equagtion c by using adams method c ------------------------------------------------------------------- external f dimension y(11),d(11),z(11,500000),b(4,11),e(11),s(11),g(11) real*8 wucha,kp1,kp2,t,tao,variable,Parameters open(26,file='variable.dat') open(28,file='Parameters.dat') open(10,file='wucha.dat') h=0.001d0 y(1)=0.01d0 y(2)=0.02d0 y(3)=0.03d0 y(4)=0.4d0 y(5)=0.2d0 y(6)=0.8d0 y(7)=0.8d0 y(8)=2.5d0 y(9)=0.9d0 y(10)=6d0 y(11)=1d0 m=500000 n=11 eps=0.00001d0 do kp1=1.50d0,1.60d0,0.1d0 do kp2=0.00d0,0.50d0,0.1d0 call adams(a,h,m,y,n,z,f,d,eps,b,e,s,g,kp1,kp2) do i=1,m t=(i-1)*h detae=((z(1,i)-z(4,i))**2.0+(z(2,i)+z(5,i))**2.0+ $ (z(3,i)-z(6,i))**2.0)**0.5 ! write(26,'(2f15.6)')t,detae ! write(28,'(6f15.6)')t,z(7,i),z(8,i),z(9,i),z(10,i),z(11,i) if(detae>50)then write(10,*)kp1,kp2,1 write(*,*)kp1,kp2,1 goto 20 elseif(detae>1.0d-3)then tao=t endif enddo if(tao<=50)then write(10,*)kp1,kp2,0 write(*,*)kp1,kp2,0 else write(10,*)kp1,kp2,1 write(*,*)kp1,kp2,1 endif 20 continue enddo enddo close(2) end subroutine f(t,y,n,d,kp1,kp2) c ---------------------------------------------------------------------- ! real*8 kp1,kp2 dimension y(n),d(n) double precision ex,ey,ez,u1,kp1,kp2 cccccccccccccc ex=y(1)-y(4) ey=y(2)+y(5) ez=y(3)-y(6) c u1=ey-ez+y(8)*(y(1)+y(4))*ex-y(10)*(y(1)+y(4))*ey+4d0*y(11)*ez c u2=0d0 cc u1=ey+y(8)*(y(1)+y(4))*ex-y(10)*(y(1)+y(4))*ey cc u2=-ex+4d0*y(11)*ex u1=ey+y(8)*(y(1)+y(4))*ex-y(10)*(y(1)+y(4))*ey+4d0*y(11)*ez u2=-ex ! do kp1=0.00d0,10.0d0,0.01d0 ! do kp2=0.00d0,10.0d0,0.01d0 ! kp1=2.0d0 ! kp2=0.5d0 cccccccccccccc d(1)=y(2)-1d0*y(1)**3+3d0*y(1)**2-y(3)+3d0 d(2)=1d0-5d0*y(1)**2-y(2) d(3)=0.006d0*(4d0*(y(1)+1.6d0)-y(3)) cccccccccccccccccccccccc d(4)=y(5)-y(7)*y(4)**3+y(8)*y(4)**2-y(6)+3d0+u1 d(5)=y(9)-y(10)*y(4)**2-y(5) d(6)=y(11)*(4d0*(y(4)+1.6d0)-y(6))+u2 cccccccccccccccccccccccccccccccccc d(7)=-kp1*kp2*y(1)**3*ex d(8)=kp1*kp2*y(1)**2*ex d(9)=kp1*kp2*ey d(10)=-kp1*kp2*y(1)**2*ey d(11)=kp1*kp2*(4d0*(y(1)+1.6d0)-y(3))*ez ! enddo ! enddo return end 高手看看 这个程序 出错了 为什么 溢出了！！  |
回复此楼» 猜你喜欢
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有8人回复
-
体制内长辈说体制内绝大部分一辈子在底层，如同你们一样大部分普通教师忙且收入低
已经有13人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有8人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有8人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有11人回复
-
为什么中国大学工科教授们水了那么多所谓的顶会顶刊，但还是做不出宇树机器人？
已经有8人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有5人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有6人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有9人回复
-
售SCI一区文章，我:8 O5 51O 54,科目齐全,可+急
已经有4人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 【求助】求fortra程序
已经有4人回复
9楼2010-09-15 09:23:21
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
2楼2010-09-14 20:22:34
maomao1210
金虫 (正式写手)
- 程序强帖: 5
- 应助: 2 (幼儿园)
- 金币: 1432.8
- 散金: 242
- 红花: 16
- 沙发: 1
- 帖子: 991
- 在线: 442.1小时
- 虫号: 253215
- 注册: 2006-05-20
- 性别: MM
- 专业: 考古理论
★
余泽成(金币+1):谢谢参与应助！ 2010-09-15 13:32:24
余泽成(金币+1):谢谢参与应助！ 2010-09-15 13:32:24
|
do kp1=1.50d0,1.60d0,0.1d0 do kp2=0.00d0,0.50d0,0.1d0 call adams(a,h,m,y,n,z,f,d,eps,b,e,s,g,kp1,kp2) do i=1,m t=(i-1)*h detae=((z(1,i)-z(4,i))**2.0+(z(2,i)+z(5,i))**2.0+ $ (z(3,i)-z(6,i))**2.0)**0.5 ! write(26,'(2f15.6)')t,detae ! write(28,'(6f15.6)')t,z(7,i),z(8,i),z(9,i),z(10,i),z(11,i) if(detae>50)then write(10,*)kp1,kp2,1 write(*,*)kp1,kp2,1 goto 20 elseif(detae>1.0d-3)then tao=t endif enddo if(tao<=50)then write(10,*)kp1,kp2,0 write(*,*)kp1,kp2,0 else write(10,*)kp1,kp2,1 write(*,*)kp1,kp2,1 endif 20 continue enddo enddo 在循环的时候，=((z(1,i)-z(4,i))**2.0+(z(2,i)+z(5,i))**2.0+ $ (z(3,i)-z(6,i))**2.0 是一个无穷大或者无穷小的数字，cvf不识别，不予计算，就提示溢出了。 gfortran 编译计算的时候，即使是NaN 也会继续计算不提示任何错误。 |
3楼2010-09-14 22:22:30
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
4楼2010-09-14 22:43:47
