声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2782|回复: 1

[Fortran] 龙格库塔法求解微分方程组求助!!!

[复制链接]
发表于 2007-5-10 11:22 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
应用四阶龙格库踏法求解微分方程组,其中微分方程组中有一项要应用数值积分计算,程序如下:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!main!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        EXTERNAL F
        DIMENSION Y(3),D(3),Z(3,11),B(3)
        DOUBLE PRECISION Y,D,Z,T,H,B
        T=0.0
        Y(1)=-1.0
        Y(2)=0.0
        Y(3)=1.0
        H=0.01
        M=3
        N=11
        CALL RKT2(T,Y,M,H,N,Z,F,D,B)
        WRITE(*,*)
        DO 10 I=1,N
          T=(I-1)*H
          WRITE(*,50) T
          WRITE(*,100) (Z(J,I),J=1,M)
          WRITE(*,*)
10        CONTINUE
50        FORMAT(1X,'T=',F7.3)
100        FORMAT(1X,'Y(1)=',D13.6,3X,'Y(2)=',D13.6,3X,'Y(3)=',D13.6)
        END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!龙格库塔法需要的外部函数!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        SUBROUTINE F(T,Y,M,D)
        DIMENSION Y(m),D(m)
        DOUBLE PRECISION Y,D,T,f1
      real*8 A0,B0
      real EPS  !龙贝格数值积分的误差
        parameter(EPS=0.0001)
      
        A0=0.D0
        B0=5
      call ROMB(A0,B0,ff,EPS,f1)

        D(1)=Y(2)+f1
        D(2)=-Y(1)
        D(3)=-Y(3)
        RETURN
        END

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!龙格-库塔法子程序!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      SUBROUTINE RKT2(T,Y,M,H,N,Z,F,D,B)
        DIMENSION Y(M),D(M),Z(M,N),A(4),B(M)
        DOUBLE PRECISION Y,D,Z,A,B,T,H,X,TT
        A(1)=H/2.0
        A(2)=A(1)
        A(3)=H
        A(4)=H
        DO 5 I=1,M
5        Z(I,1)=Y(I)
        X=T
        DO 100 J=2,N
          CALL F(T,Y,M,D)
          DO 10 I=1,M
10          B(I)=Y(I)
          DO 30 K=1,3
            DO 20 I=1,M
              Y(I)=Z(I,J-1)+A(K)*D(I)
              B(I)=B(I)+A(K+1)*D(I)/3.0
20            CONTINUE
            TT=T+A(K)
            CALL F(TT,Y,M,D)
30          CONTINUE
          DO 40 I=1,M
40          Y(I)=B(I)+H*D(I)/6.0
          DO 50 I=1,M
50          Z(I,J)=Y(I)
          T=T+H
100        CONTINUE
        T=X
        RETURN
        END

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!龙贝格数值积分子程序!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE ROMB(A,B,F,EPS,T)
        DIMENSION Y(10)
        DOUBLE PRECISION A,B,F,T,Y,H,P,S,Q
        H=B-A
        Y(1)=H*(F(A)+F(B))/2.0
        M=1
        N=1
10        P=0.0
        DO 20 I=0,N-1
20        P=P+F(A+(I+0.5)*H)
        P=(Y(1)+H*P)/2.0
        S=1.0
        DO 30 K=1,M
          S=4*S
          Q=(S*P-Y(K))/(S-1)
          Y(K)=P
          P=Q
30        CONTINUE
        IF ((ABS(Q-Y(M)).GE.EPS).AND.(M.LE.9)) THEN
          M=M+1
          Y(M)=Q
          N=N+N
          H=H/2.0
          GOTO 10
        END IF
        T=Q
        RETURN
        END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!龙贝格数值积分需要的外部函数!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      real*8 function ff(t)
        real*8 t
      
        ff=abs(sin(1/t)+cos(t))
        end


运行之后有下面的警告:Compiling Fortran...
F:\fortran\werewrrer\ewee.for
F:\fortran\werewrrer\ewee.for(37) : Warning: In the call to ROMB, actual argument #3 does not match the type and kind of the corresponding dummy argument.
      call ROMB(A0,B0,ff,EPS,f1)
----------------------^

ewee.obj - 0 error(s), 1 warning(s)

请高手给予指点...万分感谢.

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2007-5-13 07:49 | 显示全部楼层
统一一下数据类型看看
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-24 08:48 , Processed in 0.124791 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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