Next: Stiff Equations Up: The Bulirsch-Stoer Method Previous: Richardson Interpolation and Extrapolation

## Bulirsch-Stoer Step

```SUBROUTINE bsstep(y, dydx, x, htry, eps, yscal, hdid, hnext, derivs)

USE nrtype; USE nrutil, ONLY : arth, assert_eq, cumsum, iminloc, nrerror, &
outerdiff, outerprod, upper_triangle
USE nr, ONLY : mmid, pzextr
IMPLICIT NONE
REAL(sp), DIMENSION(:), INTENT(inout) :: y
REAL(sp), DIMENSION(:), INTENT(in) :: dydx, yscal
REAL(sp), INTENT(inout) :: x
REAL(sp), INTENT(in) :: htry, eps
REAL(sp), INTENT(out) :: hdid, hnext
INTERFACE
SUBROUTINE derivs(x, y, dydx)
USE nrtype
IMPLICIT NONE
real9sp), INTENT(in) :: x
REAL(sp), DIMENSION(:), INTENT(in) :: y
REAL(sp), DIMENSION(:), INTENT(out) :: dydx
END SUBROUTINE derivs
END INTERFACE
INTEGER(i4b), PARAMETER :: imax = 9, kmaxx=imax-1
REAL(sp), PARAMETER :: safe1=0.25_sp, safe2=0.7_sp, redmax=1.0e-5_sp, &
redmin=0.7_sp, tiny=1.0e-30_sp, scalmx=0.1_sp

! Bulirsch-Stoer step with monitoring of local truncation error to
! ensure accuracy and adjust stepsize. Input are athe dependent variable
! vector y and its derivative dydx at the starting value of the independent
! variable x. Also input are the stepsize that was actually accomplished,
! and hnext is the estimated next stepsize. derivs is the user-supplied
! subroutine that computes the right-hand-side derivatives. y, dydx, and
! yscal must all have the same length. Be sure to set htry on successive
! steps to the value of hnext returned from the previous step, as is the
! case if the routine is called by odeint.
! Parameters: kmaxx is the maximum row number used in the extrapolation;
! imax is the next row number; safe1 and safe2 are safety factors;
! redmax is the maximum factor used when a stepsize is reduced, redmin
! the minimum; tiny prevents division by zero; 1/scalmx is the maximum
! factor by which a stepsize can be increased.

INTEGER(i4b) :: k, km, ndum
INTEGER(i4b, DIMENSION(imax) :: nseq = (/ 2, 4, 6, 8, 10, 12, 14, 16, 18 /)
INTEGER(i4b), SAVE :: kopt, kmax
REAL(sp), DIMENSION(kmaxx, kmaxx), SAVE :: alf
REAL(sp), DIMENSION(kmaxx) :: err
REAL(sp), DIMENSION(imax), SAVE :: a
REAL(sp), SAVE :: epsold = -1.0_sp, xnew
REAL(sp) :: eps1, errmax, fact, h, red, scale, wrkmin, xest
REAL(sp), DIMENSION(SIZE(y)) :: yerr, ysav, yseq
LOGICAL(lgt) :: reduct
LOGICAL(lgt, SAVE :: first = .TRUE.

ndum = assert_eq(SIZE(y), SIZE(dydx), SIZE(yscal), 'bsstep')
IF (eps /= epsold) THEN              ! a new tolerance, reinitialize
hnext = -1.0e29_sp                ! "impossible" values
xnew=-1.0e29_sp
eps1=safe1*eps
a(:)=cumsum(nseq,1)
WHERE (upper_triangle(kmaxx, kmaxx)) alf=eps1** &
(outerdiff(a(2:), a(2:))/outerprod(arth( &
3.0_sp, 2.0_sp, kmaxx), (a(2:)-a(1)+1.0_sp)))
epsold=eps
DO kopt=2,kmaxx-1                 ! determine optimal row number for
IF (a(kopt+1) > a(kopt)*alf(kopt-1,kopt)) EXIT      ! convergence
END DO
kmax=kopt
END IF
h = htry
ysav(:) = y(:)                       ! save the starting values
IF (h /= hnext .OR. x /= xnew) THEN  ! a new stepsize or a new integration
first = .TRUE.                    ! re-establish the order window
kopt = kmax
END IF
reduct = .FALSE.
main_loop: DO
DO k = 1, kmax                    ! evaluate the sequence of modified
xnew = x+h                     ! midpoint integrations
IF (xnew == x) CALL nrerror('step size underflow in bsstep')
CALL mmid(ysav, dydx, x, h, nseq(k), yseq, derivs)
xest=(h/nseq(k))**2            !squared, since error series is even
CALL pzextr(k, xest, yseq, y, yerr) ! perform extrapolation
IF (k /= 1) THEN                  ! computer normalized error estimate
errmax=MAXVAL(ABS(yerr(:)/yscal(:)))
errmax=MAX(tiny, errmax)/eps   ! scale error relative to tolerance
km=k-1
err(km)=(errmax/safe1)**(1.0_sp/(2*km+1))
END IF
IF (k /= 1 .AND. (k >= kopt - 1 .OR. first)) THEN  ! in order window
IF (errmax < 1.0) EXIT main_loop      ! converged
IF(k == kmax .OR. k == kopt+1) THEN   ! check for possible step
red = safe2/err(km)                ! size reduction
EXIT
ELSE IF (k == kopt) THEN
IF (alf(kopt-1, kopt) < err(km)) THEN
red=1.0_sp/err(km)
EXIT
END IF
ELSE IF (kopt == kmax) THEN
IF(alf(km, kmax-1) < err(km)) THEN
red = alf(km, kmax-1) * safe2/err(km)
EXIT
END IF
ELSE IF (alf(km, kopt) < err(km)) THEN
red=alf(km,kopt=1) / err(km)
EXIT
END IF
END IF
END DO
red=MAX(MIN(red, redmin), redmax)   ! reduce step size by at least
h = h*red                           ! redmin and at most redmax
reduct = .TRUE.
END DO main_loop                       ! try again
x = xnew                               ! successful step taken
hdid=h
first=.FALSE.
kopt=1+iminloc(a(2:km+1)*MAX(err(1:km), scalmx))

! Compute optimal row for convergence and corresponding stepsize

scale=MAX(err(kopt=1), scalmx)
wrkmin=scale*a(kopt)
hnext=h/scale
IF(kopt >= k .AND. kopt /= kmax .AND. .NOT. reduct) THEN ! check for possible
fact = MAX(scale/alf(kopt-1, kopt), scalmx)           ! order increase
IF(a(kopt+1)*fact <=wrkmin) THEN                      ! but not if step
hnext = h/fact                                     ! size just reduced
kopt=kopt+1
END IF
END IF

END SUBROUTINE bsstep
```

Zdzislaw Meglicki
2001-02-26