Computational Physics - Department of Physics

3.2 Numerical Interpolation and Extrapolation 61

3.2.1 Interpolation

Let us assume that we have a set ofN+ 1 pointsy 0 =f(x 0 ),y 1 =f(x 1 ),...,yN=f(xN)where
none of thexivalues are equal. We wish to determine a polynomial of degreenso that

PN(xi) =f(xi) =yi, i= 0 , 1 ,...,N (3.7)

for our data points. If we then writePNon the form

PN(x) =a 0 +a 1 (x−x 0 )+a 2 (x−x 0 )(x−x 1 )+···+aN(x−x 0 )...(x−xN− 1 ), (3.8)

then Eq. (3.7) results in a triangular system of equations

a 0 =f(x 0 )

a 0 + a 1 (x 1 −x 0 ) =f(x 1 )

a 0 +a 1 (x 2 −x 0 )+a 2 (x 2 −x 0 )(x 2 −x 1 ) =f(x 2 )

... ... ... ...

.

The coefficientsa 0 ,...,aNare then determined in a recursive way, starting witha 0 ,a 1 ,....
The classic of interpolation formulae was created by Lagrange and is given by

PN(x) =

N

∑

i= 0

∏

k 6 =i

x−xk

xi−xk

yi. (3.9)

If we have just two points (a straight line) we get

P 1 (x) =

x−x 0

x 1 −x 0

y 1 +

x−x 1

x 0 −x 1

y 0 ,

and with three points (a parabolic approximation) we have

P 3 (x) =

(x−x 0 )(x−x 1 )

(x 2 −x 0 )(x 2 −x 1 )

y 2 +

(x−x 0 )(x−x 2 )

(x 1 −x 0 )(x 1 −x 2 )

y 1 +

(x−x 1 )(x−x 2 )

(x 0 −x 1 )(x 0 −x 2 )

y 0

and so forth. It is easy to see from the above equations that whenx=xiwe have thatf(x) =
f(xi)It is also possible to show that the approximation error (or rest term) is given by the
second term on the right hand side of

f(x) =PN(x)+

ωN+ 1 (x)f(N+^1 )(ξ)

(N+ 1 )!

. (3.10)

The functionωN+ 1 (x)is given by

ωN+ 1 (x) =aN(x−x 0 )...(x−xN),

andξ=ξ(x)is a point in the smallest interval containing all interpolation pointsxjandx.
The program we provide below is however based on divided differences. The recipe is quite
simple. If we takex=x 0 in Eq. (3.8), we then have obviously thata 0 =f(x 0 ) =y 0. Movinga 0
over to the left-hand side and dividing byx−x 0 we have

f(x)−f(x 0 )

x−x 0

=a 1 +a 2 (x−x 1 )+···+aN(x−x 1 )(x−x 2 )...(x−xN− 1 ),

where we hereafter omit the rest term