next up previous index
Next: Bulirsch-Stoer Step Up: The Bulirsch-Stoer Method Previous: The Discussion

Richardson Interpolation and Extrapolation

SUBROUTINE polint(xa, ya, x, y, dy)

  USE nrtype; USE nrutil, ONLY : assert_eq, iminloc, nrerror
  IMPLICIT NONE
  REAL(sp), DIMENSION(:), INTENT(in) :: xa, ya
  REAL(sp), INTENT(in) :: x
  REAL(sp), INTENT(out) :: y, dy

  ! Given arrays xa and ya of length N, and given a value x, this routine
  ! returns a value y, and an error estimate dy. If P(x) is the polynomial
  ! of degree N - 1 such that P(x a_i) = y a_i, i = 1, ..., N, then the
  ! returned value y = P(x).

  INTEGER(i4b) :: m, n, ns
  REAL(sp), DIMENSION(SIZE(xa)) :: c, d, den, ho

  n = assert_eq(SIZE(xa), SIZE(ya), 'polint')
  c = ya                          ! Initialize the tableau of c's and d's
  d = ya
  ho = xa - x
  ns = iminloc(ABS(x - xa))       ! Find index ns of closest table entry
  y = ya(ns)                      ! Initial approximation to y.
  ns = ns - 1
  DO m = 1, n - 1                        ! For each column of the tableau
     den(1:n-m) = ho(1:n-m) - ho(1+m:n)  ! we loop over c's and d's and
     IF (ANY(den(1:n-m) == 0.0)) &       ! update them
          CALL nrerror('polint: calculation failure')

     ! This error can occur only if two input xa's are (to within roundoff)
     ! identical.

     den(1:n - m) = (c(2:n-m+1) - d(1:n-m))/den(1:n-m)
     d(1:n-m) = ho(1+m:n) * den(1:n-m)   ! Here c's and d's get updated
     c(1:n-m) = ho(1:n-m) * den(1:n-m)
     IF (2 * ns < n-m) THEN       ! After each column in the tableau is
        dy=c(ns+1)                ! completed decide, which correction
     ELSE                         ! c or d we add to y. We take the
        dy=d(ns)                  ! straightest line through the tableau
        ns=ns-1                   ! to its apex. The partial approximations
     END IF                       ! are thus centred on x. The last dy
     y = y+dy                     ! is the measure of error.
  END DO

END SUBROUTINE polint

SUBROUTINE ratint(xa, ya, x, y, dy)

  USE nrtype; USE nrutil, ONLY : assert_eq, iminloc, nrerror
  IMPLICIT NONE
  REAL(sp), DIMENSION(:), INTENT(in) :: xa, ya
  REAL(sp), INTENT(in) :: x
  REAL(sp), INTENT(out) :: y, dy

  ! Given arrays xa and ya of length N, and given a value of x, this routine
  ! returns a value of y and an accuracy estimate dy. The value returned is 
  ! that of the diagonal rational function, evaluated at x, that passes
  ! through the N points (xa_i, ya_i), i = 1... N.

  INTEGER(i4b) :: m, n, ns
  REAL(sp), DIMENSION(SIZE(xa)) :: c, d, dd, h, t
  REAL(sp), PARAMETER :: tiny=1.0e-25_sp

  n = assert_eq(SIZE(xa), SIZE(ya), 'ratint')
  h = xa - x
  ns = iminloc(ABS(h))
  y = ya(ns)
  IF (x == xa(ns)) THEN
     dy = 0.0
     RETURN
  END IF
  c = ya
  d = ya + tiny               ! The tiny is needed to prevent 0/0
  ns = ns - 1
  DO m=1, n-1
     t(1:n-m) = (xa(1:n-m)-x) * d(1:n-m)/h(1+m:n)   ! h will never be 0
     dd(1:n-m) = t(1:n-m) - c(2:n-m+1)
     IF (ANY(dd(1:n-m) == 0.0)) &                   ! interpolating function
          CALL nrerror('failure in ratint')         ! has a pole here
     dd(1:n-m) = (c(2:n-m+1) - d(1:n-m))/dd(1:n-m)
     d(1:n-m) = c(2:n-m+1) * dd(1:n-m)
     c(1:n-m) = t(1:n-m) * dd(1:n-m)
     IF(2*ns < n-m) THEN
        dy = c(ns+1)
     ELSE
        dy = d(ns)
        ns = ns-1
     END IF
     y=y+dy
  END DO

END SUBROUTINE ratint

SUBROUTINE pzextr(iest, xest, yest, yz, dy)

  USE nrtype; USE nrutil, ONLY : assert_eq, nrerror
  IMPLICIT NONE
  INTEGER(i4b), ikntent(in) :: iest
  REAL(sp), INTENT(in) :: xest
  REAL(sp), DIMENSION(:), INTENT(in) :: yest
  REAL(sp), DIMENSION(:), INTENT(out) :: yz, dy

  ! Use polynomial extrapolation to evaluate N functions at x = 0 by fitting
  ! a polynomial to a sequence of estimates with progressively smaller
  ! values x = xest, and corresponding function vectors yest. This call 
  ! is number iest in the sequence of calls. Extrapolated function values
  ! are output as yz, and their estimated error is output as dy. yest, yz,
  ! and dy are arrays of length N.

  INTEGER(i4b), PARAMETER :: iest_max = 16
  INTEGER(i4b) :: j, nv
  INTEGER(i4b), SAVE :: nvold = -1
  REAL(sp) :: delta, f1, f2
  REAL(sp), DIMENSION(SIZE(yz)) :: d, tmp, q
  REAL(sp), DIMENSION(iest_max), SAVE :: x
  REAL(sp), DIMENSION(:,:), ALLOCATABLE, SAVE :: qcol

  nv = assert_eq(SIZE(yz), SIZE(yest), SIZE(dy), 'pzextr')
  IF (iest > iest_max) &
       CALL nrerror('pzextr: probable misuse, too much extrapolation')
  IF(nv /= nvold) THEN      ! set up internal storage
     IF(ALLOCATED(qcol)) DEALLOCATE(qcol)
     ALLOCATE(qcol(nv, iest_max))
     nvold=nv
  END IF
  x(iest) = xest            ! save current independent variable
  dy(:) = yest(:)
  yz(:) = yest(:)
  IF (iest == 1) THEN       ! store first estimate in first column
     qcol(:,1)=yest(:)      
  ELSE
     d(:) = yest(:)
     DO j=1, iest-1
        delta=1.0_sp/(x(iest-j)-xest)
        f1=xest*delta
        f2=x(iest-j)*delta
        q(:) = qcol(:,j)    ! propagate tableau 1 diagonal more
        qcol(:,j)=dy(:)
        tmp(:)=d(:)-q(:)
        dy(:)=f1*tmp(:)
        d(:)=f2*tmp(:)
        yz(:)=yz(:)+dy(:)
     END DO
     qcol(:,iest)=dy(:)
  END IF

END SUBROUTINE pzextr

SUBROUTINE rzextr(iest, xest, yest, yz, dy)

  USE nrtype; USE nrutil, ONLY : assert_eq, nrerror
  IMPLICIT NONE
  INTEGER(i4b), INTENT(in) :: iest
  REAL(sp), INTENT(in) :: xest
  REAL(sp), DIMENSION(:), INTENT(in) :: yest
  REAL(sp), DIMENSION(:), INTENT(out) :: yz, dy

  ! Exact substitute for pzextr, but uses diagonal rational function 
  ! extrapolation instead of polynomial extrapolation.

  INTEGER(i4b), PARAMETER :: iest_max = 16
  INTEGER(i4b) :: k, nv
  INTEGER(i4b), SAVE :: nvold = -1
  REAL(sp), DIMENSION(SIZE(yz)) :: yy, v, c, b, b1, ddy
  REAL(sp), DIMENSION(:,:), ALLOCATABLE, SAVE :: d
  REAL(sp), DIMENSION(iest_max), SAVE :: fx, x

  nv = assert_eq(SIZE(yz), SIZE(dy), SIZE(yest), 'rzextr')
  IF(iest > iest_max) &
       CALL nrerror('rzextr: probable misuse, too much extrapolation')
  IF (nv /= nvold) THEN
     IF(ALLOCATED(d)) DEALLOCATE(d)
     ALLOCATE(d(nv, iest_max))
     nvold=nv
  END IF
  x(iest)=xest                 ! save current independent variable
  IF(iest == 1) THEN
     yz=yest
     d(:,1)=yest
     dy=yest
  ELSE
     fx(2:iest)=x(iest-1:1:-1)/xest
     yy=yest                   ! evaluate next diagonal in tableau
     v=d(1:nv, 1)
     c=yy
     d(1:nv, 1)=yy
     DO k=2, iest
        b1=fx(k)*v
        b=b1-c
        WHERE(b/=0.0)
           b=(c-v)/b
           ddy=c*b
           c=b1*b
        ELSEWHERE              ! care needed to avoid division by 0
           ddy=v
        END WHERE
        IF (k/=iest) v=d(1:nv, k)
        d(1:nv, k) = ddy
        yy=yy+ddy
     END DO
     dy=ddy
     yz=yy
  END IF

END SUBROUTINE rzextr



Zdzislaw Meglicki
2001-02-26