L-20 MCS 471 Friday 22 February 2002: Lagrange and Neville Interpolation

This worksheet illustrates some of the concepts we discussed in class.

1. The Interpolation Problem

We show that the problem of finding an interpolating polynomial through (x[i],f[i]), i=0,1,..,n is equivalent to solving a linear system.

> n := 3;

The interpolating polynomial has unknown coefficients a[i]:

> p := sum(a['i']*z^'i','i'=0..n);

The interpolation conditions are

> for i from 0 to n do

> eqs[i] := subs(z=x[i],p) = f[i];

> od;

Now we extract the matrix of this linear system

> unknowns := [seq(a[i],i=0..n)];

> listeqs := convert(eqs,list):

> A:= linalg[genmatrix](listeqs,unknowns);

and compute the determinant of the matrix (in factored form):

> factor(linalg[det](A));

We see that the determinant is different from zero provided x[i] /= x[j], for i /= j.

2. Lagrange Interpolation

One solution to the interpolation problem is to take a special form of the interpolating polynomial:

> q := sum(l['i']*f['i'],'i'=0..n);

where l[i] is defined as follows:

> for i from 0 to n do

> l[i] := product((z-x['j'])/(x['i']-x['j']),'j'=0..i-1)*product((z-x['j'])/(x['i']-x['j']),'j'=i+1..n);

> od;

This is how our interpolating polynomial looks like:

> q;

We turn the expression for q into a function for evaluation:

> funq := unapply(q,z):

Checking the interpolation conditions:

> for i from 0 to n do

> funq(x[i]);

> od;

3. Neville Interpolation

One shortcoming of Lagrange interpolation is that everything has to be recalculated when an additional data point is used. Here we present the idea behind Neville interpolation.

Above we constructed q as the interpolating polynomial through (x[i],f[i]), i=0,1,..,n. Let us construct the interpolating polynomial through (x[i],f[i]), i = 1,2,..,n+1.

> q1 := sum(k['i']*f['i'],'i'=1..n+1);

where k is defined in similar fashion as the l[i] above:

> for i from 1 to n+1 do
k[i] := product((z-x['j'])/(x['i']-x['j']),'j'=1..i-1)*product((z-x['j'])/(x['i']-x['j']),'j'=i+1..n+1);
od;

Then the interpolating polynomial looks like

> q1;

As q interpolates through (x[i],f[i]), i = 0,1,..,n, and q1 through (x[i],f[i]), i=1,2,..,n+1, we now claim that the following polynomial interpolates through all n+2 points:

> qq1 := (z-x[n+1])/(x[0]-x[n+1])*'q' + (z-x[0])/(x[n+1]-x[0])*'q1';

Or fully evaluated:

> eqq1 := map(t->eval(t),qq1);

To check the claim we create a function and evaluate at the interpolation points:

> funqq1 := unapply(eqq1,z):

> for i from 0 to n+1 do

> simplify(funqq1(x[i]));

> od;

This leads to a general interpolation method. To be continued...