Lagrange Interpolations (Lec 10)

yukita@k.hosei.ac.jp

[Graphics:HTMLFiles/index_1.gif]

Interpolation (or approximation)  is understood as  a well chosen embedding of a finite dimensional vector space,  a control parameter space, into  a  function space. You can't grasp the mathematical substance of the theory without having this view point.

•Lagrange Polynomial Interpolation

Problem:  Given  parameter knots {t _ 0, t _ 1, ...., t _ n} and data points P _ 0(x _ 0, y _ 0), ..., P _ n(x _ n, y _ n),  find a curve  that satisfies the condition  P(t _ i)=P _ i.

We include information about  the set of knots  knots={knots[0], ...,knots[n]}.  Integer  n  indicates the last index or the length of the knots plus one. Integer n is also the degree of Lagrange polonomials that are to be generated.

Definition:

L[knots_,n_][i_][t_] :=
Product[ (t - knots[j])/(knots[i] - knots[j]),
         {j,0,i-1} ] *
Product[ (t - knots[j])/(knots[i] - knots[j]),
         {j,i+1,n} ]

Let's see how the definition works in the case degree n=3.

Clear[knots]

Array[knots, 4, 0]

{knots[0], knots[1], knots[2], knots[3]}

We would like a parallel assignment like
{knots[0],knots[1],knots[2],knots[3]}={0,1,2,3}.
To do this in batch, you might consider the following input.
Array[knots,4,0]=Range[0,3]
However, this will produce an error. To avoid it, we must evaluate the LHS prior to the assignment.

Evaluate[Array[knots, 4, 0]] = Range[0, 3]

{0, 1, 2, 3}

Let us check the parallel assignment was safely done.

Map[knots, Range[0, 3]]

{0, 1, 2, 3}

L[knots,3][0][t]

1/6 (1 - t) (2 - t) (3 - t)

Map[ L[knots,3][0], Range[0,3]]

{1, 0, 0, 0}

L[{1,2,3,4}][1][t]

L[{1, 2, 3, 4}][1][t]

Map[ L[knots,3][1], Range[0,3]]

{0, 1, 0, 0}

Let's do it in batch.

Table[ L[knots,3][i][j],
       {i,0,3},{j,0,3}] // MatrixForm

( 1   0   0   0 )    0   1   0   0    0   0   1   0    0   0   0   1

Another way to produce the same matrix is

Outer[ L[knots, 3][#1][#2] &, Range[0, 3], Range[0, 3]] // MatrixForm

( 1   0   0   0 )    0   1   0   0    0   0   1   0    0   0   0   1

Show[ GraphicsArray[
        Table[Plot[{L[knots,3][2i+j][t],1},
                   {t,-0.5,3.5},
                   PlotRange->{-2,2},
                   Ticks->{Automatic,{-1,0,1}},
                   PlotLabel->"L[..]["<>
                               ToString[2i+j]<>
                              "]",
                   DisplayFunction->Identity],
             {i,0,1},{j,0,1}] ],
      DisplayFunction->$DisplayFunction];

[Graphics:HTMLFiles/index_27.gif]

We consider the two linear mappings.

interp : R^(n + 1)-->  R[x]
      
(y _ 0, ···, y _ n) |->  f(x)=Underoverscript[∑, i = 0, arg3] y _ i  L _ (knots, i) (x)
(* pseudo code, don't evaluate it *)

where R[x] is a ring of polynomials over the real number field R.

proj : R[x] --> R^(n + 1)
         f(x) |-> ( f(x _ 0), ··· , f(x _ n) )
(* pseudo code, don't evaluate this *)

What we have shown is

proj o interp : R^(n + 1)-->R^(n + 1)
(* pseudo code, don't evaluate this *)

is the identity map.

Let's look at a concrete example.

•Example 1 -- interpolating  eight  points

data={{0,0},{1,3},{3,3},{5,-2},{6,-2},{7,+1},{9,4},{10,4}};

Clear[cps, knots]

Let us prepare the control points.

Array[cps, 8, 0]

{cps[0], cps[1], cps[2], cps[3], cps[4], cps[5], cps[6], cps[7]}

Evaluate[Array[cps, 8, 0]] = data ;

We define the set of knots.

Array[knots, 8, 0]

{knots[0], knots[1], knots[2], knots[3], knots[4], knots[5], knots[6], knots[7]}

Evaluate[Array[knots, 8, 0]] = Range[0, 7] ;

P[t_]=Sum[ cps[i] L[knots,7][i][t], {i,0,7} ]//Expand

{(349 t)/210 - (433 t^2)/180 + (79 t^3)/30 - (155 t^4)/144 + (49 t^5)/240 - (13 t^6)/720 + t^7 ... (568 t)/35 + (3173 t^2)/72 - (1717 t^3)/48 + (475 t^4)/36 - (299 t^5)/120 + (17 t^6)/72 - t^7/112}

g1=ParametricPlot[Evaluate[P[t]], {t,0,7},
                  DisplayFunction->Identity,
                  PlotLabel->"Lagrange"];

g0=Graphics[ Table[ Circle[P[j],0.1 ], {j,0,7} ]];

Show[{g1,g0}, DisplayFunction->$DisplayFunction ];

[Graphics:HTMLFiles/index_47.gif]

The above interpolation is not so beautiful, is it?

•Example  2 -- interpolating a sine curve

theXs={-Pi/4, 0, Pi/4, Pi/2, 3Pi/4, Pi, 5Pi/4 };

data=Thread[{ theXs , Sin[theXs] }]

{{-π/4, -1/2^(1/2)}, {0, 0}, {π/4, 1/2^(1/2)}, {π/2, 1}, {(3 π)/4, 1/2^(1/2)}, {π, 0}, {(5 π)/4, -1/2^(1/2)}}

Clear[cps, knots]

Let us prepare the control points.

Array[cps, Length[theXs], 0]

{cps[0], cps[1], cps[2], cps[3], cps[4], cps[5], cps[6]}

Evaluate[Array[cps, Length[theXs], 0]] = data ;

We define the set of knots.

Array[knots, Length[theXs], 0]

{knots[0], knots[1], knots[2], knots[3], knots[4], knots[5], knots[6]}

Evaluate[Array[knots, Length[theXs], 0]] = theXs ;

The interpolation is given as

P[t_]=Sum[ cps[i] L[knots,Length[theXs]-1][i][t], {i,0,Length[theXs]-1} ];

ParametricPlot[Evaluate[P[t]], {t,-Pi/2, 3Pi/2},
               Ticks->{ theXs, Automatic } ];

[Graphics:HTMLFiles/index_56.gif]

It looks like a sine curve. For comparison we give here the linear interpolation plot.

Show[ Graphics[ Line[ data ],
                Axes->True,
                Ticks->{ theXs, {-1,1}},
                PlotRange->{{-1.8,4.8},{-1,1}}] ];

[Graphics:HTMLFiles/index_57.gif]

•Exercises

•Problem 1

Plot the Lagrange interpolation with the control points in  Example 1. Use the  functional tricks suggested in the hints below.

•Hints

Map[t - # &, Range[0, 9]]

{t, -1 + t, -2 + t, -3 + t, -4 + t, -5 + t, -6 + t, -7 + t, -8 + t, -9 + t}

Drop[Map[t - # &, Range[0, 9]], {3 + 1}]

{t, -1 + t, -2 + t, -4 + t, -5 + t, -6 + t, -7 + t, -8 + t, -9 + t}

Fold[Times, 1, Drop[Map[t - # &, Range[0, 9]], {3 + 1}]]

(-9 + t) (-8 + t) (-7 + t) (-6 + t) (-5 + t) (-4 + t) (-2 + t) (-1 + t) t

You should not do the following. This is just for inspiration.

Clear[L, knots, cps]

Evaluate[Array[knots, 10, 0]] = Range[0, 9]

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

Clear[L]

Array[knots, 10, 0]

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

L[knots, 9][3][t_] =  Fold[Times, 1, Drop[Map[t - # &, Array[knots, 10, 0]], {3 + 1}]]/ Fold[Times, 1, Drop[Map[knots[3] - # &, Array[knots, 10, 0]], {3 + 1}]]

((-9 + t) (-8 + t) (-7 + t) (-6 + t) (-5 + t) (-4 + t) (-2 + t) (-1 + t) t)/4320

Map[L[knots, 9][3], Range[0, 9]]

{0, 0, 0, 1, 0, 0, 0, 0, 0, 0}

•Solution

Converted by Mathematica  (June 3, 2003)