PHYS 251 - Introduction to Computer Techniques in Physics

Polynomial Interpolation

Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence preceding the mathematical presentation of the equation. They should not be read as if they are part of the sentence. When reference is made to an equation, the word "equation" will be spelled out and not abbreviated. Mathematical symbols are not in all cases precise. Greek symbols will show as such on a Windows operating system, but not a Unix or Mac system.

Polynomial Interpolation

The justification for approximating with polynomials rest on a theorem by Weierstrass in 1885, to which

"...every continuous function in an interval (a,b) can be represented in that interval to any desired accuracy by a polynomial."

One can also show that the following propositions are true for a polynomial

1. nth differences of an nth-degree polynomial are constant for an arithmetic progression (equal interval) of the independent variable,
2. if nth differences of a tabulated function for equal interval data are constant, the unknown function is a polynomial of degree n.

Example

Let f(x) = x² + 3x + 2. Forming differences in a horizontal difference table, we have the table shown below.

xi	fi(x)	D1fi(x)	D2fi(x)	D3fi(x)	D4fi(x)
0.000000	2.000000
1.000000	6.000000	4.000000
2.000000	12.000000	6.000000	2.000000
3.000000	20.000000	8.000000	2.000000	0.000000
4.000000	30.000000	10.000000	2.000000	0.000000
5.000000	42.000000	12.000000	2.000000	0.000000
6.000000	56.000000	14.000000	2.000000	0.000000

Since the function is indeed a 2nd degree polynomial the 2nd differences are constant.

Polynomial Interpolation Retaining Second Differences

If we restrict our consideration to the retention of second differences and ignore all higher-order differences, there are four useful polynomial approximations that can be derived for equal-interval interpolation.

Newton's Forward Interpolation: Given values of a function f(x) at a finite number of discrete points (n),

x0,	x1,	x2,	x3,	...,	xn-1,	xn
f(x0) = f0,	f(x1) = f1,	f(x2) = f2,	f(x3) = f3,	...,	f(xn-1) = fn-1,	f(xn) = fn

and letting u = [(x - x_i) / (x_i+1 - x_i)], then Newton's forward interpolation formula can be expressed using difference notation as

$p(x)\; =\; p(x$₀ + uDx)

or substituting the difference values of the function, we have

$p(x)\; =\; f$₀ + uD₁f₁ + [u(u-1)/2!]D₂f₂ + [u(u-1)(u-2)/3!]D₃f₃ +...+ [u(u-1)(u-2)...(u-n+1)/n!]D_nf_n +...

which can be simplified by retaining only second order differences and substituting the functional values for the differences (Eq. 1)

$p(x)\; =\; f$₀ + u ( f₁ - f₀ ) + [ u ( u - 1 ) / 2 ] ( f₂ - 2 f₁ + f₀ )

Since Newton's forward interpolation formula uses values of the function at points ahead in the table, it can be used at the beginning of the tabulated values, but not at the end.

Newton's Backward Interpolation: Given values of a function f(x) at a finite number of discrete points (n), and letting u = [(x - x_i) / (x_i - x_i-1)], then Newton's backward interpolation formula can be expressed using difference notation as

$p(x)\; =\; p(x$₀ + uDx)

or substituting the difference values of the function, we have

$p(x)\; =\; f$_n + uD₁f_n + [u(u+1)/2!]D₂f_n + [u(u+1)(u+2)/3!]D₃f_n +...+ [u(u+1)(u+2)...(u+n-1)/n!]D_nf_n +...

which can be simplified by retaining only second order differences and substituting the functional values for the differences (Eq. 2)

$p(x)\; =\; f$₀ + u ( f₀ - f_-1 ) + [ u ( u + 1 ) / 2 ] ( f₀ - 2 f_-1 + f_-2 )

which is useful for interpolating near the end of the tabulated functions.

Stirling's Central-Difference Interpolation: This interpolation formula is based on a diagonal difference table rather than a horizontal difference table. Without going in to the distinctions between the two types of difference tables, we can write down an interpolation formula including only second differences that is similar to Newton's interpolation formulas. Thus given values of a function f(x) at a finite number of discrete points (n), and letting u = [(x - x_i) / (x_i+1 - x_i)], Stirling's central-difference interpolation formula can be expressed using functional values as (Eq. 3)

$p(x)\; =\; f$₀ + ( u / 2 ) ( f₁ - f_-1 ) + [ u² / 2 ] ( f₁ - 2 f₀ + f_-1 )

which also uses three successive points to approximate the function as does Newton.

Bessel's Central-Difference Interpolation: This interpolation formula is also based on a diagonal difference table rather than a horizontal difference table. We shall also write down an interpolation formula including only second differences that is similar to Newton's and Stirling's interpolation formulas. Thus given values of a function f(x) at a finite number of discrete points (n), and letting u = [(x - x_i) / (x_i+1 - x_i)], Bessel's central-difference interpolation formula can be expressed using functional values as (Eq. 4)

$p(x)\; =\; f$₀ + u ( f₁ - f₀ ) + [ u ( u - 1 ) / 4 ] ( f₂ - f₁ - f₀ + f_-1 )

which uses four successive points rather than three to approximate the function. Note that in equations 1, 2, 3, and 4, there is a zero-order term, a first-order term, and a second-order term. The zero- and first-order terms constitute essential a linear interpolation, while the second-order term is a correction term to that linear interpolation.

Example

Find f(x) = exp(-x), at the point x = 1.7565, given the following values of the function.

i	xi	fi
-2	1.730000000	0.177284100
-1	1.740000000	0.175520401
0	1.750000000	0.173773944
1	1.760000000	0.172044864
2	1.770000000	0.170332989

Forming the value of u for the four equidistance interpolation formulas, we have for Newton's forward, Stirling's, and Bessel's interpolation formulas

$u\; =\; (\; x\; -\; x$₀ ) / ( x₁ - x₀ ) = ( 1.7565 - 1.75 ) / ( 1.76 - 1.75 ) = 0.0065 / 0.01 = + 0.65 ,

Newton's Forward Interpolation:

$u\; (\; u\; -\; 1\; )\; /\; 2\; =\; 0.65\; (\; 0.65\; -\; 1\; )\; /\; 2\; =\; -\; 0.11375\; ,$

Stirling's Central-Difference Interpolation:

$u\; /\; 2\; =\; 0.65\; /\; 2\; =\; 0.325\; ,$

$u2/\; 2\; =\; (\; 0.65\; )2/\; 2\; =\; 0.21125\; ,$

Bessel's Central-Difference Interpolation:

$u\; (\; u\; -\; 1\; )\; /\; 4\; =\; 0.65\; (\; 0.65\; -\; 1\; )\; /\; 4\; =\; -\; 0.056875\; .$

For Newton's backward formula, the value of u is given by

$u\; =\; (\; x\; -\; x$₁ ) / ( x₁ - x₀ ) = ( 1.7565 - 1.76 ) / ( 1.76 - 1.75 ) = - 0.0035 / 0.01 = - 0.35 ,

Newton's Backward Interpolation:

$u\; (\; u\; +\; 1\; )\; /\; 2\; =\; -\; 0.35\; (\; -\; 0.35\; +\; 1\; )\; /\; 2\; =\; -\; 0.11375\; .$

The following table evaluates the three terms that appear in each of the four interpolation formula. Note that in order to apply Newton's backward formula, we must re-label the points so that i = -1, 0, 1 become i = -2, -1, 0.

Interpolation Formula	Zero-Order Term	First-Order Term	Second-Order Term	f(u)
Newton's Forward	0.173773944	-0.001123902	-0.000001957	0.172648085
Newton's Backward	0.172044864	+0.000605178	-0.000001977	0.172648065
Stirling's Formula	0.173773944	-0.001129550	+0.000003671	0.172648065
Bessel's Formula	0.173773944	-0.001123902	-0.000001967	0.172648075
Exact Solution				0.172648076

It is evident that Bessel's interpolation formula has given us the most accurate result. This is inpart true because Bessel's formula is known to be more accurate when u lies between 0.25 and 0.75.

Non-Equidistant Interpolation

A frequently encounter interpolation problem is one in which the data is not equally spaced in the independent variable. This leads us to consider Lagrangian interpolation in the following document on Non-Equidistant Interpolation.