Numerical methods. Note 13
[~Home~] Interpolation and Extrapolation. Numerical Methods. Note~ « 13 »

Polynomial interpolation. Lagrange interpolating polynomial Pn-1(x) is the polynomial of the order n-1 that passes through the n given points (xi,yi), i=1..n

P(x) = ∑i yi ∏k≠i (x-xk)/(xi-xk)~.
Might show erratic behaviour for higher degree polynomials. Not useful for n larger than, say, 5.

Rational function interpolation. Here the rational function is a ration of two polynomials Pm(x) and ,Qk(x),

R(x) = Pm(x)/Qk(x) .
It has potentially larger range of conversion than polynomial. Still dangerous for higher orders.

Spline interpolation. At the interval between the tabulated points i and i+1 the function is approximated by a polynomial:

y[i](x) = yi + bi (x-xi) + ci (x-xi)2 + di (x-xi)3 + ...~.
The coefficients bi, ci, di,... are calculated from the continuity equations for the function and its derivatives at the boundary
y[i-1](xi)=y[i](xi),~ y[i-1]'(xi)=y[i]'(xi),~ y[i-1]''(xi)=y[i]''(xi),~...
For the quadratic spline the equations become
yi+bih+cih2 = yi+1 (continuity of the function at the i+1st point),
bi + 2cih = bi+1 (continuity of the derivative)
The most popular spline is the cubic spline. Here are the fortran routines from netlib.org: spline.f - for constructing the spline (calculating coefficients b,c,d) and seval.f - for the very interpolation.

Integration and differentiation of tabulated functions. One of the ways is first to interpolate the tabulated data and then integrate or differentiate the intepolating function.

Problems

  1. Make a function locate(x,xt[]) which, given an ordered array xt[] and a number x, finds index j such that xt[j]~<~x~<~xt[j+1].
  2. Make a function pinterp(x,xt[],yt[],k) which does a polynomial interpolation of a table (xt[],yt[]) at a point x using k nearest tabulated points.
  3. Make a function for the quadratic spline interpolation.
  4. Make a Lagrange polynomial interpolation and a spline interpolation of the following tabulated function (notice that i changes from ZERO (!) to 10):
    xi = i + sin(i)/2, yi = i + cos(i2), i=0,10
    Make a plot of the result.
"Copyleft" © 2004 D.V.Fedorov (fedorov at phys dot au dot dk)