************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE BSSTEP(NV,DYDX,HTRY,EPS,YSCAL,DERIVS, C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C in C $ X,Y, C ^^^ C inout C $ HDID,HNEXT) C ^^^^^^^^^^ C out C Bulirsch-Stoer stepper routine for numerical integration of C differential equations. C NV dependent variables Y and their derivatives DYDX are input at X C along with the attempted stepsize HTRY, the required accuracy EPS C and the scaling vector YSCAL. DERIVS is a routine for evaluating C the right-hand sides. C X and Y are replaced with their output values. C C HDID is the accomplished stepsize, HNEXT the next one. C C Based on ideas from NUMERICAL RECIPES. (The Second Edition C discusses more involved techniques. They are not used here C because tests showed that for this application the difference in C computing time was negligible; thus I opted for the shorter and C clearer code.) C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C named constants as array bounds INTEGER NMAX, IMAX PARAMETER (NMAX=10, IMAX=11) ************************************************************************ C input variables from the calling routine (parameter list) C number of dependent variables INTEGER NV C derivatives, attempted stepsize, accuracy, scaling vector, C independent and dependent variables DOUBLE PRECISION DYDX(NV), HTRY, EPS, YSCAL(NV), X, Y(NV) ************************************************************************ C output variables to the calling routine (parameter list) C accomplished stepsize, next stepsize DOUBLE PRECISION HDID, HNEXT ************************************************************************ C output to called routines (parameter list) C all called routines C modified midpoint method, polynomial extrapolation, right-hand C sides EXTERNAL MMM, POLEX, DERIVS C named constants C number of previous estimates extrapolation should use INTEGER NUSE PARAMETER (NUSE=7) C variables C sequence of number of steps INTEGER NSEQ(IMAX) C x, y, stepsize and derivatives for modified mitpoint, x for C polynomial extrapolation DOUBLE PRECISION XSAV, YSAV(NMAX), DYSAV(NMAX), H, XEST ************************************************************************ C output from called routines (parameter list) C output vector from polynomial extrapolation and its error DOUBLE PRECISION YSEQ(NMAX), YERR(NMAX) ************************************************************************ C variables in common blocks C to let routines higher up now if something goes wrong INTEGER MSTEP COMMON /MERROR/ MSTEP C save until the next call SAVE /MERROR/ ************************************************************************ C internally used quantities C named constants C ONE mustn't necessarily be `mathematically equivalent' to 1.0; C it's a parameter like the others---one could choose another value. C Has nothing to do with making it easier to code DOUBLE PRECISION C constants. DOUBLE PRECISION ONE, SHRINK, GROW PARAMETER (ONE=1.D0, SHRINK=.95D0, GROW=1.2D0) C variables C DO loops INTEGER I, J C for error checking DOUBLE PRECISION ERRMAX ************************************************************************ ************************************************************************ C initialisation of variables C internally used quantities DATA NSEQ /2,4,6,8,12,16,24,32,48,64,96/ ************************************************************************ ************************************************************************ C run-time initialisation MSTEP = 0 H = HTRY XSAV = X ************************************************************************ ************************************************************************ C calculations C save the starting values DO 100 I = 1, NV YSAV(I) = Y(I) DYSAV(I) = DYDX(I) 100 CONTINUE C GOTO 1 allows one to repeat the main part of the routine if the C step fails. 1 CONTINUE C evaluate sequence of modified midpoint integrations DO 200 I = 1, IMAX C call the routine for modified midpoint method CALL MMM(XSAV,NV,YSAV,DYSAV,H,NSEQ(I),DERIVS, YSEQ) C series error is even, so squared XEST = (H/NSEQ(I))**2 C extrapolation with polynomials CALL POLEX(NV,XEST,YSEQ,I,NUSE, Y,YERR) C monitor local truncation error ERRMAX = 0.0D0 DO 12 J = 1, NV ERRMAX = MAX(ERRMAX,ABS(YERR(J)/YSCAL(J))) 12 CONTINUE C scale accuracy relative to tolerance ERRMAX = ERRMAX/EPS C has the step converged? IF (ERRMAX .LT. ONE) THEN C set the return variables X = X + H HDID = H IF (I. EQ. NUSE) THEN HNEXT = H*SHRINK ELSE IF (I .EQ. NUSE - 1) THEN HNEXT = H*GROW ELSE HNEXT = (H*NSEQ(NUSE - 1))/NSEQ(I) ENDIF ************************************************************************ C normal return RETURN ENDIF ************************************************************************ C get this far if the step fails 200 CONTINUE C Here it doesn't matter, with IMAX=11 and NUSE=7; C if one changes these values, one should think about whether C or not the integer division in the exponent is really wanted. C (In practice, one will probably never even get here.) C Anyway, this reduces the step size for another try if the step C fails. H = 0.25D0*H/2**((IMAX - NUSE)/2) C If H is too small, then set a flag and get out. IF (X + H. EQ. X) THEN MSTEP = 1 ************************************************************************ C abnormal return RETURN END IF ************************************************************************ C repeat if one gets this far GOTO 1 ************************************************************************ ************************************************************************ END