请教用Fortran编写程序求解微分方程组
通过在论坛里学习了这么久，使我对于分岔有了一些初步的认识，并且可以在建立微分方程组后用matlab编写出求解方程及绘制具体图形。但是随着时间的推移，我越发觉得matlab的优点是后处理，前期计算非常缓慢，尤其是遇到循环，速度更是让人难以忍受（之前的一个敏感性分析，算了1000次，2个小时）。在做分岔图的时候也是如此。所以想采用科学计算能力更强的fortran进行前期计算，现在也正在学，但渴望提高快一些。我希望论坛里用fortran并且搞分岔稳定方面的朋友能够教我一下，我不是一个聪明人，但还算努力，苦于没有老师，研究这个方面的师兄们也早就不在了。
请大家提供一个duffing系统或者Rossler系统或者其他的任何一个微分方程组求解，保存数据的fortran语言例子，让我好好学习一下。
Chen系统的就可以，a，b，c参数任意可
dx/dt=a*(y-x)
dy/dt=(c-a)*x+c*y-x*z
dz/dt=x*y-b*z
如果搞清楚了，我还会同大家分享经验（我真的希望能够这样），谢谢！ 怎么突然想改弦更张了，fortran语言我不熟悉，等做出来了给我们也介绍下经验:victory:
回复楼主
你是想求解微分方程吗？Fortran有现成的龙格库塔积分程序，如果需要可以回复，我只是不知道你现在还需要不，呵呵回复 板凳 lixiaoju 的帖子
谢谢楼上，我找到了一本算法集，但是现在暂时不需要用fortran计算了。等再进行微分方程求解时，再向你请教。 这个就是用四阶龙格库塔方法积分微分方程，徐世良有本书，专门有这方面的编程程序，你可以找找啊。fortran用起来确实怪方便的。 谢谢楼上的建议！ lixiaoju 发表于 2010-7-2 09:43 static/image/common/back.gif你是想求解微分方程吗？Fortran有现成的龙格库塔积分程序，如果需要可以回复，我只是不知道你现在还需要不， ...
你好，我想要有关一阶微分方程组及求解的Fortran程序。你有么 笨与呆无关 发表于 2013-10-30 15:32
你好，我想要有关一阶微分方程组及求解的Fortran程序。你有么
一般fortran的函数库中就有各种龙哥库塔法的函数可以调用，实在没有的话可以参考下面的程序
经典四阶龙哥库塔法 SUBROUTINE runge4(N,X,Y,K,H,get_dydx)
C-----THIS SUBROUTINE ADVANCES THE SOLUTION TO THE SYSTEM OF N ODE'S
C-----DY(i)/DX = f(X,Y(1),Y(2),,.....Y(N)), i=1,...,N
C-----USING THE FOURTH ORDER RUNGE-KUTTA METHOD.
C-----N = Number of differential equations in the system. N<=10
C-----X = INITIAL VALUE OF THE INDEPENDENT VARIABLE; IT IS INCREASED
C----- BY THE SUBROUTINE.
C-----Y = INITIAL VALUE OF THE DEPENDENT VARIABLE(S); THE ROUTINE
C----- OVERWRITES THIS WITH THE NEW VALUE(S).
C-----H = time step
C-----K = number of steps to advance the equations
C-----get_dydx:subroutineget_dydx(X,Y,DYDX)
IMPLICIT REAL*8(A-H,O-Z)
EXTERNAL get_dydx
DIMENSION Y(10),YS(10),YSS(10),YSSS(10),T1(10),T2(10),
1 T3(10),T4(10)
C-----THE MAIN LOOP
C print*, X, Y(1), Y(2), Y(3), Y(4), Y(5)
DO 1 I=1,K
C-----TEMPORARY ARRAYS ARE NEEDED FOR THE FUNCTIONS TO SAVE THEM
C-----FOR THE FINAL CORRECTOR STEP.
C-----FIRST (HALF STEP) PREDICTOR
call get_dydx(X,Y,T1)
DO 2 J=1,N
YS(J) = Y(J) + .5 * H * T1(J)
2 CONTINUE
C-----SECOND STEP (HALF STEP CORRECTOR)
X = X + .5 * H
call get_dydx(X,YS,T2)
DO 3 J=1,N
YSS(J) = Y(J) + .5 * H * T2(J)
3 CONTINUE
C-----THIRD STEP (FULL STEP MIDPOINT PREDICTOR)
call get_dydx(X,YSS,T3)
DO 4 J=1,N
YSSS(J) = Y(J) + H * T3(J)
4 CONTINUE
C-----FINAL STEP (SIMPSON'S RULE CORRECTOR)
X = X + .5 * H
call get_dydx(X,YSSS,T4)
DO 5 J=1,N
Y(J) = Y(J) + (H/6.)*(T1(J) + 2.*(T2(J)+T3(J)) + T4(J))
5 CONTINUE
1 CONTINUE
RETURN
END
四五阶龙哥库塔法
SUBROUTINE RKF45(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK)
C
C FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD
C
C WRITTEN BY H.A.WATTS AND L.F.SHAMPINE
C SANDIA LABORATORIES
C ALBUQUERQUE,NEW MEXICO
C
C RKF45 IS PRIMARILY DESIGNED TO SOLVE NON-STIFF AND MILDLY STIFF
C DIFFERENTIAL EQUATIONS WHEN DERIVATIVE EVALUATIONS ARE INEXPENSIVE.
C RKF45 SHOULD GENERALLY NOT BE USED WHEN THE USER IS DEMANDING
C HIGH ACCURACY.
C
C ABSTRACT
C
C SUBROUTINERKF45INTEGRATES A SYSTEM OF NEQN FIRST ORDER
C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM
C DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN))
C WHERE THE Y(I) ARE GIVEN AT T .
C TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TOUT BUT IT
C CAN BE USED AS A ONE-STEP INTEGRATOR TO ADVANCE THE SOLUTION A
C SINGLE STEP IN THE DIRECTION OF TOUT.ON RETURN THE PARAMETERS IN
C THE CALL LIST ARE SET FOR CONTINUING THE INTEGRATION. THE USER HAS
C ONLY TO CALL RKF45 AGAIN (AND PERHAPS DEFINE A NEW VALUE FOR TOUT).
C ACTUALLY, RKF45 IS AN INTERFACING ROUTINE WHICH CALLS SUBROUTINE
C RKFS FOR THE SOLUTION.RKFS IN TURN CALLS SUBROUTINEFEHL WHICH
C COMPUTES AN APPROXIMATE SOLUTION OVER ONE STEP.
C
C RKF45USES THE RUNGE-KUTTA-FEHLBERG (4,5)METHOD DESCRIBED
C IN THE REFERENCE
C E.FEHLBERG , LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS WITH STEPSIZE
C CONTROL , NASA TR R-315
C
C THE PERFORMANCE OF RKF45 IS ILLUSTRATED IN THE REFERENCE
C L.F.SHAMPINE,H.A.WATTS,S.DAVENPORT, SOLVING NON-STIFF ORDINARY
C DIFFERENTIAL EQUATIONS-THE STATE OF THE ART ,
C SANDIA LABORATORIES REPORT SAND75-0182 ,
C TO APPEAR IN SIAM REVIEW.
C
C
C THE PARAMETERS REPRESENT-
C F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES YP(I)=DY(I)/DT
C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C Y(*) -- SOLUTION VECTOR AT T
C T -- INDEPENDENT VARIABLE
C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED
C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR LOCAL
C ERROR TEST. AT EACH STEP THE CODE REQUIRES THAT
C ABS(LOCAL ERROR) .LE. RELERR*ABS(Y) + ABSERR
C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS
C IFLAG -- INDICATOR FOR STATUS OF INTEGRATION
C WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO RKF45 WHICH IS
C NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED
C AT LEAST3+6*NEQN
C IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL TO
C RKF45 WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE
C DIMENSIONED AT LEAST5
C
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C IMPORTANT!!!!SOME OF THE VARIABLES IN RKF45 NEED TO BE
C DOUBLE PRECISION. THEY ARE:
C
C THE ARRAY "Y"
C T
C TOUT
C RELERR
C ABSERR
C THE ARRAY "WORK"
C
C
CFIRST CALL TO RKF45
C
C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAYS
C IN THE CALL LIST- Y(NEQN) , WORK(3+6*NEQN) , IWORK(5),
C DECLARE F IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE F(T,Y,YP) AND
C INITIALIZE THE FOLLOWING PARAMETERS-
C
C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED.(NEQN .GE. 1)
C Y(*) -- VECTOR OF INITIAL CONDITIONS
C T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE
C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED.
C T=TOUT IS ALLOWED ON THE FIRST CALL ONLY, IN WHICH CASE
C RKF45 RETURNS WITH IFLAG=2 IF CONTINUATION IS POSSIBLE.
C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES
C WHICH MUST BE NON-NEGATIVE. RELERR MUST BE A VARIABLE WHILE
C ABSERR MAY BE A CONSTANT. THE CODE SHOULD NORMALLY NOT BE
C USED WITH RELATIVE ERROR CONTROL SMALLER THAN ABOUT 1.E-8 .
C TO AVOID LIMITING PRECISION DIFFICULTIES THE CODE REQUIRES
C RELERR TO BE LARGER THAN AN INTERNALLY COMPUTED RELATIVE
C ERROR PARAMETER WHICH IS MACHINE DEPENDENT. IN PARTICULAR,
C PURE ABSOLUTE ERROR IS NOT PERMITTED. IF A SMALLER THAN
C ALLOWABLE VALUE OF RELERR IS ATTEMPTED, RKF45 INCREASES
C RELERR APPROPRIATELY AND RETURNS CONTROL TO THE USER BEFORE
C CONTINUING THE INTEGRATION.
C IFLAG -- +1,-1INDICATOR TO INITIALIZE THE CODE FOR EACH NEW
C PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG=-1
C ONLY WHEN ONE-STEP INTEGRATOR CONTROL IS ESSENTIAL. IN THIS
C CASE, RKF45 ATTEMPTS TO ADVANCE THE SOLUTION A SINGLE STEP
C IN THE DIRECTION OF TOUT EACH TIME IT IS CALLED. SINCE THIS
C MODE OF OPERATION RESULTS IN EXTRA COMPUTING OVERHEAD, IT
C SHOULD BE AVOIDED UNLESS NEEDED.
C
C
COUTPUT FROM RKF45
C
C Y(*) -- SOLUTION AT T
C T -- LAST POINT REACHED IN INTEGRATION.
C IFLAG = 2 -- INTEGRATION REACHED TOUT. INDICATES SUCCESSFUL RETUR
C AND IS THE NORMAL MODE FOR CONTINUING INTEGRATION.
C =-2 -- A SINGLE SUCCESSFUL STEP IN THE DIRECTION OF TOUT
C HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING
C INTEGRATION ONE STEP AT A TIME.
C = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE RELATIVE ERROR
C TOLERANCE WAS TOO SMALL. RELERR HAS BEEN INCREASED
C APPROPRIATELY FOR CONTINUING.
C = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN
C 3000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS
C IS APPROXIMATELY 500 STEPS.
C = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION
C VANISHED MAKING A PURE RELATIVE ERROR TEST
C IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE.
C USING THE ONE-STEP INTEGRATION MODE FOR ONE STEP
C IS A GOOD WAY TO PROCEED.
C = 6 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED
C ACCURACY COULD NOT BE ACHIEVED USING SMALLEST
C ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR
C TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE
C ATTEMPTED.
C = 7 -- IT IS LIKELY THAT RKF45 IS INEFFICIENT FOR SOLVING
C THIS PROBLEM. TOO MUCH OUTPUT IS RESTRICTING THE
C NATURAL STEPSIZE CHOICE. USE THE ONE-STEP INTEGRATOR
C MODE.
C = 8 -- INVALID INPUT PARAMETERS
C THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS
C SATISFIED - NEQN .LE. 0
C T=TOUTANDIFLAG .NE. +1 OR -1
C RELERR OR ABSERR .LT. 0.
C IFLAG .EQ. 0OR.LT. -2OR.GT. 8
C WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO INTEREST
C TO THE USER BUT NECESSARY FOR SUBSEQUENT CALLS.
C WORK(1),...,WORK(NEQN) CONTAIN THE FIRST DERIVATIVES
C OF THE SOLUTION VECTOR Y AT T. WORK(NEQN+1) CONTAINS
C THE STEPSIZE H TO BE ATTEMPTED ON THE NEXT STEP.
C IWORK(1) CONTAINS THE DERIVATIVE EVALUATION COUNTER.
C
C
CSUBSEQUENT CALLS TO RKF45
C
C SUBROUTINE RKF45 RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE
C THE INTEGRATION. IF THE INTEGRATION REACHED TOUT, THE USER NEED ONL
C DEFINE A NEW TOUT AND CALL RKF45 AGAIN. IN THE ONE-STEP INTEGRATOR
C MODE (IFLAG=-2) THE USER MUST KEEP IN MIND THAT EACH STEP TAKEN IS
C IN THE DIRECTION OF THE CURRENT TOUT. UPON REACHING TOUT (INDICATED
C BY CHANGING IFLAG TO 2),THE USER MUST THEN DEFINE A NEW TOUT AND
C RESET IFLAG TO -2 TO CONTINUE IN THE ONE-STEP INTEGRATOR MODE.
C
C IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS TO
C CONTINUE (IFLAG=3,4 CASES), HE JUST CALLS RKF45 AGAIN. WITH IFLAG=3
C THE RELERR PARAMETER HAS BEEN ADJUSTED APPROPRIATELY FOR CONTINUING
C THE INTEGRATION. IN THE CASE OF IFLAG=4 THE FUNCTION COUNTER WILL
C BE RESET TO 0 AND ANOTHER 3000 FUNCTION EVALUATIONS ARE ALLOWED.
C
C HOWEVER,IN THE CASE IFLAG=5, THE USER MUST FIRST ALTER THE ERROR
C CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE INTEGRATION CAN
C PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED.
C
C ALSO,IN THE CASE IFLAG=6, IT IS NECESSARY FOR THE USER TO RESET
C IFLAG TO 2 (OR -2 WHEN THE ONE-STEP INTEGRATION MODE IS BEING USED)
C AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH BEFORE THE
C INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, EXECUTION WILL
C BE TERMINATED. THE OCCURRENCE OF IFLAG=6 INDICATES A TROUBLE SPOT
C (SOLUTION IS CHANGING RAPIDLY,SINGULARITY MAY BE PRESENT) AND IT
C OFTEN IS INADVISABLE TO CONTINUE.
C
C IF IFLAG=7 IS ENCOUNTERED, THE USER SHOULD USE THE ONE-STEP
C INTEGRATION MODE WITH THE STEPSIZE DETERMINED BY THE CODE OR
C CONSIDER SWITCHING TO THE ADAMS CODES DE/STEP,INTRP. IF THE USER
C INSISTS UPON CONTINUING THE INTEGRATION WITH RKF45, HE MUST RESET
C IFLAG TO 2 BEFORE CALLING RKF45 AGAIN. OTHERWISE,EXECUTION WILL BE
C TERMINATED.
C
C IF IFLAG=8 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS
C THE INVALID INPUT PARAMETERS ARE CORRECTED.
C
C IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN INFORMATION
C REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, WORK AND IWORK
C SHOULD NOT BE ALTERED.
C
C
INTEGER NEQN,IFLAG,IWORK(5)
DOUBLE PRECISION Y(NEQN),T,TOUT,RELERR,ABSERR,WORK(1)
C
EXTERNAL F
C
INTEGER K1,K2,K3,K4,K5,K6,K1M
C
C
C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY
C
K1M=NEQN+1
K1=K1M+1
K2=K1+NEQN
K3=K2+NEQN
K4=K3+NEQN
K5=K4+NEQN
K6=K5+NEQN
C
C THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG
C CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE
C ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER,
C HE MUST USE RKFS DIRECTLY.
C
CALL RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK(1),WORK(K1M),
1 WORK(K1),WORK(K2),WORK(K3),WORK(K4),WORK(K5),WORK(K6),
2 WORK(K6+1),IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5))
C
RETURN
END
SUBROUTINE RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,YP,H,F1,F2,F3,
1 F4,F5,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,KFLAG)
C
C FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD
C
C
C RKFS INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL
C EQUATIONS AS DESCRIBED IN THE COMMENTS FOR RKF45 .
C THE ARRAYS YP,F1,F2,F3,F4,AND F5 (OF DIMENSION AT LEAST NEQN) AND
C THE VARIABLES H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,AND KFLAG ARE USED
C INTERNALLY BY THE CODE AND APPEAR IN THE CALL LIST TO ELIMINATE
C LOCAL RETENTION OF VARIABLES BETWEEN CALLS. ACCORDINGLY, THEY
C SHOULD NOT BE ALTERED. ITEMS OF POSSIBLE INTEREST ARE
C YP - DERIVATIVE OF SOLUTION VECTOR AT T
C H- AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP
C NFE- COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION EVALUATIONS
C
C
LOGICAL HFAILD,OUTPUT
C
INTEGERNEQN,IFLAG,NFE,KOP,INIT,JFLAG,KFLAG
DOUBLE PRECISIONY(NEQN),T,TOUT,RELERR,ABSERR,H,YP(NEQN),
1F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN),SAVRE,
2SAVAE
C
EXTERNAL F
C
DOUBLE PRECISIONA,AE,DT,EE,EEOET,ESTTOL,ET,HMIN,REMIN,RER,S,
1SCALE,TOL,TOLN,U26,EPSP1,EPS,YPK
C
INTEGERK,MAXNFE,MFLAG
C
DOUBLE PRECISIONDABS,DMAX1,DMIN1,DSIGN
C
CREMIN IS THE MINIMUM ACCEPTABLE VALUE OF RELERR.ATTEMPTS
CTO OBTAIN HIGHER ACCURACY WITH THIS SUBROUTINE ARE USUALLY
CVERY EXPENSIVE AND OFTEN UNSUCCESSFUL.
C
DATA REMIN/1.D-12/
C
C
C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER
C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE.
C AS SET, THIS CORRESPONDS TO ABOUT 500 STEPS.
C
DATA MAXNFE/3000/
C
C
C CHECK INPUT PARAMETERS
C
C
IF (NEQN .LT. 1) GO TO 10
IF ((RELERR .LT. 0.0D0).OR.(ABSERR .LT. 0.0D0)) GO TO 10
MFLAG=IABS(IFLAG)
IF ((MFLAG .GE. 1).AND.(MFLAG .LE. 8)) GO TO 20
C
C INVALID INPUT
10 IFLAG=8
RETURN
C
C IS THIS THE FIRST CALL
20 IF (MFLAG .EQ. 1) GO TO 50
C
C CHECK CONTINUATION POSSIBILITIES
C
IF ((T .EQ. TOUT) .AND. (KFLAG .NE. 3)) GO TO 10
IF (MFLAG .NE. 2) GO TO 25
C
C IFLAG = +2 OR -2
IF (KFLAG .EQ. 3) GO TO 45
IF (INIT .EQ. 0) GO TO 45
IF (KFLAG .EQ. 4) GO TO 40
IF ((KFLAG .EQ. 5).AND.(ABSERR .EQ. 0.0D0)) GO TO 30
IF ((KFLAG .EQ. 6).AND.(RELERR .LE. SAVRE).AND.
1 (ABSERR .LE. SAVAE)) GO TO 30
GO TO 50
C
C IFLAG = 3,4,5,6,7 OR 8
25 IF (IFLAG .EQ. 3) GO TO 45
IF (IFLAG .EQ. 4) GO TO 40
IF ((IFLAG .EQ. 5) .AND. (ABSERR .GT. 0.0D0)) GO TO 45
C
C INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO
C THE INSTRUCTIONS PERTAINING TO IFLAG=5,6,7 OR 8
30 STOP
C
C RESET FUNCTION EVALUATION COUNTER
40 NFE=0
IF (MFLAG .EQ. 2) GO TO 50
C
C RESET FLAG VALUE FROM PREVIOUS CALL
45 IFLAG=JFLAG
IF (KFLAG .EQ. 3) MFLAG=IABS(IFLAG)
C
C SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT
C INPUT CHECKING
50 JFLAG=IFLAG
KFLAG=0
C
C SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS
SAVRE=RELERR
SAVAE=ABSERR
C
CCOMPUTE MACHINE EPSILON
C
EPS = 1.0D0
53 EPS = EPS/2.0D0
EPSP1 = EPS + 1.0D0
IF (EPSP1 .GT. 1.0D0) GO TO 53
C
C RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS
C 2*EPS+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING
C FROM IMPOSSIBLE ACCURACY REQUESTS
C
RER=2.0D0*EPS+REMIN
IF (RELERR .GE. RER) GO TO 55
C
C RELATIVE ERROR TOLERANCE TOO SMALL
RELERR=RER
IFLAG=3
KFLAG=3
RETURN
C
55 U26=26.0D0*EPS
C
DT=TOUT-T
C
IF (MFLAG .EQ. 1) GO TO 60
IF (INIT .EQ. 0) GO TO 65
GO TO 80
C
C INITIALIZATION --
C SET INITIALIZATION COMPLETION INDICATOR,INIT
C SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP
C EVALUATE INITIAL DERIVATIVES
C SET COUNTER FOR FUNCTION EVALUATIONS,NFE
C EVALUATE INITIAL DERIVATIVES
C SET COUNTER FOR FUNCTION EVALUATIONS,NFE
C ESTIMATE STARTING STEPSIZE
C
60 INIT=0
KOP=0
C
A=T
CALL F(A,Y,YP)
NFE=1
IF (T .NE. TOUT) GO TO 65
IFLAG=2
RETURN
C
C
65 INIT=1
H=DABS(DT)
TOLN=0.
DO 70 K=1,NEQN
TOL=RELERR*DABS(Y(K))+ABSERR
IF (TOL .LE. 0.) GO TO 70
TOLN=TOL
YPK=DABS(YP(K))
IF (YPK*H**5 .GT. TOL) H=(TOL/YPK)**0.2D0
70 CONTINUE
IF (TOLN .LE. 0.0D0) H=0.0D0
H=DMAX1(H,U26*DMAX1(DABS(T),DABS(DT)))
JFLAG=ISIGN(2,IFLAG)
C
C
C SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT
C
80 H=DSIGN(H,DT)
C
C TEST TO SEE IF RKF45 IS BEING SEVERELY IMPACTED BY TOO MANY
C OUTPUT POINTS
C
IF (DABS(H) .GE. 2.0D0*DABS(DT)) KOP=KOP+1
IF (KOP .NE. 100) GO TO 85
C
C UNNECESSARY FREQUENCY OF OUTPUT
KOP=0
IFLAG=7
RETURN
C
85 IF (DABS(DT) .GT. U26*DABS(T)) GO TO 95
C
C IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND RETURN
C
DO 90 K=1,NEQN
90 Y(K)=Y(K)+DT*YP(K)
A=TOUT
CALL F(A,Y,YP)
NFE=NFE+1
GO TO 300
C
C
C INITIALIZE OUTPUT POINT INDICATOR
C
95 OUTPUT= .FALSE.
C
C TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION,
C SCALE THE ERROR TOLERANCES
C
SCALE=2.0D0/RELERR
AE=SCALE*ABSERR
C
C
C STEP BY STEP INTEGRATION
C
100 HFAILD= .FALSE.
C
C SET SMALLEST ALLOWABLE STEPSIZE
C
HMIN=U26*DABS(T)
C
C ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT.
C LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE AND
C THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE.
C
DT=TOUT-T
IF (DABS(DT) .GE. 2.0D0*DABS(H)) GO TO 200
IF (DABS(DT) .GT. DABS(H)) GO TO 150
C
C THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE
C OUTPUT POINT
C
OUTPUT= .TRUE.
H=DT
GO TO 200
C
150 H=0.5D0*DT
C
C
C
C CORE INTEGRATOR FOR TAKING A SINGLE STEP
C
C THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW IN
C COMPUTING THE ERROR TOLERANCE FUNCTION ET.
C TO AVOID PROBLEMS WITH ZERO CROSSINGS,RELATIVE ERROR IS MEASURED
C USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION AT THE
C BEGINNING AND END OF A STEP.
C THE ERROR ESTIMATE FORMULA HAS BEEN GROUPED TO CONTROL LOSS OF
C SIGNIFICANCE.
C TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED
C TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T.
C PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO
C SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE
C CHATTERING ON PROBLEMS HAVING DISCONTINUITIES.
C TO PREVENT UNNECESSARY FAILURES, THE CODE USES 9/10 THE STEPSIZE
C IT ESTIMATES WILL SUCCEED.
C AFTER A STEP FAILURE, THE STEPSIZE IS NOT ALLOWED TO INCREASE FOR
C THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE EFFICIENT ON
C PROBLEMS HAVING DISCONTINUITIES AND MORE EFFECTIVE IN GENERAL
C SINCE LOCAL EXTRAPOLATION IS BEING USED AND EXTRA CAUTION SEEMS
C WARRANTED.
C
C
C TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS.
C IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H
C
200 IF (NFE .LE. MAXNFE) GO TO 220
C
C TOO MUCH WORK
IFLAG=4
KFLAG=4
RETURN
C
C ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H
C
220 CALL FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,F1)
NFE=NFE+5
C
C COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR ESTIMATES
C AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE ERROR IS
C MEASURED WITH RESPECT TO THE AVERAGE OF THE MAGNITUDES OF THE
C SOLUTION AT THE BEGINNING AND END OF THE STEP.
C
EEOET=0.0D0
DO 250 K=1,NEQN
ET=DABS(Y(K))+DABS(F1(K))+AE
IF (ET .GT. 0.0D0) GO TO 240
C
C INAPPROPRIATE ERROR TOLERANCE
IFLAG=5
RETURN
C
240 EE=DABS((-2090.0D0*YP(K)+(21970.0D0*F3(K)-15048.0D0*F4(K)))+
1 (22528.0D0*F2(K)-27360.0D0*F5(K)))
250 EEOET=DMAX1(EEOET,EE/ET)
C
ESTTOL=DABS(H)*EEOET*SCALE/752400.0D0
C
IF (ESTTOL .LE. 1.0D0) GO TO 260
C
C
C UNSUCCESSFUL STEP
C REDUCE THE STEPSIZE , TRY AGAIN
C THE DECREASE IS LIMITED TO A FACTOR OF 1/10
C
HFAILD= .TRUE.
OUTPUT= .FALSE.
S=0.1D0
IF (ESTTOL .LT. 59049.0D0) S=0.9D0/ESTTOL**0.2D0
H=S*H
IF (DABS(H) .GT. HMIN) GO TO 200
C
C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE
IFLAG=6
KFLAG=6
RETURN
C
C
C SUCCESSFUL STEP
C STORE SOLUTION AT T+H
C AND EVALUATE DERIVATIVES THERE
C
260 T=T+H
DO 270 K=1,NEQN
270 Y(K)=F1(K)
A=T
CALL F(A,Y,YP)
NFE=NFE+1
C
C
C CHOOSE NEXT STEPSIZE
C THE INCREASE IS LIMITED TO A FACTOR OF 5
C IF STEP FAILURE HAS JUST OCCURRED, NEXT
C STEPSIZE IS NOT ALLOWED TO INCREASE
C
S=5.0D0
IF (ESTTOL .GT. 1.889568D-4) S=0.9D0/ESTTOL**0.2D0
IF (HFAILD) S=DMIN1(S,1.0D0)
H=DSIGN(DMAX1(S*DABS(H),HMIN),H)
C
C END OF CORE INTEGRATOR
C
C
C SHOULD WE TAKE ANOTHER STEP
C
IF (OUTPUT) GO TO 300
IF (IFLAG .GT. 0) GO TO 100
C
C
C INTEGRATION SUCCESSFULLY COMPLETED
C
C ONE-STEP MODE
IFLAG=-2
RETURN
C
C INTERVAL MODE
300 T=TOUT
IFLAG=2
RETURN
C
END
SUBROUTINE FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,S)
C
C FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD
C
C FEHL INTEGRATES A SYSTEM OF NEQN FIRST ORDER
C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM
C DY(I)/DT=F(T,Y(1),---,Y(NEQN))
C WHERE THE INITIAL VALUES Y(I) AND THE INITIAL DERIVATIVES
C YP(I) ARE SPECIFIED AT THE STARTING POINT T. FEHL ADVANCES
C THE SOLUTION OVER THE FIXED STEP H AND RETURNS
C THE FIFTH ORDER (SIXTH ORDER ACCURATE LOCALLY) SOLUTION
C APPROXIMATION AT T+H IN ARRAY S(I).
C F1,---,F5 ARE ARRAYS OF DIMENSION NEQN WHICH ARE NEEDED
C FOR INTERNAL STORAGE.
C THE FORMULAS HAVE BEEN GROUPED TO CONTROL LOSS OF SIGNIFICANCE.
C FEHL SHOULD BE CALLED WITH AN H NOT SMALLER THAN 13 UNITS OF
C ROUNDOFF IN T SO THAT THE VARIOUS INDEPENDENT ARGUMENTS CAN BE
C DISTINGUISHED.
C
C
external F
INTEGERNEQN
DOUBLE PRECISIONY(NEQN),T,H,YP(NEQN),F1(NEQN),F2(NEQN),
1F3(NEQN),F4(NEQN),F5(NEQN),S(NEQN)
C
DOUBLE PRECISIONCH
INTEGERK
C
CH=H/4.0D0
DO 221 K=1,NEQN
221 F5(K)=Y(K)+CH*YP(K)
CALL F(T+CH,F5,F1)
C
CH=3.0D0*H/32.0D0
DO 222 K=1,NEQN
222 F5(K)=Y(K)+CH*(YP(K)+3.0D0*F1(K))
CALL F(T+3.0D0*H/8.0D0,F5,F2)
C
CH=H/2197.0D0
DO 223 K=1,NEQN
223 F5(K)=Y(K)+CH*(1932.0D0*YP(K)+(7296.0D0*F2(K)-7200.0D0*F1(K)))
CALL F(T+12.0D0*H/13.0D0,F5,F3)
C
CH=H/4104.0D0
DO 224 K=1,NEQN
224 F5(K)=Y(K)+CH*((8341.0D0*YP(K)-845.0D0*F3(K))+
1 (29440.0D0*F2(K)-32832.0D0*F1(K)))
CALL F(T+H,F5,F4)
C
CH=H/20520.0D0
DO 225 K=1,NEQN
225 F1(K)=Y(K)+CH*((-6080.0D0*YP(K)+(9295.0D0*F3(K)-
1 5643.0D0*F4(K)))+(41040.0D0*F1(K)-28352.0D0*F2(K)))
CALL F(T+H/2.0D0,F1,F5)
C
C COMPUTE APPROXIMATE SOLUTION AT T+H
C
CH=H/7618050.0D0
DO 230 K=1,NEQN
230 S(K)=Y(K)+CH*((902880.0D0*YP(K)+(3855735.0D0*F3(K)-
1 1371249.0D0*F4(K)))+(3953664.0D0*F2(K)+
2 277020.0D0*F5(K)))
C
RETURN
END 四五阶龙哥库塔法使用范例
CSAMPLE PROGRAM FOR RKF45
C
C
EXTERNAL ORBIT
DOUBLE PRECISION T,Y(4),TOUT,RELERR,ABSERR
DOUBLE PRECISION TFINAL,TPRINT,ECC,ALFA,ALFASQ,WORK(27)
INTEGER IWORK(5), IFLAG, NEQN
DOUBLE PRECISION DSQRT
COMMON ALFASQ
ECC = 0.25
ALFA = 3.141592653589/4.0
ALFASQ = ALFA*ALFA
NEQN = 4
T = 0.0
Y(1) = 1.0 - ECC
Y(2) = 0.0
Y(3) = 0.0
Y(4) = ALFA*DSQRT( (1.0 + ECC)/(1.0 - ECC) )
RELERR = 1.0D-9
ABSERR = 0.0
TFINAL = 12.0
TPRINT = 0.5
IFLAG = 1
TOUT = T
10 CALL RKF45(ORBIT,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK)
WRITE(1,11) T, Y(1), Y(2)
GO TO (80,20,30,40,50,60,70,80), IFLAG
20 TOUT = T + TPRINT
IF (T .LT. TFINAL) GO TO 10
STOP
30 WRITE(1,31) RELERR,ABSERR
GO TO 10
40 WRITE(1,41)
GO TO 10
50 ABSERR = 1.0D-9
WRITE(1,31) RELERR,ABSERR
GO TO 10
60 RELERR = 10.0*RELERR
WRITE(1,31) RELERR,ABSERR
IFLAG = 2
GO TO 10
70 WRITE(1,71)
IFLAG = 2
GO TO 10
80 WRITE(1,81)
STOP
C
11 FORMAT(F5.1, 2F15.9)
31 FORMAT(17H TOLERANCES RESET, 2E12.3)
41 FORMAT(11H MANY STEPS)
71 FORMAT(12H MUCH OUTPUT)
81 FORMAT(14H IMPROPER CALL)
END
SUBROUTINE ORBIT (T, Y, YP)
DOUBLE PRECISION T, Y(4), YP(4), R, ALFASQ
DOUBLE PRECISION DSQRT
COMMON ALFASQ
R = Y(1)*Y(1) + Y(2)*Y(2)
R = R*DSQRT(R)/ALFASQ
YP(1) = Y(3)
YP(2) = Y(4)
YP(3) = -Y(1)/R
YP(4) = -Y(2)/R
RETURN
END
页:
[1]