3.2 Linear Algebra3.4 Differentiation

§3.3 Interpolation

Contents

§3.3(i) Lagrange Interpolation

The nodes or abscissas z_{k} are real or complex; function values are f_{k}=f(z_{k}). Given n+1 distinct points z_{k} and n+1 corresponding function values f_{k}, the Lagrange interpolation polynomial is the unique polynomial P_{n}(z) of degree not exceeding n such that P_{n}(z_{k})=f_{k}, k=0,1,\dots,n. It is given by

3.3.1P_{n}(z)=\sum _{{k=0}}^{n}\ell _{k}(z)f_{k}=\sum _{{k=0}}^{n}\frac{\omega _{{n+1}}(z)}{(z-z_{k})\omega _{{n+1}}^{{\prime}}(z_{k})}f_{k},

where

3.3.2
\ell _{k}(z)={\sideset{}{{}^{{\prime}}}{\prod}_{{j=0}}^{n}}\frac{z-z_{j}}{z_{k}-z_{j}},
\ell _{k}(z_{j})=\delta _{{k,j}}.

Here the prime signifies that the factor for j=k is to be omitted, \delta _{{k,j}} is the Kronecker symbol, and \omega _{{n+1}} is the nodal polynomial

3.3.3\omega _{{n+1}}(z)=\prod _{{k=0}}^{n}(z-z_{k}).

With an error term the Lagrange interpolation formula for f is given by

3.3.4f(z)=\sum _{{k=0}}^{n}\ell _{{k}}(z)f_{k}+R_{n}(z).

If f, x (=z), and the nodes x_{k} are real, and f^{{(n+1)}} is continuous on the smallest closed interval I containing x,x_{0},x_{1},\dots,x_{n}, then the error can be expressed

3.3.5R_{n}(x)=\frac{f^{{(n+1)}}(\xi)}{(n+1)!}\,\omega _{{n+1}}(x),

for some \xi\in I. If f is analytic in a simply-connected domain D1.13(i)), then for z\in D,

3.3.6R_{n}(z)=\frac{\omega _{{n+1}}(z)}{2\pi i}\int _{C}\frac{f(\zeta)}{(\zeta-z)\omega _{{n+1}}(\zeta)}d\zeta,

where C is a simple closed contour in D described in the positive rotational sense and enclosing the points z,z_{1},z_{2},\dots,z_{n}.

§3.3(ii) Lagrange Interpolation with Equally-Spaced Nodes

The (n+1)-point formula (3.3.4) can be written in the form

3.3.7f_{t}=f(x_{0}+th)=\sum _{{k=n_{0}}}^{{n_{1}}}A_{k}^{n}f_{k}+R_{{n,t}},n_{0}<t<n_{1},

where the nodes x_{k}=x_{0}+kh (h>0) and function f are real,

3.3.8
n_{0}=-\tfrac{1}{2}(n-\sigma),
n_{1}=\tfrac{1}{2}(n+\sigma),
3.3.9\sigma=\tfrac{1}{2}(1-(-1)^{n}),

and A_{k}^{n} are the Lagrangian interpolation coefficients defined by

3.3.10A_{k}^{n}=\frac{(-1)^{{n_{1}+k}}}{(k-n_{0})!\,(n_{1}-k)!(t-k)}\,\prod _{{m=n_{0}}}^{{n_{1}}}(t-m).

The remainder is given by

3.3.11R_{{n,t}}=R_{n}(x_{0}+th)=\frac{h^{{n+1}}}{(n+1)!}f^{{(n+1)}}(\xi)\prod _{{k=n_{0}}}^{{n_{1}}}(t-k),

where \xi is as in §3.3(i).

Let c_{n} be defined by

3.3.12c_{n}=\frac{1}{(n+1)!}\max\prod _{{k=n_{0}}}^{{n_{1}}}\left|t-k\right|,

where the maximum is taken over t-intervals given in the formulas below. Then for these t-intervals,

3.3.13\left|R_{{n,t}}\right|\leq c_{n}h^{{n+1}}\left|f^{{(n+1)}}(\xi)\right|.

Linear Interpolation

3.3.14f_{t}=(1-t)f_{0}+tf_{1}+R_{{1,t}},0<t<1,
3.3.15c_{1}=\tfrac{1}{8},0<t<1.

Three-Point Formula

3.3.16f_{t}=\sum _{{k=-1}}^{1}A_{k}^{2}f_{k}+R_{{2,t}},\left|t\right|<1,
3.3.17
A_{{-1}}^{2}=\tfrac{1}{2}t(t-1),
A_{0}^{2}=1-t^{2},
A_{1}^{2}=\tfrac{1}{2}t(t+1),
3.3.18c_{2}=1/(9\sqrt{3})=0.0641\ldots,\left|t\right|<1.

Four-Point Formula

3.3.19f_{t}=\sum _{{k=-1}}^{2}A_{k}^{3}f_{k}+R_{{3,t}},-1<t<2,
3.3.20
A_{{-1}}^{3}=-\tfrac{1}{6}t(t-1)(t-2),
A_{0}^{3}=\tfrac{1}{2}(t^{2}-1)(t-2),
A_{1}^{3}=-\tfrac{1}{2}t(t+1)(t-2),
A_{2}^{3}=\tfrac{1}{6}t(t^{2}-1),
3.3.21c_{3}=\begin{cases}\tfrac{3}{128}=0.0234\ldots,&0<t<1,\\
\tfrac{1}{24}=0.0416\ldots,&-1<t<0\mbox{ or }1<t<2.\\
\end{cases}

Five-Point Formula

3.3.22f_{t}=\sum _{{k=-2}}^{2}A_{k}^{4}f_{k}+R_{{4,t}},\left|t\right|<2,
3.3.23
A_{{-2}}^{4}=\tfrac{1}{24}t(t^{2}-1)(t-2),
A_{{-1}}^{4}=-\tfrac{1}{6}t(t-1)(t^{2}-4),
A_{0}^{4}=\tfrac{1}{4}(t^{2}-1)(t^{2}-4),
A_{1}^{4}=-\tfrac{1}{6}t(t+1)(t^{2}-4),
A_{2}^{4}=\tfrac{1}{24}t(t^{2}-1)(t+2),
3.3.24c_{4}=\begin{cases}0.0118\ldots,&\left|t\right|<1,\\
0.0302\ldots,&1<\left|t\right|<2.\\
\end{cases}

Six-Point Formula

3.3.25f_{t}=\sum _{{k=-2}}^{3}A_{k}^{5}f_{k}+R_{{5,t}},-2<t<3,
3.3.26
A_{{-2}}^{5}=-\tfrac{1}{120}t(t^{2}-1)(t-2)(t-3),
A_{{-1}}^{5}=\tfrac{1}{24}t(t-1)(t^{2}-4)(t-3),
A_{0}^{5}=-\tfrac{1}{12}(t^{2}-1)(t^{2}-4)(t-3),
A_{1}^{5}=\tfrac{1}{12}t(t+1)(t^{2}-4)(t-3),
A_{2}^{5}=-\tfrac{1}{24}t(t^{2}-1)(t+2)(t-3),
A_{3}^{5}=\tfrac{1}{120}t(t^{2}-1)(t^{2}-4),
3.3.27c_{5}=\begin{cases}0.00488\ldots,&0<t<1,\\
0.00701\ldots,&-1<t<0\mbox{ or }1<t<2,\\
0.0234\ldots,&-2<t<-1\mbox{ or }2<t<3.\end{cases}

Seven-Point Formula

3.3.28f_{t}=\sum _{{k=-3}}^{3}A_{k}^{6}f_{k}+R_{{6,t}},\left|t\right|<3,
3.3.29
A_{{-3}}^{6}=\tfrac{1}{720}t(t^{2}-1)(t-3)(t^{2}-4),
A_{{-2}}^{6}=-\tfrac{1}{120}t(t^{2}-1)(t-2)(t^{2}-9),
A_{{-1}}^{6}=\tfrac{1}{48}t(t-1)(t^{2}-4)(t^{2}-9),
A_{0}^{6}=-\tfrac{1}{36}(t^{2}-1)(t^{2}-4)(t^{2}-9),
A_{1}^{6}=\tfrac{1}{48}t(t+1)(t^{2}-4)(t^{2}-9),
A_{2}^{6}=-\tfrac{1}{120}t(t^{2}-1)(t+2)(t^{2}-9),
A_{3}^{6}=\tfrac{1}{720}t(t^{2}-1)(t+3)(t^{2}-4),
3.3.30c_{6}=\begin{cases}0.00245\ldots,&\left|t\right|<1,\\
0.00459\ldots,&1<\left|t\right|<2,\\
0.0190\ldots,&2<\left|t\right|<3.\\
\end{cases}

Eight-Point Formula

3.3.31f_{t}=\sum _{{k=-3}}^{4}A_{k}^{7}f_{k}+R_{{7,t}},-3<t<4,
3.3.32
A_{{-3}}^{7}=-\tfrac{1}{5040}t(t^{2}-1)(t-3)(t-4)(t^{2}-4),
A_{{-2}}^{7}=\tfrac{1}{720}t(t^{2}-1)(t-2)(t-4)(t^{2}-9),
A_{{-1}}^{7}=-\tfrac{1}{240}t(t-1)(t-4)(t^{2}-4)(t^{2}-9),
A_{0}^{7}=\tfrac{1}{144}(t^{2}-1)(t-4)(t^{2}-4)(t^{2}-9),
A_{1}^{7}=-\tfrac{1}{144}t(t+1)(t-4)(t^{2}-4)(t^{2}-9),
A_{2}^{7}=\tfrac{1}{240}t(t^{2}-1)(t+2)(t-4)(t^{2}-9),
A_{3}^{7}=-\tfrac{1}{720}t(t^{2}-1)(t+3)(t-4)(t^{2}-4),
A_{4}^{7}=\tfrac{1}{5040}t(t^{2}-1)(t^{2}-4)(t^{2}-9),
3.3.33c_{7}=\begin{cases}0.00106\ldots,&0<t<1,\\
0.00139\ldots,&-1<t<0\mbox{ or }1<t<2,\\
0.00321\ldots,&-2<t<-1\mbox{ or }2<t<3,\\
0.0158\ldots,&-3<t<-2\mbox{ or }3<t<4.\\
\end{cases}

§3.3(iii) Divided Differences

The divided differences of f relative to a sequence of distinct points z_{0},z_{1},z_{2},\dots are defined by

3.3.34
f=f_{0},
f=({[z_{1}]}f-{[z_{0}]}f)/(z_{1}-z_{0}),
f=({[z_{1},z_{2}]}f-{[z_{0},z_{1}]}f)/(z_{2}-z_{0}),

and so on. Explicitly, the divided difference of order n is given by

3.3.35[z_{0},z_{1},\dots,z_{n}]f=\sum _{{k=0}}^{n}\left(\ifrac{f(z_{k})}{\prod _{{\substack{0\leq j\leq n\\
j\neq k}}}(z_{k}-z_{j})}\right).

If f and the z_{k} (=x_{k}) are real, and f is n times continuously differentiable on a closed interval containing the x_{k}, then

3.3.36[x_{0},x_{1},\dots,x_{n}]f=\frac{f^{{(n)}}(\xi)}{n!}

and again \xi is as in §3.3(i). If f is analytic in a simply-connected domain D, then for z\in{D},

3.3.37[z_{0},z_{1},\dots,z_{n}]f=\frac{1}{2\pi i}\int _{C}\frac{f(\zeta)}{\omega _{{n+1}}(\zeta)}d\zeta,

where \omega _{{n+1}}(\zeta) is given by (3.3.3), and C is a simple closed contour in {D} described in the positive rotational sense and enclosing z_{0},z_{1},\dots,z_{n}.

§3.3(iv) Newton’s Interpolation Formula

This represents the Lagrange interpolation polynomial in terms of divided differences:

3.3.38f(z)=[z_{0}]f+(z-z_{0})[z_{0},z_{1}]f+(z-z_{0})(z-z_{1})[z_{0},z_{1},z_{2}]f+\cdots+(z-z_{0})(z-z_{1})\cdots(z-z_{{n-1}})[z_{0},z_{1},\dots,z_{n}]f+R_{n}(z).

The interpolation error R_{n}(z) is as in §3.3(i). Newton’s formula has the advantage of allowing easy updating: incorporation of a new point z_{{n+1}} requires only addition of the term with [z_{0},z_{1},\dots,z_{{n+1}}]f to (3.3.38), plus the computation of this divided difference. Another advantage is its robustness with respect to confluence of the set of points z_{0},z_{1},\dots,z_{n}. For example, for k+1 coincident points the limiting form is given by [z_{0},z_{0},\dots,z_{0}]f=f^{{(k)}}(z_{0})/k!.

§3.3(v) Inverse Interpolation

In this method we interchange the roles of the points z_{k} and the function values f_{k}. It can be used for solving a nonlinear scalar equation f(z)=0 approximately. Another approach is to combine the methods of §3.8 with direct interpolation and §3.4.

Example

To compute the first negative zero a_{1}=-2.33810\; 7410\ldots of the Airy function f(x)=\mathop{\mathrm{Ai}\/}\nolimits\!\left(x\right)9.2). The inverse interpolation polynomial is given by

3.3.39x(f)=[f_{0}]x+(f-f_{0})[f_{0},f_{1}]x+(f-f_{0})(f-f_{1})[f_{0},f_{1},f_{2}]x;

compare (3.3.38). With x_{0}=-2.2, x_{1}=-2.3, x_{2}=-2.4, we obtain

3.3.40x=-2.2+1.44011\; 1973(f-0.09614\; 53780)+0.08865\; 85832\*(f-0.09614\; 53780)(f-0.02670\; 63331),

and with f=0 we find that x=-2.33823\; 2462, with 4 correct digits. By using this approximation to x as a new point, x_{3}=x, and evaluating [f_{0},f_{1},f_{2},f_{3}]x=1.12388\; 6190, we find that x=-2.33810\; 7409, with 9 correct digits.

For comparison, we use Newton’s interpolation formula (3.3.38)

3.3.41f(x)=0.09614\; 53780+0.69439\; 0 4495(x+2.1)-0.03007\; 14275(x+2.2)(x+2.3),

with the derivative

3.3.42f^{{\prime}}(x)=0.55906\; 90257-0.06014\; 28550x,

and compute an approximation to a_{1} by using Newton’s rule (§3.8(ii)) with starting value x=-2.5. This gives the new point x_{3}=-2.33934\; 0 514. Then by using x_{3} in Newton’s interpolation formula, evaluating [x_{0},x_{1},x_{2},x_{3}]f=-0.26608\; 28233 and recomputing f^{{\prime}}(x), another application of Newton’s rule with starting value x_{3} gives the approximation x=2.33810\; 7373, with 8 correct digits.

§3.3(vi) Other Interpolation Methods

For Hermite interpolation, trigonometric interpolation, spline interpolation, rational interpolation (by using continued fractions), interpolation based on Chebyshev points, and bivariate interpolation, see Bulirsch and Rutishauser (1968), Davis (1975, pp. 27–31), and Mason and Handscomb (2003, Chapter 6). These references also describe convergence properties of the interpolation formulas.

For interpolation of a bounded function f on \Real the cardinal function of f is defined by

3.3.43C(f,h)(x)=\sum _{{k=-\infty}}^{\infty}f(kh)S(k,h)(x),

where

3.3.44S(k,h)(x)=\frac{\mathop{\sin\/}\nolimits\!\left(\pi(x-kh)/h\right)}{\pi(x-kh)/h},

is called the Sinc function. For theory and applications see Stenger (1993, Chapter 3).