Polynomial Interpolation

Sergey Fomel   02/19/2002

Example Function

An example function for studing polynomial interpolation is

f[x_] := 1/x

fplot  =  
        Plot[f[x], {x, 1, 6}, PlotStyle -> {Thickness[0.01]}, 
          Frame -> True, PlotRange -> {0, 1.1}];

We select three nodes and find the interpolating polynomial for them.

nodes  =  {2, 4, 5};

tbl  =  Transpose[{nodes, Map[f, nodes]}];

tblplot  =  
        ListPlot[tbl, PlotRange -> {{1, 6}, {0, 1.1}}, 
          Frame -> True, PlotStyle -> PointSize[0.02]];

Mathematica can interpolate points using Interpolation function:

dsk[x_] := Disk[x, {0.05, 0.011}]

<< Graphics`Legend`

intplot = 
        Plot[{f[x], Interpolation[tbl, InterpolationOrder -> 2][
              x]}, {x, 1, 6}, Frame -> True, PlotRange -> {0, 1.1}, 
          PlotStyle -> {Dashing[{0.01}], Dashing[{0}]}, 
          PlotLegend -> {"True", "Interpolation"}, 
          LegendPosition -> {0, 0}, 
          Epilog -> Map[dsk, tbl]];

The second-order interpolation polynomial does not reconstruct the
function exactly but fits the input data.

Lagrange Interpolation

The Lagrange form of the interpolation polynomial is

Lk[x_, xk_, k_] := 
      Product[(x - xk[[i]])/(xk[[k]] - xk[[i]]), {i, 1, 
            k - 1}]*
  Product[(x - 
                xk[[i]])/(xk[[k]] - xk[[i]]), {i, k + 1, 
            Length[xk]}]

Lag[x_] := 
      Sum[f[nodes[[k]]]  Lk[x, nodes, k], {k, 
          Length[nodes]}]

Lag[x]

Expand[%]

plotlk[k_] := 
      Plot[Lk[x, nodes, k], {x, 1, 6}, PlotStyle -> {Thickness[0.01]}, 
        DisplayFunction -> Identity]

<< Graphics`Arrow`

Show[Table[plotlk[k], {k, Length[nodes]}], 
        Frame -> 
          True, 
PlotLabel -> "Lagrange \
Polynomials", DisplayFunction -> $DisplayFunction, 
        Epilog -> {Arrow[{2, 2}, {1.5, Lk[1.5, nodes, 1] + 0.1}], 
            Text["L1", {2, 2}, {-1.5, 0}], 
            Arrow[{3.5, 2.2}, {3.5, Lk[3.5, nodes, 2] + 0.1}], 
            Text["L2", {3.5, 2.2}, {0, -1}], 
            Arrow[{4.5, 2}, {5, Lk[5, nodes, 3] + 0.1}], 
            Text["L3", {4.5, 2}, {1, -1}]}];

L1(x) is one at x1 and zero at x2 and x3. Likewise, L2(x) is one at x2 and zero at x1 and x3. L3(x) is one at x3 and zero at x1 and x2.

Putting it all together:

Plot[Lag[x], {x, 1, 6}, PlotRange -> {0, 1.1}, 
        Frame -> True, Epilog -> Map[dsk, tbl], 
        PlotLabel -> "Lagrange Interpolation"];

Newton Interpolation

The Newton form of the interpolation polynomial is

Nk[x_, xk_, k_] := Product[x - xk[[i]], {i, k - 1}]

divdif[k_, 0] :=  f[nodes[[k]]]

divdif[k_, 
        n_] := (divdif[k + 1, n - 1] - 
            divdif[k, n - 1])/(nodes[[k + n]] - 
            nodes[[k]])

New[x_] := 
      Sum[divdif[1, k - 1] Nk[x, nodes, k], {k, 1, 
          Length[nodes]}]

New[x]

Simplify[New[x] - Lag[x]]

plotnk[k_] := 
      Plot[Nk[x, nodes, k], {x, 1, 6}, PlotStyle -> {Thickness[0.01]}, 
        DisplayFunction -> Identity]

Show[Table[plotnk[k], {k, Length[nodes]}], Frame -> True, 
        PlotLabel -> "Newton Polynomials", 
        DisplayFunction -> $DisplayFunction, 
        Epilog -> {Arrow[{2, 3.5}, {2, Nk[2, nodes, 1] + 0.2}], 
            Text["N1", {2, 3.5}, {0, -1}], 
            Arrow[{3.5, 4}, {4, Nk[4, nodes, 2] + 0.2}], 
            Text["N2", {3.5, 4}, {0, -1}], 
            Arrow[{4.5, 6}, {5.45, Nk[5.5, nodes, 3]}], 
            Text["N3", {4.5, 6}, {1, 0}]}];

Putting it all together:

Plot[New[x], {x, 1, 6}, PlotRange -> {0, 1.1}, 
        Frame -> True, Epilog -> Map[dsk, tbl], 
        PlotLabel -> "Newton Interpolation"];

Neville Interpolation

The Neville form of the interpolation polynomial is defined by recursion

Nev[x_, k_, 0] := f[nodes[[k]]]

Nev[x_, k_, 
        n_] := (Nev[x, k + 1, n - 1] (x - nodes[[k]]) + 
            Nev[x, k, 
                n - 1] (nodes[[k + n]] - x))/(nodes[[k + 
                  n]] - nodes[[k]])

Nev[x, 1, 2]

Simplify[Nev[x, 1, 2] - Lag[x]]

plotnv[k_, n_] := 
      Plot[Nev[x, k, n], {x, 1, 6}, 
        DisplayFunction -> Identity]

Show[Table[plotnv[k, 0], {k, Length[nodes]}], Frame -> True, 
        PlotLabel -> "Neville Polynomials, Order 0", 
        DisplayFunction -> $DisplayFunction, PlotRange -> {0, 1.1}, 
        Epilog -> Map[dsk, tbl]];

Show[Table[plotnv[k, 1], {k, Length[nodes] - 1}], Frame -> True, 
        PlotLabel -> "Neville Polynomials, Order 1", 
        DisplayFunction -> $DisplayFunction, PlotRange -> {0, 1.1}, 
        Epilog -> Map[dsk, tbl]];

Show[plotnv[1, 2], Frame -> True, 
        PlotLabel -> "Neville Interpolation", 
        DisplayFunction -> $DisplayFunction, PlotRange -> {0, 1.1}, 
        Epilog -> Map[dsk, tbl]];