| 查看: 515 | 回复: 1 | |||
| 当前主题已经存档。 | |||
shaoww木虫 (小有名气)
|
[交流]
[求助]用fortran数值解常微分方程遇到的问题（金币20有劳版主）
|
||
|
要解的是两个二阶的常微分方程组成的方程组，先降阶化为四个方程，然后用积分一步的变步长龙格－库塔法，得到的结果不理想，用matlab模拟的结果是振荡衰减的，但用fortran得到的只有振荡，而无衰减。初次学习用fortran写程序，不知道是程序写的问题还是解微分方程的方法的问题，朋友们指点指点 !积分一步的变步长龙格-库塔法解二阶常微分方程组=============== module parameters DOUBLE PRECISION c,k0,r13,r33,no,ne,er,ei,e0,id,A0,dm,kz,kzi,kx1,kx2,kx3,a,p,b,q, A3x,A3z, A2x,B2x,A2z,B2z,A1x,C1x,A1z,C1z end module parameters ！============================================= PROGRAM main use parameters implicit none integer, parameter::L=60000,M=4; !m为方程组个数，L为步数 EXTERNAL F DIMENSION Y(M),D(M),U(M),V(M),G(M),E(M) DOUBLE PRECISION:: Y,D,U,V,G,E,Z,T,H !Y初值，T初始时刻，H步长 REAL EPS !eps为精度 integer I,J !循环节控制 OPEN(1,STATUS='UNKNOWN',FILE='f:\program files\matlab\work\SPRSP\output\ GRKT-2 \output.TXT') !========付给每个参数值========================== c=6.624D-40 k0=0.992918032108D7 r13=6.7D-11 r33=1.34D-9 no=2.3117D0 ne=2.2987D0 er=-18.0D0 ei=0.7D0 e0=1.0D0 id=1.0D0 A0=300.0D0 dm=4.0D-8 kz=1.070D7 kzi=(0.5*k0**2*ei*e0**2)/(kz*((e0+er)**2+ei**2)) kx1=sqrt((k0*ne)**2-kz**2) kx2=sqrt(kz**2-k0**2*er-kzi**2) kx3=sqrt(kz**2-k0**2*e0) a=2.0*(k0**2)*(ne**4)*r33*c p=(k0**2)*(ne**2)-kz**2 b=2.0*(k0**2)*(no**4)*r13*c q=(k0**2)*(no**2)-kz**2 A3x=A0 A3z=kx3*A3x/kz A2x=0.5*(e0/er-kx3/kx2)*exp(-kx2*dm)*A3x B2x=0.5*(e0/er+kx3/kx2)*exp(kx2*dm)*A3x A2z=0.5*(1-e0*kx2/(er*kx3))*exp(-kx2*dm)*A3z B2z=0.5*(1+e0*kx2/(er*kx3))*exp(kx2*dm)*A3z A1x=er*(A2x+B2x)/(ne**2) C1x=no**2*kx2*(A2x-B2x)/(ne**2) A1z=A2z+B2z C1z=-kx1**2*er*(A2z-B2z)/(ne**2*kx2) !设置程序的初值 T=0.0D0 H=0.00000001D0 Y(1)=A1x Y(2)=C1x Y(3)=A1z Y(4)=C1z EPS=1.0D-04 !======================================= DO 10 I=1,L CALL GRKT2(T,H,Y,M,F,EPS,D,U,V,G,E) T=T+H WRITE(1,*) T,Y(1),Y(3) WRITE(1,*) 10 CONTINUE END PROGRAM main !=========================================== SUBROUTINE F(T,Y,M,D) use parameters INTEGER M DOUBLE PRECISION:: T,Y(M),D(M) !二阶方程组降阶化为四个方程，y(1),y(3)表示所需要求的f，y(2),y(4)表示f' D(1)=Y(2) D(2)=-a*y(1)*(y(1)*y(2)+y(3)*y(4))/(id+y(1)**2.0+y(3)**2.0)-p*y(1) D(3)=Y(4) D(4)=-b*y(3)*(y(1)*y(2)+y(3)*y(4))/(id+y(1)**2.0+y(3)**2.0)-q*y(3) ENDSUBROUTINE !=============================== SUBROUTINE GRKT2(T,H,Y,M,F,EPS,D,U,V,G,E) DIMENSION Y(M),D(M),A(4),U(M),V(M),G(M),E(M) DOUBLE PRECISION Y,D,A,U,V,G,T1,H1,X,HH,E HH=H N=1 P=1+EPS X=T DO 5 I=1,M 5 V(I)=Y(I) 10 IF (P.GE.EPS) THEN A(1)=HH/2.0 A(2)=A(1) A(3)=HH A(4)=HH DO 20 I=1,M G(I)=Y(I) Y(I)=V(I) 20 CONTINUE DT=H/N T=X DO 100 J=1,N CALL F(T,Y,M,D) DO 30 I=1,M E(I)=Y(I) 30 U(I)=Y(I) DO 50 K=1,3 DO 40 I=1,M Y(I)=E(I)+A(K)*D(I) U(I)=U(I)+A(K+1)*D(I)/3.0 40 CONTINUE TT=T+A(K) CALL F(T,Y,M,D) 50 CONTINUE DO 60 I=1,M 60 Y(I)=U(I)+HH*D(I)/6.0 T=T+DT 100 CONTINUE P=0.0 DO 110 I=1,M Q=ABS(Y(I)-G(I)) IF (Q.GT.P)P=Q 110 CONTINUE HH=HH/2.0 N=N+N GOTO 10 END IF T1=X RETURN END [ Last edited by shaoww on 2008-4-16 at 21:00 ] |
回复此楼» 猜你喜欢
怎么看青基中了没有啊
已经有6人回复
-
国社科又开始会评了，不知道这次命运如何
已经有11人回复
-
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有3人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有3人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
有多少人是今天查系统知道结果的？
已经有19人回复
-
基金未中，这种答复是模板吗？
已经有7人回复
-
哪位高人中了，把查询到的截图贴出来让我看看，让我长长见识
已经有7人回复
-
面上意见出来了
已经有23人回复
shaoww
木虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 3678.3
- 红花: 1
- 帖子: 147
- 在线: 41.6小时
- 虫号: 369505
- 注册: 2007-05-12
- 性别: MM
- 专业: 光学
2楼2008-04-18 19:16:52
投票: