subroutine bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,DERIVS) implicit double precision (a-h,o-z) c c Bulirsch-Stoer step with monitoring of local truncation error to c ensure accuracy and adjust stepsize. c INPUTS: c Y(NV) dependant variable vector of length NV c X value of independent variable c DYDX(NV) deriv of Y w.r.t. X at the starting value of X c HTRY first try at the step size in X c EPS required accuracy c YSCAL(NV) vector that scales EPS c c on OUTPUT: c Y(NV) values at end of step c X value of X at end of step c HDID size of step in X c HNEXT estimated next step size c c subroutine DERIVS supplies the derivatives of Y(NV) with respect to X c PARAMETER (NMAX=50,KMAXX=8,IMAX=KMAXX+1,SAFE1=0.25d0,SAFE2=0.7d0, 1 REDMAX=1.d-5,REDMIN=0.7d0,TINY=1.d-30,SCALMX=0.1d0) DIMENSION Y(NV),DYDX(NV),YSCAL(NV),YERR(NMAX), 1 a(IMAX),alf(KMAXX,KMAXX),err(KMAXX), 2 ysav(NMAX),yseq(NMAX),nseq(IMAX) LOGICAL first,reduct SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew EXTERNAL derivs DATA first/.true./,epsold/-1.d0/ DATA NSEQ/2,4,6,8,10,12,14,16,18/ c a new tolerance, so reinitialize if (eps.ne.epsold) then hnext=-1.d29 xnew=-1.d29 eps1=SAFE1*eps c compute work coefficients A(k) a(1)=nseq(1)+1 do 11 k=1,KMAXX a(k+1)=a(k)+nseq(k+1) 11 continue c compute alfa(k,q) do 13 iq=2,KMAXX do 12 k=1,iq-1 alf(k,iq)=eps1**((a(k+1)-a(iq+1))/ 1 ((a(iq+1)-a(1)+1.d0)*(2*k+1))) 12 continue 13 continue epsold=eps c Determine the optimal row number for convergence do 14 kopt=2,KMAXX-1 if (a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt)) goto 1 14 continue 1 kmax=kopt endif H=HTRY C save the starting value of x,y(nv) DO 15 I=1,NV YSAV(I)=Y(I) 15 CONTINUE c A new stepsize or a new integration; re-establish the order window if (h.ne.hnext.or.x.ne.xnew) then first=.true. kopt=kmax endif reduct=.false. C EVALUATE THE SEQUENCE OF MODIFIED MIDPOINT INTEGRATIONS 2 DO 17 k=1,kmax xnew=x+h if (xnew.eq.x) STOP 'step size underflow in BSSTEP' CALL MMID(ysav,dydx,nv,x,h,nseq(k),yseq,DERIVS) C SQUARED BELOW SINCE ERROR SERIES IS EVEN XEST=(H/NSEQ(k))**2 C PERFORM EXTRAPOLATION CALL RZEXTR(k,XEST,YSEQ,Y,YERR,NV) c CALL PZEXTR(k,XEST,YSEQ,Y,YERR,NV) C Compute normalized error estimate eps(k) if (k.ne.1) then ERRMAX=tiny DO 16 i=1,NV ERRMAX=DMAX1(ERRMAX,DABS(YERR(i)/YSCAL(i))) 16 CONTINUE C SCALE error RELATIVE TO TOLERANCE ERRMAX=ERRMAX/EPS km=k-1 err(km)=(errmax/SAFE1)**(1.d0/(2*km+1)) endif if (k.ne.1.and.(k.ge.kopt-1.or.first))then c converged IF (ERRMAX.LT.1.d0) goto 4 C check for possible stepsize reduction if (k.eq.kmax.or.k.eq.kopt+1) then red=SAFE2/err(km) goto 3 else if (k.eq.kopt) then if (alf(kopt-1,kopt).lt.err(km))then red=1.d0/err(km) goto 3 endif else if (kopt.eq.kmax) then if (alf(km,kmax-1).lt.err(km))then red=alf(km,kmax-1)*SAFE2/err(km) goto 3 endif else if (alf(km,kopt).lt.err(km))then red=alf(km,kopt-1)/err(km) goto 3 endif endif 17 continue c Reduce stepsize by at least REDMIN and at most REDMAX 3 red=dmin1(red,REDMIN) red=dmax1(red,REDMAX) h=h*red reduct=.true. c try again goto 2 c Succesful step taken! 4 x=xnew hdid=h first=.false. wrkmin=1.d35 c compute optimal row for convergence and corresponding stepsize do 18 kk=1,km fact=dmax1(err(kk),SCALMX) work=fact*a(kk+1) if (work.lt.wrkmin) then scale=fact wrkmin=work kopt=kk+1 endif 18 continue hnext=h/scale c Check for possible order increase but not if stepsize was just reduced if (kopt.ge.k.and.kopt.ne.kmax.and..not.reduct) then fact=dmax1(scale/alf(kopt-1,kopt),SCALMX) if (a(kopt+1)*fact.le.wrkmin)then hnext=h/fact kopt=kopt+1 endif endif return end C SUBROUTINE MMID(Y,DYDX,NVAR,XS,HTOT,NSTEP,YOUT,DERIVS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C c Modified midpoint step. Dependent variable vector Y(NVAR) and c its derivative vector DYDX(NVAR) are input at XS. Also input is c HTOT, the total step to be made, and NSTEP, the number of substeps c to be used. The output is returned as YOUT(NVAR), which need not be c a distinct array from Y(NVAR); if it is distinct, however, then Y(NVAR) c and DYDX(NVAR) are returned undamaged c PARAMETER(NMAX=50) external DERIVS DIMENSION Y(NVAR),DYDX(NVAR),YOUT(NVAR),YM(NMAX),YN(NMAX) C STEPSIZE THIS TRIP H=HTOT/NSTEP DO 11 I=1,NVAR YM(I)=Y(I) C FIRST STEP YN(I)=Y(I)+H*DYDX(I) 11 CONTINUE X=XS+H C WILL USE YOUT FOR TEMPORARY STORAGE OF DERIVATIVES CALL DERIVS(X,YN,YOUT) H2=2.D0*H C GENERAL STEP DO 13 N=2,NSTEP DO 12 I=1,NVAR SWAP=YM(I)+H2*YOUT(I) YM(I)=YN(I) YN(I)=SWAP 12 CONTINUE X=X+H CALL DERIVS(X,YN,YOUT) 13 CONTINUE C LAST STEP DO 14 I=1,NVAR YOUT(I)=0.5d0*(YM(I)+YN(I)+H*YOUT(I)) 14 CONTINUE RETURN END c subroutine pzextr(iest,xest,yest,yz,dy,nv) implicit double precision (a-h,o-z) c c Uses polynomial extrapolation to evaluate nv functions at x=0 by fitting a c polynomial to a sequence of estimates with progressively smaller values c x=xest, and corresponding function vectors yest(nv). This call is number c iest in the sequence of calls. Extrapolated function values are output c as yz(nv), and their estimated error is output as dy(nv). c parameter(IMAX=13,NMAX=50) dimension dy(nv),yest(nv),yz(nv),d(NMAX),qcol(NMAX,IMAX), 1 x(IMAX) SAVE qcol,x c Save current independent variable x(iest)=xest do 11 j=1,nv dy(j)=yest(j) yz(j)=yest(j) 11 continue c Store first estimate in first column if (iest.eq.1) then do 12 j=1,nv qcol(j,1)=yest(j) 12 continue else do 13 j=1,nv d(j)=yest(j) 13 continue do 15 k1=1,iest-1 delta=1.d0/(x(iest-k1)-xest) f1=xest*delta f2=x(iest-k1)*delta c Propagate tableau 1 diagonal more do 14 j=1,nv q=qcol(j,k1) qcol(j,k1)=dy(j) delta=d(j)-q dy(j)=f1*delta d(j)=f2*delta yz(j)=yz(j)+dy(j) 14 continue 15 continue do 16 j=1,nv qcol(j,iest)=dy(j) 16 continue endif return END c SUBROUTINE RZEXTR(IEST,XEST,YEST,YZ,DY,NV) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Uses diagonal rational function extrapolation to evaluate NV functions c at X=0 by fitting a diagonal rational functino to a sequence of estimates c with progressively smaller values X=XEST, and corresponding function c vectors YEST. This call is number IEST in the sequence of calls. c Extrapolated function values are output as YZ, and their estimated error is c output as DY. PARAMETER (IMAX=13,NMAX=50) C MAXIMUM EXPECTED VALUE OF NUSE IS NCOL, OF NV IS NMAX, OF IEST IS IMAX DIMENSION dy(nv),yest(nv),yz(nv),d(NMAX,IMAX),fx(IMAX),x(IMAX) SAVE d,x C SAVE CURRENT INDEPENDENT VARIABLE X(IEST)=XEST IF (IEST.EQ.1) THEN DO 11 J=1,NV YZ(J)=YEST(J) D(J,1)=YEST(J) DY(J)=YEST(J) 11 CONTINUE ELSE do 12 k=1,iest-1 fx(k+1)=x(iest-k)/xest 12 continue C EVALUATE NEXT DIAGONAL IN TABLEAU DO 14 J=1,NV YY=YEST(J) V=D(J,1) C=YY D(J,1)=YY DO 13 K=2,iest B1=FX(K)*V B=B1-C IF (B.NE.0.D0) THEN B=(C-V)/B DDY=C*B C=B1*B C CARE NEEDED TO AVOID DIVISION BY ZERO... ELSE DDY=V ENDIF if (k.ne.iest) v=d(j,k) D(J,K)=DDY YY=YY+DDY 13 CONTINUE DY(J)=DDY YZ(J)=YY 14 CONTINUE ENDIF RETURN END