用4阶龙格－库塔法解微分方程

xj2070

MATLAB provides a package tool funfun, which, besides
the two classical Runge-Kutta Fehlberg methods, RK23 (second-order and
third-order pair) and RK45 (fourth-order and fifth-order pair), also implements
other methods suitable for solving stiff problems. These methods are
derived from BDF methods and are included in the MATLAB
program ode15s,

so you needn't programming it yourself
----------------------------------------------------------------------------------------------
by the way,the code of the fourth- orde RK algorthim in fortran is presented
as follows:
------------------------------------------------------------------------------------------------

SUBROUTINE XRKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend)
PARAMETER(Big=1E13)
INTEGER I,J,K,M,Mend
REAL A,B,T,Tol,Y,Y0
REAL Br,H,Err,Hmin,Hmax,K1,K2,K3,K4,K5,K6,TJ,YJ
REAL A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,S
REAL A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5
REAL M1,M3,M4,M5,M6,Y1,Y2,Y3,Y4,Y5,Y6,Ygood,Ynew
DIMENSION T(1:1000),Y(1:1000)
EXTERNAL F
DATA A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,
A A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5
B/ 0.25, 0.25, 0.375, 0.09375, 0.28125, 0.9230769231,
C 0.8793809741, -3.277196177, 3.320892126, 1.0,
D 2.032407407, -8.0, 7.173489279, -0.2058966862, 0.5,
E -0.2962962963, 2.0, -1.381676413, 0.4529727096, -0.275,
F 0.002777777778, -0.02994152047, -0.02919989367, 0.02,
G 0.03636363636, 0.1157407407, 0.5489278752, 0.535331384, -0.2/
H=(B-A)/M
Hmin=H/1024
Hmax=H*64
T(0)=A
Y(0)=Y0
T(0)=A
J=0
TJ=A
Br=B-0.00001*ABS(B)
10 IF (T(J).LT.B) THEN
IF ((T(J)+H).GT.Br) H=B-T(J)
TJ=T(J)
YJ=Y(J)
Y1=YJ
K1=H*F(TJ,Y1)
Y2=YJ+B2*K1
IF (Big.LT.ABS(Y2)) RETURN
K2=H*F(TJ+A2*H,Y2)
Y3=YJ+B3*K1+C3*K2
IF (Big.LT.ABS(Y3)) RETURN
K3=H*F(TJ+A3*H,Y3)
Y4=YJ+B4*K1+C4*K2+D4*K3
IF (Big.LT.ABS(Y4)) RETURN
K4=H*F(TJ+A4*H,Y4)
Y5=YJ+B5*K1+C5*K2+D5*K3+E5*K4
IF (Big.LT.ABS(Y5)) RETURN
K5=H*F(TJ+A5*H,Y5)
Y6=YJ+B6*K1+C6*K2+D6*K3+E6*K4+F6*K5
IF (Big.LT.ABS(Y6)) RETURN
K6=H*F(TJ+A6*H,Y6)
Err=ABS(R1*K1+R3*K3+R4*K4+R5*K5+R6*K6)
Ynew=YJ+N1*K1+N3*K3+N4*K4+N5*K5
Err=ABS(Err)
IF (Err.LT.Tol.OR.H.LE.2*Hmin) THEN
Y(J+1)=Ynew
IF ((TJ+H).GT.Br) THEN
T(J+1)=B
ELSE
T(J+1)=TJ+H
ENDIF
J=J+1
TJ=T(J)
ENDIF
IF (Err.EQ.0) THEN
S=0
ELSE
S=0.84*(Tol*H/Err)**0.25
ENDIF
IF (S.LT.0.75.AND.H.GT.(2*Hmin)) THEN
H=H/2
ENDIF
IF (S.GT.1.50.AND.(2*H).LT.Hmax) THEN
H=H*2
ENDIF
IF (Big.LT.ABS(Y(J)).OR.MaxM.EQ.J) RETURN
Mend=J
IF (B.GT.T(J)) THEN
M=J+1
ELSE
M=J
ENDIF
WRITE(9,1000) T(J),Y(J)
GOTO 10
ENDIF
RETURN
1000 FORMAT(5X,2F18.9)
END