龙格库塔法求解微分方程组求助!!!
应用四阶龙格库踏法求解微分方程组,其中微分方程组中有一项要应用数值积分计算,程序如下:!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!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)
请高手给予指点...万分感谢. 统一一下数据类型看看
页:
[1]