Lecture 15: Function and Derivative Interpolation

1.Interpolation in function and derivative values - Hermite interpolation

In addition to sampling of the function f:ℝ→ℝ, information on the derivatives of f might also be available, as in the data set

𝒟'={(xi,yi=f(xi),yi'=f'(xi)),i=0,1,…,n}. (1)

The extended data set can again be interpolated by a polynomial, this time of degree 2n+1 given in the monomial, Lagrange or Newton form.

Monomial form of interpolating polynomial.
Using the monomial basis

ℳ2n+1(t)=[ 1 t t2 t3 … t2n+1 ],

the interpolating polynomial is

p(t)=ℳ2n+1(t)⁡𝒂=[ 1 t … t2n+1 ][ 0 a1 ⋮ a2n+1 ]=a0+a1t+⋯+a2n+1⁡t2n+1,

with derivative

p'(t)=a1+2a2t+⋯+(2n+1)a2n+1t2n.

The above suggests constructing a basis set of monomials and their derivatives

ℳ2n+1'(t)=[ 1 t t2 t3 … t2n+1 0 1 2t 3t2 … (2n+1)⁡t2n ],

to allow setting the function p(xi)=yi, and derivative conditions p(xi)=yi'. The columns of ℳ2n+1'(t) are linearly independent since

α[ tj j⁡tj-1 ]+β[ tk k⁡tk-1 ]=0,

can only be satisfied for all t by α=β=0.

∘Sampling at 𝒙∈ℝ(n+1) gives 𝑴=ℳ2n+1'(𝒙)∈ℝ(2n+2)×(2n+2), e.g., for n=2, the matrix is

𝑴=( 1 x0 x02 x03 x04 x05 1 x1 x12 x13 x14 x15 1 x2 x22 x23 x24 x25 0 1 2⁢x0 3⁢x02 4⁢x03 5⁢x04 0 1 2⁢x1 3⁢x12 4⁢x13 5⁢x14 0 1 2⁢x2 3⁢x22 4⁢x23 5⁢x24 ),

is obtained.

For general n, 𝑴 is of full rank for distinct sample points with a determinant reminiscent of that of the Vandermonde matrix

det(𝑴)=∏0⩽i<j⩽n(xi-xj)4.

The interpolation conditions lead to the linear system

𝑴⁡𝒂=[ 𝒚 𝒚' ],

whose solution requires 𝒪([2(n+1)]3/3) operations. An error formula is again obtained by repeated application of Rolle's theorem, i.e., for p interpolant of data set 𝒟', ∃ξt∈(x0,xn) such that

f(t)-p(t)=f(2n+2)(ξt)(2n+2)!⁡∏j=0n(t-xj)2.

The above approach generalizes to higher-order derivatives, e.g., for

𝒟''={(xi,yi=f(xi),yi'=f'(xi),yi''=f''(xi)),i=0,1,…,n}, (2)

the basis set is

ℳ3n+2''(t)=[ 1 t t2 t3 … t3n+2 0 1 2t 3t2 … (3n+2)⁡t3n+1 0 0 2 6t … (3n+2)(3n+1)t3n ],

with interpolant

p(t)=ℳ⁡3n+2''(t)𝒂,

with 𝒂∈ℝ3(n+1) determined by solving

𝑴⁡𝒂=[ 𝒚 𝒚' 𝒚'' ]

with 𝑴=ℳ3n+2''(𝒙)∈ℝ(3n+3)×(3n+3), and error formula

f(t)-p(t)=f(3n+3)(ξt)(3n+3)!⁡∏j=0n(t-xj)3.

Lagrange form of interpolating polynomial.
As in the function value interpolation case, a basis set that evaluates to an identity matrix when sampled at 𝒙∈ℝn+1 is obtained by L⁡U-factorization of the sampled monomial matrix

ℳ2n+1'(𝒙)=𝑴=𝑳⁡⁡𝑼=𝑰⁡𝑳⁡𝑼=ℒ2n+1'(𝒙)⁡𝑳⁡𝑼,

that for arbitrary t leads to the basis set

ℒ2n+1'(t)=ℳ2n+1'(t)⁡𝑼-1⁡𝑳-1=[ a0(t) a1(t) … an(t) b0(t) b1(t) … bn(t) a0'(t) a1'(t) … an'(t) b0'(t) b1'(t) … bn'(t) ].

The interpolating polynomial of data set 𝒟'={(xi,yi=f(xi),yi'=f'(xi)),i=0,1,…,n} is

p(t)=∑i=0nyi⁡ai(t)+∑i=0nyi'⁡bi(t),

where the basis functions can be expressed in terms of the Lagrange polynomials

ℓi(t)=∏j=0n't-xjxi-xj,

as

ai(t)=[1-2(t-xi)ℓi'(xi)]ℓi2(t),⁡⁡⁡bi(t)=(t-xi)ℓi2(t),

and have the properties

ai(xj)=δi⁡j,⁡⁡ai'(xj)=0,⁡⁡bi(xj)=0,⁡⁡bi'(xj)=δi⁡j⁡.

∘As, an example, consider the L⁡U-factorization of matrix 𝑴=ℳ2n+1'(𝒙) for n=1

( 1 0 0 0 1 1 0 0 0 -1x0-x1 1 0 0 -1x0-x1 -1 1 )=( 1 0 0 0 1 1 0 0 0 1x1-x0 1 0 0 1x1-x0 -1 1 )( 1 x0 x02 x03 0 x1-x0 x12-x02 x13-x03 0 0 x0-x1 2⁢x02-x1⁢x0-x12 0 0 0 (x0-x1)2 ).

∘The factor inverses are

𝑳-1=( 1 0 0 0 -1 1 0 0 1x1-x0 1x0-x1 1 0 2x1-x0 2x0-x1 1 1 ),⁡𝑼-1=( 1 x0x0-x1 x0⁢x1x0-x1 -x02⁢x1(x0-x1)2 0 1x1-x0 -x0+x1x0-x1 x0⁢(x0+2⁢x1)(x0-x1)2 0 0 1x0-x1 -2⁢x0+x1(x0-x1)2 0 0 0 1(x0-x1)2 )

∘The functions that result

{[1-2(t-x0)1x0-x1](t-x1⁢x0-x1)2,[1-2(t-x1)1x1-x0](t-x0⁢x1-x0)2,(t-x0)(t-x1⁢x0-x1)2,(t-x1)(t-x0⁢x1-x0)2},

are indeed expressed in terms of ℓi(t) as

{[1-2(t-x0)ℓ0'(x0)]ℓ02(t),[1-2(t-x1)ℓ1'(x1)]ℓ02(t),(t-x0)ℓ02(t),(t-x1)ℓ12(t)}.

The procedure can be extended to derivatives of arbitrary order, e.g., the data set 𝒟'' is interpolated by

p(t)=∑i=0nyi⁡ai(t)+∑i=0nyi'⁡bi(t)+∑i=0nyi''⁡ci(t),

where the Lagrange basis polynomials are given as

ℒ3n+2''(t)=ℳ3n+2''(t)⁡𝑼-1⁡𝑳-1=[ a0(t) … an(t) b0(t) … bn(t) c0(t) … cn(t) a0'(t) … an'(t) b0'(t) … bn'(t) c0'(t) … cn'(t) a0''(t) … an''(t) b0''(t) … bn''(t) c0''(t) … cn''(t) ].

Newton form of interpolating polynomial.
As before, inverting only one factor of the ℳ2n+1'(t)=ℒ2n+1'(t)⁡𝑳⁡𝑼 mapping yields a triangular basis set 𝒮'(t)=[ s0(t) s1(t) s2(t) … ]

ℳ2n+1'(t)𝑼-1=𝒮2n+1'(t)⁡.

∘The first six basis polynomials obtained for n=2 are

{1,t-x0x1-x0,(t-x0)⁢(t-x1)(x2-x0)(x2-x1),(t-x0)⁢(t-x1)⁢(t-x2)(x2-x0)⁢(x1-x0),(t-x0)2⁢(t-x1)⁢(t-x2)(x1-x0)2⁢(x1-x2),(t-x0)2⁢(t-x1)2⁢(t-x2)(x2-x0)2⁢(x2-x1)2}.

∘A closer link to divided difference and differential calculus is obtained by permuting rows of 𝑴, e.g., for n=2

𝑷⁡𝑴=( 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 )( 1 x0 x02 x03 x04 x05 1 x1 x12 x13 x14 x15 1 x2 x22 x23 x24 x25 0 1 2⁢x0 3⁢x02 4⁢x03 5⁢x04 0 1 2⁢x1 3⁢x12 4⁢x13 5⁢x14 0 1 2⁢x2 3⁢x22 4⁢x23 5⁢x24 )=( 1 x0 x02 x03 x04 x05 0 1 2⁢x0 3⁢x02 4⁢x03 5⁢x04 1 x1 x12 x13 x14 x15 0 1 2⁢x1 3⁢x12 4⁢x13 5⁢x14 1 x2 x22 x23 x24 x25 0 1 2⁢x2 3⁢x22 4⁢x23 5⁢x24 ).

The first six basis polynomials thus obtained are

{1,⁡t-x0,⁡⁡(t-x0)2(x1-x0)2,⁡⁡(t-x0)2⁢(t-x1)(x1-x0)2,⁡⁡(t-x0)​2⁢(t-x1)2(x2-x0)2⁢(x2-x1)2,(t-x0)2⁢(t-x1)2⁢(t-x2)(x2-x0)2⁢(x2-x1)2}.

and upon rescaling generalize to the basis set

𝒩2n+1'(t)=[ n0(t) n1(t) … n2n+1(t) ],

with

n2k(t)=∏j=0k-1(t-xj)2,n2k+1(t)=(t-xk)n2k(t),k=0,1,…,n

known as the Newton basis set with repetitions.

The interpolating polynomial in Newton divided difference form is

p(t)=[y0]+[y0,y0](t-x0)+[y1,y0,y0](t-x0)2+⋯+[yn,yn,…,y0,y0](t-x0)2…(t-xn-1)2(t-xn).

Divided difference with repeated values are replaced by their, limits, i.e., the derivatives

[yk,yk]=limxk-1→xkyk-yk-1xk-xk-1=yk'⁡.

The forward substitution can again be organized in a table.

i xi [yi] [yi,yi-1] [yi,yi-1,yi-2] 0 x0 y0 - - 0 x0 y0 y0' 1 x1 y1 (-y0') 1 x1 y1 y1' (y1'-) 2 x2 y2 (-y1') ⋱ 2 x2 y2 y2' ⋮ ⋮ ⋮ ⋮ n xn yn yn' (yn'-) …

Table 1. Table of repeated divided differences. The Newton basis coefficients are the diagonal terms.

Interpolation of data containing higher derivatives, or differing orders of derivative information at each node are poissible. For multiple repeated values arising in the limit xi+k→xi of sample points xi⩽xi+1⩽⋯⩽xi+k the divided difference is determined by

[yi+k,yi+k-1,…,yi]={ xi<xi+k yi(k)k! xi=xi+k ..