About the Project
3 Numerical MethodsAreas

§3.3 Interpolation


§3.3(i) Lagrange Interpolation

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

3.3.1 Pn(z)=k=0nk(z)fk=ωn+1(z)k=0nω~kz-zkfk=k=0nω~kfkz-zk/k=0nω~kz-zk,


3.3.2 k(z) =j=0nz-zjzk-zj,
k(zj) =δk,j.

Here the prime signifies that the factor for j=k is to be omitted, δk,j is the Kronecker symbol, and ωn+1(z) is the nodal polynomial

3.3.3 ωn+1(z)=k=0n(z-zk),

and the weights ω~k are

3.3.3_1 ω~k=1ωn+1(zk)=j=0n1zk-zj.

The final expression in (3.3.1) is the Barycentric form of the Lagrange interpolation formula. It is a direct consequence of the identity

3.3.3_2 1=k=0nk(z)=ωn+1(z)k=0nω~kz-zk,

and according to Berrut and Trefethen (2004) it is the most efficient representation of Pn(z).

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

3.3.4 f(z)=k=0nk(z)fk+Rn(z).

If f, x (=z), and the nodes xk are real, and f(n+1) is continuous on the smallest closed interval I containing x,x0,x1,,xn, then the error can be expressed

3.3.5 Rn(x)=f(n+1)(ξ)(n+1)!ωn+1(x),

for some ξI. If f is analytic in a simply-connected domain D1.13(i)), then for zD,

3.3.6 Rn(z)=ωn+1(z)2πiCf(ζ)(ζ-z)ωn+1(ζ)dζ,

where C is a simple closed contour in D described in the positive rotational sense and enclosing the points z,z1,z2,,zn.

§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.7 ft=f(x0+th)=k=n0n1Aknfk+Rn,t,

where the nodes xk=x0+kh (h>0) and function f are real,

3.3.8 n0 =-12(n-σ),
n1 =12(n+σ),
3.3.9 σ=12(1-(-1)n),

and Akn are the Lagrangian interpolation coefficients defined by

3.3.10 Akn=(-1)n1+k(k-n0)!(n1-k)!(t-k)m=n0n1(t-m).

The remainder is given by

3.3.11 Rn,t=Rn(x0+th)=hn+1(n+1)!f(n+1)(ξ)k=n0n1(t-k),

where ξ is as in §3.3(i).

Let cn be defined by

3.3.12 cn=1(n+1)!maxk=n0n1|t-k|,

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

3.3.13 |Rn,t|cnhn+1|f(n+1)(ξ)|.

Linear Interpolation

3.3.14 ft =(1-t)f0+tf1+R1,t,
3.3.15 c1 =18,

Three-Point Formula

3.3.16 ft=k=-11Ak2fk+R2,t,
3.3.17 A-12 =12t(t-1),
A02 =1-t2,
A12 =12t(t+1),
3.3.18 c2=1/(93)=0.0641,

Four-Point Formula

3.3.19 ft=k=-12Ak3fk+R3,t,
3.3.20 A-13 =-16t(t-1)(t-2),
A03 =12(t2-1)(t-2),
A13 =-12t(t+1)(t-2),
A23 =16t(t2-1),
3.3.21 c3={3128=0.0234,0<t<1,124=0.0416,-1<t<0 or 1<t<2.

Five-Point Formula

3.3.22 ft=k=-22Ak4fk+R4,t,
3.3.23 A-24 =124t(t2-1)(t-2),
A-14 =-16t(t-1)(t2-4),
A04 =14(t2-1)(t2-4),
A14 =-16t(t+1)(t2-4),
A24 =124t(t2-1)(t+2),
3.3.24 c4={0.0118,|t|<1,0.0302,1<|t|<2.

Six-Point Formula

3.3.25 ft=k=-23Ak5fk+R5,t,
3.3.26 A-25 =-1120t(t2-1)(t-2)(t-3),
A-15 =124t(t-1)(t2-4)(t-3),
A05 =-112(t2-1)(t2-4)(t-3),
A15 =112t(t+1)(t2-4)(t-3),
A25 =-124t(t2-1)(t+2)(t-3),
A35 =1120t(t2-1)(t2-4),
3.3.27 c5={0.00488,0<t<1,0.00701,-1<t<0 or 1<t<2,0.0234,-2<t<-1 or 2<t<3.

Seven-Point Formula

3.3.28 ft=k=-33Ak6fk+R6,t,
3.3.29 A-36 =1720t(t2-1)(t-3)(t2-4),
A-26 =-1120t(t2-1)(t-2)(t2-9),
A-16 =148t(t-1)(t2-4)(t2-9),
A06 =-136(t2-1)(t2-4)(t2-9),
A16 =148t(t+1)(t2-4)(t2-9),
A26 =-1120t(t2-1)(t+2)(t2-9),
A36 =1720t(t2-1)(t+3)(t2-4),
3.3.30 c6={0.00245,|t|<1,0.00459,1<|t|<2,0.0190,2<|t|<3.

Eight-Point Formula

3.3.31 ft=k=-34Ak7fk+R7,t,
3.3.32 A-37 =-15040t(t2-1)(t-3)(t-4)(t2-4),
A-27 =1720t(t2-1)(t-2)(t-4)(t2-9),
A-17 =-1240t(t-1)(t-4)(t2-4)(t2-9),
A07 =1144(t2-1)(t-4)(t2-4)(t2-9),
A17 =-1144t(t+1)(t-4)(t2-4)(t2-9),
A27 =1240t(t2-1)(t+2)(t-4)(t2-9),
A37 =-1720t(t2-1)(t+3)(t-4)(t2-4),
A47 =15040t(t2-1)(t2-4)(t2-9),
3.3.33 c7={0.00106,0<t<1,0.00139,-1<t<0 or 1<t<2,0.00321,-2<t<-1 or 2<t<3,0.0158,-3<t<-2 or 3<t<4.

§3.3(iii) Divided Differences

The divided differences of f relative to a sequence of distinct points z0,z1,z2, are defined by

3.3.34 [z0]f =f0,
[z0,z1]f =([z1]f-[z0]f)/(z1-z0),
[z0,z1,z2]f =([z1,z2]f-[z0,z1]f)/(z2-z0),

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

3.3.35 [z0,z1,,zn]f=k=0n(f(zk)/0jnjk(zk-zj)).

If f and the zk (=xk) are real, and f is n times continuously differentiable on a closed interval containing the xk, then

3.3.36 [x0,x1,,xn]f=f(n)(ξ)n!

and again ξ is as in §3.3(i). If f is analytic in a simply-connected domain D, then for zD,

3.3.37 [z0,z1,,zn]f=12πiCf(ζ)ωn+1(ζ)dζ,

where ωn+1(ζ) is given by (3.3.3), and C is a simple closed contour in D described in the positive rotational sense and enclosing z0,z1,,zn.

§3.3(iv) Newton’s Interpolation Formula

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

3.3.38 f(z)=[z0]f+(z-z0)[z0,z1]f+(z-z0)(z-z1)[z0,z1,z2]f++(z-z0)(z-z1)(z-zn-1)[z0,z1,,zn]f+Rn(z).

The interpolation error Rn(z) is as in §3.3(i). Newton’s formula has the advantage of allowing easy updating: incorporation of a new point zn+1 requires only addition of the term with [z0,z1,,zn+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 z0,z1,,zn. For example, for k+1 coincident points the limiting form is given by [z0,z0,,z0]f=f(k)(z0)/k!.

§3.3(v) Inverse Interpolation

In this method we interchange the roles of the points zk and the function values fk. 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.


To compute the first negative zero a1=-2.33810 7410 of the Airy function f(x)=Ai(x)9.2). The inverse interpolation polynomial is given by

3.3.39 x(f)=[f0]x+(f-f0)[f0,f1]x+(f-f0)(f-f1)[f0,f1,f2]x;

compare (3.3.38). With x0=-2.2, x1=-2.3, x2=-2.4, we obtain

3.3.40 x=-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, x3=x, and evaluating [f0,f1,f2,f3]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.41 f(x)=0.09614 53780+0.69439 04495(x+2.1)-0.03007 14275(x+2.2)(x+2.3),

with the derivative

3.3.42 f(x)=0.55906 90257-0.06014 28550x,

and compute an approximation to a1 by using Newton’s rule (§3.8(ii)) with starting value x=-2.5. This gives the new point x3=-2.33934 0514. Then by using x3 in Newton’s interpolation formula, evaluating [x0,x1,x2,x3]f=-0.26608 28233 and recomputing f(x), another application of Newton’s rule with starting value x3 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 the cardinal function of f is defined by

3.3.43 C(f,h)(x)=k=-f(kh)S(k,h)(x),


3.3.44 S(k,h)(x)=sin(π(x-kh)/h)π(x-kh)/h,

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