| 查看: 1716 | 回复: 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 高手看看 这个程序 出错了 为什么 溢出了！！  |
回复此楼» 猜你喜欢
308求调剂
已经有4人回复
-
NSFC申报书里申请人简历中代表性论著还需要在申报书最后的附件里面再上传一遍吗
已经有14人回复
-
材料与化工一志愿南昌大学327求调剂推荐
已经有6人回复
-
化学调剂0703
已经有7人回复
-
327求调剂
已经有11人回复
-
调剂
已经有8人回复
-
梁成伟老师课题组欢迎你的加入
已经有7人回复
-
伙伴们，祝我生日快乐吧
已经有24人回复
-
中科院材料273求调剂
已经有3人回复
-
材料工程专硕274一志愿211求调剂
已经有5人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 【求助】求fortra程序
已经有4人回复
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
maomao1210
金虫 (正式写手)
- 程序强帖: 5
- 应助: 2 (幼儿园)
- 金币: 1432.8
- 散金: 242
- 红花: 16
- 沙发: 1
- 帖子: 991
- 在线: 442.1小时
- 虫号: 253215
- 注册: 2006-05-20
- 性别: MM
- 专业: 考古理论
snoopyzhao ，嘿嘿，
5楼2010-09-14 22:57:57
6楼2010-09-15 08:54:13
maomao1210
金虫 (正式写手)
- 程序强帖: 5
- 应助: 2 (幼儿园)
- 金币: 1432.8
- 散金: 242
- 红花: 16
- 沙发: 1
- 帖子: 991
- 在线: 442.1小时
- 虫号: 253215
- 注册: 2006-05-20
- 性别: MM
- 专业: 考古理论
7楼2010-09-15 09:07:29
8楼2010-09-15 09:22:15
9楼2010-09-15 09:23:21
maomao1210
金虫 (正式写手)
- 程序强帖: 5
- 应助: 2 (幼儿园)
- 金币: 1432.8
- 散金: 242
- 红花: 16
- 沙发: 1
- 帖子: 991
- 在线: 442.1小时
- 虫号: 253215
- 注册: 2006-05-20
- 性别: MM
- 专业: 考古理论
10楼2010-09-15 09:48:19
