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