SUBROUTINE RKQS(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C c Fifth-order Runge-Kutta step with monitoring of local truncation error c to ensure accuracy and adjust step size c c INPUTS: c Y(N) dependent variable vector at start of step c DYDX(N) derivative vector at start of step c X independent variable at start of step c HTRY stepsize to be attempted c EPS accuracy c YSCAL(N) vector against which the error is scaled c c OUTPUTS c Y(N) dependent variable vector at end of step c X independent variable at end of step c HDID step size used c HNEXT first guess at next step size c INTEGER n,NMAX DOUBLE PRECISION eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n) EXTERNAL derivs PARAMETER (NMAX=50) INTEGER i DOUBLE PRECISION errmax,h,xnew,yerr(NMAX),ytemp(NMAX),SAFETY, 1 PGROW,PSHRINK,ERRCON PARAMETER (SAFETY=0.9d0,PGROW=-0.2d0,PSHRINK=-0.25d0, 1 ERRCON=1.89d-4) C SET STEPSIZE TO THE INITIAL TRIAL VALUE H=HTRY c take a step 1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs) C EVALUATE ACCURACY ERRMAX=0.0 DO 11 i=1,N ERRMAX=MAX(ERRMAX,ABS(yerr(i)/yscal(I))) 11 CONTINUE C SCALE RELATIVE TO REQUIRED TOLERANCE ERRMAX=ERRMAX/EPS IF (ERRMAX.GT.1.d0) THEN C TRUNCATION ERROR TOO LARGE, REDUCE STEPSIZE H=SAFETY*H*(ERRMAX**PSHRINK) if (h.lt.0.1d0*h) then h=0.1d0*h endif xnew=x+h if (xnew.eq.x) STOP 'stepsize underflow in RKQS' GOTO 1 ELSE C STEP SUCCEEDED IF (ERRMAX.GT.ERRCON) THEN HNEXT=SAFETY*H*(ERRMAX**PGROW) ELSE HNEXT=5.d0*H ENDIF hdid=h x=x+h DO 12 I=1,N Y(I)=YTEMP(I) 12 CONTINUE RETURN endif END C SUBROUTINE RKCK(Y,DYDX,N,X,H,YOUT,YERR,DERIVS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C given values for N variables Y(N) and their derivatives DYDX(N) known c at X, use the fifth-order Cash-Karp Runge-Kutta method to advance the c solution over an interval H and return the incremented variables as YOUT(N) c Also return an estimate of the local truncation error in YOUT using the c embedded fourth-order method. The user supplies the subroutine c DERIVS(X,Y,YOUT) which returns derrivatives DYDX(N) at X c INTEGER n,NMAX DOUBLE PRECISION h,x,dydx(n),y(n),yerr(n),yout(n) EXTERNAL derivs PARAMETER (NMAX=50) DOUBLE PRECISION ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX), 1 ak6(NMAX),ytemp(NMAX),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43, 2 B51,B52,B53,B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3, 3 DC4,DC5,DC6 PARAMETER(A2=0.2d0,A3=0.3d0,A4=0.6d0,A5=1.d0,A6=0.875d0,B21=0.2d0, 1 B31=3.d0/40.d0,B32=9.d0/40.d0,B41=0.3d0,B42=-0.9d0,B43=1.2d0, 2 B51=-11.d0/54.d0,B52=2.5d0,B53=-70.d0/27.d0,B54=35.d0/27.d0, 3 B61=1631.d0/55296.d0,B62=175.d0/512.d0,B63=575.d0/13824.d0, 4 B64=44275.d0/110592.d0,B65=253.d0/4096.d0,C1=37.d0/378.d0, 5 C3=250.d0/621.d0,C4=125.d0/594.d0,C6=512.d0/1771.d0, 6 DC1=C1-2825.d0/27648.d0,DC3=C3-18575.d0/48384.d0, 7 DC4=C4-13525.d0/55296.d0,DC5=-277.d0/14336.d0,DC6=C6-0.25d0) C FIRST STEP DO 11 I=1,N ytemp(I)=Y(I)+B21*h*DYDX(I) 11 CONTINUE C SECOND STEP call derivs(X+A2*h,ytemp,ak2) DO 12 I=1,N ytemp(I)=Y(I)+H*(B31*dydx(i)+B32*ak2(i)) 12 CONTINUE C THIRD STEP call derivs(X+A3*h,ytemp,ak3) DO 13 I=1,N ytemp(I)=Y(I)+H*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i)) 13 CONTINUE C FOURTH STEP call derivs(X+A4*H,ytemp,ak4) do 14 i=1,n ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+ 1 B54*ak4(i)) 14 continue C FIFTH STEP call derivs(x+a5*h,ytemp,ak5) do 15 i=1,n ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)+ 1 B64*ak4(i)+B65*ak5(i)) 15 continue c SIXTH STEP call derivs(x+A6*h,ytemp,ak6) C ACCUMULATE INCREMENTS WITH PROPER WEIGHTS DO 16 I=1,N YOUT(I)=Y(I)+H*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+ 1 C6*ak6(i)) 16 continue c Estimate error as difference between fourth and fifth order methods. do 17 i=1,n yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i) 1 +DC6*ak6(i)) 17 continue RETURN END c