Computational Physics - Department of Physics

62 3 Numerical differentiation and interpolation

f(N+^1 )(ξ)

(N+ 1 )!

(x−x 1 )(x−x 2 )...(x−xN).

The quantity

f 0 x=

f(x)−f(x 0 )

x−x 0

,

is a divided difference of first order. If we then takex=x 1 , we have thata 1 =f 01. Movinga 1
to the left again and dividing byx−x 1 we obtain

f 0 x−f 01

x−x 1

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

and the quantity

f 01 x=

f 0 x−f 01

x−x 1

,

is a divided difference of second order. We note that the coefficient

a 1 =f 01 ,

is determined fromf 0 xby settingx=x 1. We can continue along this line and define the divided
difference of orderk+ 1 as

f 01 ...kx=

f 01 ...(k− 1 )x−f 01 ...(k− 1 )k

x−xk , (3.11)

meaning that the corresponding coefficientakis given by

ak=f 01 ...(k− 1 )k.

With these definitions we see that Eq. (3.10) can be rewrittenas

f(x) =a 0 +∑

k= 1

N f 01 ...k(x−x 0 )...(x−xk− 1 )+

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

(N+ 1 )!.

If we replacex 0 ,x 1 ,...,xkin Eq. (3.11) withxi+ 1 ,xi+ 2 ,...,xk, that is we count fromi+ 1 tok
instead of counting from 0 tokand replacexwithxi, we can then construct the following
recursive algorithm for the calculation of divided differences

fxixi+ 1 ...xk=

fxi+ 1 ...xk−fxixi+ 1 ...xk− 1

xk−xi

.

Assuming that we have a table with function values(xj,f(xj) =yj)and need to construct the
coefficients for the polynomialPN(x). We can then view the last equation by constructing the
following table for the case whereN= 3.

x 0 y 0

fx 0 x 1

x 1 y 1 fx 0 x 1 x 2

fx 1 x 2 fx 0 x 1 x 2 x 3

x 2 y 2 fx 1 x 2 x 3

fx 2 x 3

x 3 y 3

.

The coefficients we are searching for will then be the elements along the main diagonal.
We can understand this algorithm by considering the following. First we construct the unique
polynomial of order zero which passes through the pointx 0 ,y 0. This is justa 0 discussed above.