Interpolation and Polynomial Approximation

師大

Interpolation and Polynomial Approximation

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

November 18, 2007

師大

Outline

1 Interpolation and the Lagrange polynomial

2 Divided Differences

3 Hermite Interpolation

4 Cubic spline interpolation

師大

Introduction

Question

In reviewing these data, how to provide a reasonable estimate of the population, say, in 1965, or even in 2010.

師大

Interpolation

Suppose we do not know the function f , but a few information (data) about f , now we try to compute a function g that approximates f . Theorem (Weierstrass Approximation Theorem)

Suppose that f is defined and continuous on [a, b]. For any ε > 0, there exists a polynomial P (x), defined on [a, b], with the property that

|f (x) − P (x)| < ε, for all x in [a, b].

Reason for using polynomial

1 They uniformly approximate continuous function (Weierstrass Theorem)

2 The derivatives and indefinite integral of a polynomial are easy to determine and are also polynomials.

師大

Interpolation and the Lagrange polynomial

Property

The linear function passing through (x₀, f (x₀)) and (x1, f (x1)) is unique.

Let

L₀(x) = x − x₁

x0− x₁, L₁(x) = x − x₀ x1− x₀ and

P (x) = L₀(x)f (x₀) + L₁(x)f (x₁).

Then

P (x0) = f (x0), P (x1) = f (x1).

師大

Question

How to find the polynomial of degree n that passes through

(x₀, f (x₀)), . . . , (x_n, f (x_n))?

Theorem

If (xi, yi), xi, yi ∈ R, i = 0, 1, . . . , n, are n + 1 distinct pairs of data point, then there is a unique polynomial Pn of degree at most n such that

Pn(xi) = yi, (0 ≤ i ≤ n). (1)

師大

Proof of Existence (by mathematical induction)

1 The theorem clearly holds for n = 0 (only one data point (x₀, y₀)) since one may choose the constant polynomial P₀(x) = y₀ for all x.

2 Assume that the theorem holds for n ≤ k, that is, there is a polynomial Pk, deg(Pk) ≤ k, such that yi= Pk(xi), for 0 ≤ i ≤ k.

3 Next we try to construct a polynomial of degree at most k + 1 to interpolate (xi, yi), 0 ≤ i ≤ k + 1. Let

P_k+1(x) = P_k(x) + c(x − x₀)(x − x₁) · · · (x − x_k), where

c = y_k+1− P_k(x_k+1)

(x_k+1− x₀)(x_k+1− x₁) · · · (x_k+1− x_k). Since xi are distinct, the polynomial P_k+1(x) is well-defined and deg(P_k+1) ≤ k + 1. It is easy to verify that

P_k+1(xi) = yi, 0 ≤ i ≤ k + 1.

師大

Proof of Uniqueness

Suppose there are two such polynomials Pn and Qn satisfying (1). Define S_n(x) = P_n(x) − Q_n(x).

Since both deg(P_n) ≤ n and deg(Q_n) ≤ n, deg(S_n) ≤ n. Moreover Sn(xi) = Pn(xi) − Qn(xi) = yi− y_i= 0,

for 0 ≤ i ≤ n. This means that Sn has at least n + 1 zeros, it therefore must be S_n= 0. Hence P_n= Q_n.

師大

Idea

Construct polynomial P (x) with deg(P ) ≤ n as

P (x) = Ln,0(x)f (x0) + · · · + Ln,n(x)f (xn), where

1 L_n,k(x) are polynomial with degree n for 0 ≤ k ≤ n.

2 L_n,k(x_k) = 1 and L_n,k(xi) = 0 for i 6= k.

Then

P (xk) = f (xk) for k = 0, 1, . . . , n.

1 deg(L_n,k) = n and L_n,k(x_i) = 0 for i 6= k:

L_n,k(x) = c_k(x − x₀)(x − x₁) · · · (x − x_k−1)(x − x_k+1) · · · (x − x_n)

2 L_n,k(x_k) = 1:

Ln,k(x) = (x − x0)(x − x1) · · · (x − xk−1)(x − xk+1) · · · (x − xn) (x_k− x₀)(x_k− x₁) · · · (x_k− x_k−1)(x_k− x_k+1) · · · (x_k− x_n)

師大

Theorem

If x0, x1, . . . , xn are n + 1 distinct numbers and f is a function whose values are given at these numbers, then a unique polynomial of degree at most n exists with

f (x_k) = P (x_k), for k = 0, 1, . . . , n.

This polynomial is given by P (x) =

n

X

k=0

f (x_k)L_n,k(x)

which is called the nth Lagrange interpolating polynomial.

Note that we will write L_n,k(x) simply as L_k(x) when there is no confusion as to its degree.

師大

Example

Given the following 4 data points,

x_i 0 1 3 5

y_i 1 2 6 7

find a polynomial in Lagrange form to interpolate these data.

Solution: The interpolating polynomial in the Lagrange form is P3(x) = L0(x) + 2L1(x) + 6L2(x) + 7L3(x) with L0(x) = (x − 1)(x − 3)(x − 5)

(0 − 1)(0 − 3)(0 − 5) = − 1

15(x − 1)(x − 3)(x − 5), L₁(x) = (x − 0)(x − 3)(x − 5)

(1 − 0)(1 − 3)(1 − 5) = 1

8x(x − 3)(x − 5), L₂(x) = (x − 0)(x − 1)(x − 5)

(3 − 0)(3 − 1)(3 − 5) = − 1

12x(x − 1)(x − 5), L3(x) = (x − 0)(x − 1)(x − 3)

(5 − 0)(5 − 1)(5 − 3) = 1

40x(x − 1)(x − 3).

師大

Question

What’s the error involved in approximating f (x) by the interpolating polynomial P (x)?

Theorem Suppose

1 x₀, . . . , x_n are distinct numbers in [a, b],

2 f ∈ Cⁿ⁺¹[a, b].

Then, ∀x ∈ [a, b], ∃ ξ(x) ∈ (a, b) such that f (x) = P (x) +fⁿ⁺¹(ξ(x))

(n + 1)! (x − x0)(x − x1) · · · (x − xn). (2) Proof: If x = x_k, for any k = 0, 1, . . . , n, then f (x_k) = P (x_k) and (2) is satisfied.

師大

If x 6= x_k, for all k = 0, 1, . . . , n, define

g(t) = f (t) − P (t) − [f (x) − P (x)]

n

Y

i=0

t − xi

x − x_i.

Since f ∈ Cⁿ⁺¹[a, b] and P ∈ C^∞[a, b], it follows that g ∈ Cⁿ⁺¹[a, b].

Since

g(x_k) = [f (x_k) − P (x_k)] − [f (x) − P (x)]

n

Y

i=0

x_k− x_i x − x_i = 0, and

g(x) = [f (x) − P (x)] − [f (x) − P (x)]

n

Y

i=0

x − xi

x − x_i = 0,

it implies that g is zero at x, x₀, x₁, . . . , x_n. By the Generalized Rolle’s Theorem, ∃ ξ ∈ (a, b) such that gⁿ⁺¹(ξ) = 0. That is

師大

0 = g⁽ⁿ⁺¹⁾(ξ) (3)

= f⁽ⁿ⁺¹⁾(ξ) − P⁽ⁿ⁺¹⁾(ξ) − [f (x) − P (x)] dⁿ⁺¹ dtⁿ⁺¹

" _n Y

i=0

(t − xi) (x − xi)

#

t=ξ

.

Since deg(P ) ≤ n, it implies that P⁽ⁿ⁺¹⁾(ξ) = 0. On the other hand, Qn

i=0[(t − x_i)/(x − x_i)] is a polynomial of degree (n + 1), so

n

Y

i=0

(t − x_i) (x − x_i) =

 1

Qn

i=0(x − x_i)



tⁿ⁺¹+ (lower-degree terms in t),

and

dⁿ⁺¹ dtⁿ⁺¹

" _n Y

i=0

(t − x_i) (x − xi)

#

= (n + 1)!

Qn

i=0(x − xi).

師大

Equation (3) becomes

0 = f⁽ⁿ⁺¹⁾(ξ) − [f (x) − P (x)] (n + 1)!

Qn

i=0(x − xi), i.e.,

f (x) = P (x) + f⁽ⁿ⁺¹⁾(ξ) (n + 1)!

n

Y

i=0

(x − xi).

Example

1 Goal: Prepare a table for the function f (x) = e^x for x ∈ [0, 1].

2 x_j+1− x_j = h for j = 0, 1, . . . , n − 1.

3 What should h be for linear interpolation to give an absolute error of at most 10⁻⁶?

師大

Suppose x ∈ [x_j, x_j+1]. Equation (2) implies that

|f (x) − P (x)| =     

f⁽²⁾(ξ)

2! (x − x_j)(x − x_j+1)     

= |f⁽²⁾(ξ)|

2! |x − x_j||x − x_j+1|

= e^ξ

2|x − jh||x − (j + 1)h|

≤ 1

2

 max

ξ∈[0,1]e^ξ

 

xj≤x≤xmax_j+1|x − jh||x − (j + 1)h|



= e

2



xj≤x≤xmax_j+1|x − jh||x − (j + 1)h|

 .

師大

Let

g(x) = (x − jh)(x − (j + 1)h), for jh ≤ x ≤ (j + 1)h.

Then

xj≤x≤xmax_j+1|g(x)| =     g



j + 1 2

 h

   

= h² 4 . Consequently,

|f (x) − P (x)| ≤ eh²

8 ≤ 10⁻⁶, which implies that

h < 1.72 × 10⁻³.

Since n = (1 − 0)/h must be an integer, one logical choice for the step size is h = 0.001.

師大

Difficulty for the Lagrange interpolation

1 If more data points are added to the interpolation problem, all the cardinal functions Lk have to be recalculated.

2 We shall now derive the interpolating polynomials in a manner that uses the previous calculations to greater advantage.

Definition

1 f is a function defined at x₀, x₁, . . . , x_n

2 Suppose that m1, m2, . . . , m_k are k distinct integers with 0 ≤ mi ≤ n for each i.

The Lagrange polynomial that interpolates f at the k points xm1, xm2, . . . , xmk is denoted Pm1,m2,...,mk(x).

師大

Theorem

Let f be defined at distinct points x₀, x₁, . . . , x_k, and 0 ≤ i, j ≤ k, i 6= j.

Then

P (x) = (x − xj)

(x_i− x_j)P0,1,...,j−1,j+1,...,k(x) − (x − xi)

(x_i− x_j)P0,1,...,i−1,i+1,...,k(x) describes the k-th Lagrange polynomial that interpolates f at the k + 1 points x0, x1, . . . , xk.

Proof: Since

deg(P0,1,...,j−1,j+1,...,k) ≤ k − 1 and

deg(P0,1,...,i−1,i+1,...,k) ≤ k − 1,

it implies that deg(P ) ≤ k. If 0 ≤ r ≤ k and r 6= i, j, then

師大

P (xr) = (xr− x_j)

(xi− x_j)P0,1,...,j−1,j+1,...,k(xr) −(xr− x_i)

(xi− x_j)P0,1,...,i−1,i+1,...,k(xr)

= (xr− x_j)

(x_i− x_j)f (x_r) − (x_r− x_i)

(x_i− x_j)f (x_r) = f (x_r).

Moreover

P (x_i) = (x_i− x_j)

(xi− x_j)P0,1,...,j−1,j+1,...,k(x_i) − (x_i− x_i)

(xi− x_j)P0,1,...,i−1,i+1,...,k(x_i)

= f (xi) and

P (xj) = (x_j− x_j)

(xi− x_j)P0,1,...,j−1,j+1,...,k(xj) −(x_j− x_i)

(xi− x_j)P0,1,...,i−1,i+1,...,k(xj)

= f (xj).

Therefore P (x) agrees with f at all points x0, x1, . . . , x_k. By the uniqueness theorem, P (x) is the k-th Lagrange polynomial that interpolates f at the k + 1 points x0, x1, . . . , xk, i.e., P ≡ P0,1,...,k.

師大

Neville’s method

The theorem implies that the Lagrange interpolating polynomial can be generated recursively. The procedure is called the Neville’s method.

1 Denote

Q_i,j = Pi−j,i−j+1,...,i−1,i.

2 Hence Qi,j, 0 ≤ j ≤ i, denotes the interpolating polynomial of degree j on the j + 1 points xi−j, xi−j+1, . . . , xi−1, xi.

3 The polynomials can be computed in a manner as shown in the following table.

x₀ P₀= Q_0,0

x₁ P₁= Q_1,0 P_0,1 = Q_1,1

x₂ P₂= Q_2,0 P_1,2 = Q_2,1 P_0,1,2= Q_2,2

x₃ P₃= Q_3,0 P_2,3 = Q_3,1 P_1,2,3= Q_3,2 P_0,1,2,3 = Q_3,3

x₄ P₄= Q_4,0 P_3,4 = Q_4,1 P_2,3,4= Q_4,2 P_1,2,3,4 = Q_4,3 P_0,1,2,3,4= Q_4,4

師大

x0 P0= Q0,0

x1 P1= Q1,0 P0,1 = Q1,1

x2 P2= Q2,0 P1,2 = Q2,1 P0,1,2= Q2,2

x3 P3= Q3,0 P2,3 = Q3,1 P1,2,3= Q3,2 P0,1,2,3 = Q3,3

x4 P4= Q4,0 P3,4 = Q4,1 P2,3,4= Q4,2 P1,2,3,4 = Q4,3 P0,1,2,3,4= Q4,4

P0,1(x) = (x − x1)P0(x) − (x − x0)P1(x)

x₀− x₁ ,

P1,2(x) = (x − x2)P1(x) − (x − x1)P2(x)

x₁− x₂ ,

...

P_0,1,2(x) = (x − x₂)P_0,1(x) − (x − x₀)P_1,2(x) x0− x₂

...

師大

Example

Compute approximate value of f (1.5) by using the following datas:

x f (x)

1.0 0.7651977 1.3 0.6200860 1.6 0.4554022 1.9 0.2818186 2.2 0.1103623

師大 1 The first-degree approximation:

Q1,1(1.5) = (x − x0)Q1,0− (x − x₁)Q0,0

x₁− x₀

= (1.5 − 1.0)Q1,0− (1.5 − 1.3)Q_0,0

1.3 − 1.0 = 0.5233449,

Q_2,1(1.5) = (x − x1)Q2,0− (x − x₂)Q1,0

x₂− x₁

= (1.5 − 1.3)Q2,0− (1.5 − 1.6)Q_1,0

1.6 − 1.3 = 0.5102968,

Q_3,1(1.5) = (x − x2)Q3,0− (x − x₃)Q2,0

x₃− x₂

= (1.5 − 1.6)Q_3,0− (1.5 − 1.9)Q_2,0

1.9 − 1.6 = 0.5132634,

Q_4,1(1.5) = (x − x₃)Q_4,0− (x − x₄)Q_3,0 x₄− x₃

= (1.5 − 1.9)Q_4,0− (1.5 − 2.2)Q_3,0

2.2 − 1.9 = 0.5104270.

師大 1 The second-degree approximation:

Q2,2(1.5) = (x − x₁)Q_2,1− (x − x₂)Q_1,1 x2− x₁

= (1.5 − 1.3)Q_2,1− (1.5 − 1.6)Q_1,1

1.6 − 1.3 = 0.5124715,

Q3,2(1.5) = (x − x₂)Q_3,1− (x − x₃)Q_2,1 x3− x₂

= (1.5 − 1.6)Q3,1− (1.5 − 1.9)Q_2,1

1.9 − 1.6 = 0.5112857,

Q4,2(1.5) = (x − x3)Q4,1− (x − x₄)Q3,1

x4− x₃

= (1.5 − 1.9)Q4,1− (1.5 − 2.2)Q_3,1

2.2 − 1.9 = 0.5137361.

師大

x f (x) 1st-deg 2nd-deg 3rd-deg 4th-deg 1.0 0.7651977

1.3 0.6200860 0.5233449

1.6 0.4554022 0.5102968 0.5124715

1.9 0.2818186 0.5132634 0.5112857 0.5118127

2.2 0.1103623 0.5104270 0.5137361 0.5118302 0.5118200 Table: Results of the higher-degree approximations

x f (x) 1st-deg 2nd-deg 3rd-deg 4th-deg 5th-deg

1.0 0.7651977

1.3 0.6200860 0.5233449

1.6 0.4554022 0.5102968 0.5124715

1.9 0.2818186 0.5132634 0.5112857 0.5118127

2.2 0.1103623 0.5104270 0.5137361 0.5118302 0.5118200

2.5 −0.0483838 0.4807699 0.5301984 0.5119070 0.5118430 0.5118277 Table: Results of adding (x₅, f (x₅))

師大

Divided Differences

Neville’s method: generate successively higher-degree polynomial approximations at a specific point.

Divided-difference methods: successively generate the polynomials themselves.

Suppose that Pn(x) is the nth Lagrange polynomial that agrees with f at distinct x₀, x₁, . . . , x_n. Express P_n(x) in the form

Pn(x) = a0+ a1(x − x0) + a2(x − x0)(x − x1) + · · · +a_n(x − x₀)(x − x₁) · · · (x − x_n−1).

Determine constant a0:

a₀ = P_n(x₀) = f (x₀) Determine a₁:

f (x0) + a1(x1− x₀) = Pn(x1) = f (x1) ⇒ a₁ = f (x1) − f (x0) x1− x₀ .

師大

The zero divided difference of the function f with respect to xi: f [xi] = f (xi)

The first divided difference of f with respect to x_i and x_i+1 is denoted and defined as

f [xi, xi+1] = f [xi+1] − f [xi] x_i+1− x_i . The second divided difference of f is defined as

f [xi, xi+1, xi+2] = f [xi+1, xi+2] − f [xi, xi+1] xi+2− x_i .

The k-th divided difference relative to xi, xi+1, . . . , x_i+k is given by f [x_i, x_i+1, . . . , x_i+k] = f [xi+1, xi+2, . . . , x_i+k] − f [xi, xi+1, . . . , x_i+k−1]

x_i+k− x_i .

師大

As might be expected from the evaluation of a0, a1, . . . , an in Pn(x), the required coefficients are given by

a_k= f [x₀, x₁, . . . , x_k], k = 0, 1, . . . , n.

Therefore Pn(x) can be expressed as P_n(x) = f [x₀] +

n

X

k=1

f [x₀, x₁, . . . , x_k](x − x₀)(x − x₁) · · · (x − x_k−1).

This formula is known as the Newton’s divided-difference formula.

xi k = 0 k = 1 k = 2 k = 3

x0 f [x0]

> f [x₀, x₁]

x₁ f [x₁] > f [x₀, x₁, x₂]

> f [x1, x2] > f [x0, x1, x2, x3] x₂ f [x₂] > f [x₁, x₂, x₃]

> f [x₂, x₃] x3 f [x3]

Table: Outline for the generation of the divided differences

師大

Example

Given the following 4 points (n = 3)

x_i 0 1 3 5

y_i 1 2 6 7

find a polynomial of degree 3 in Newton’s form to interpolate these data.

Solution:

xi k = 0 k = 1 k = 2 k = 3

0 1

> 1

1 2 > ¹₃

> 2 > −₁₂₀¹⁷

3 6 > −³₈

> ¹₂

5 7

So,

P (x) = 1 + x +1

3x(x − 1) − 17

120x(x − 1)(x − 3).

師大

Theorem

Suppose f ∈ Cⁿ[a, b] and x0, x1, . . . , xn are distinct numbers in [a, b].

Then there exists ξ ∈ (a, b) such that

f [x₀, x₁, . . . , x_n] = f⁽ⁿ⁾(ξ) n! . Proof: Define

g(x) = f (x) − Pn(x).

Since P_n(x_i) = f (x_i) for i = 0, 1, . . . , n, g has n + 1 distinct zeros in [a, b]. By the generalized Rolle’s Theorem, ∃ ξ ∈ (a, b) such that

0 = g⁽ⁿ⁾(ξ) = f⁽ⁿ⁾(ξ) − P_n⁽ⁿ⁾(ξ).

Note that

P_n⁽ⁿ⁾(x) = n!f [x0, x1, . . . , xn].

As a consequence

f [x0, x1, . . . , xn] = f⁽ⁿ⁾(ξ) n! .

師大

Let

h = x_n− x₀

n = x_i+1− x_i, i = 0, 1, . . . , n − 1.

Then each x_i= x₀+ ih, i = 0, 1, . . . , n. For any x ∈ [a, b], write x = x0+ sh, s ∈ R.

Hence x − xi = (s − i)h and Pn(x) = f [x0] +

n

X

k=1

f [x0, x1, . . . , x_k](x − x0)(x − x1) · · · (x − x_k−1)

= f (x0) +

n

X

k=1

f [x0, x1, . . . , xk](s − 0)h(s − 1)h · · · (s − k + 1)h

= f (x0) +

n

X

k=1

f [x0, x1, . . . , xk]s(s − 1) · · · (s − k + 1)h^k

= f (x₀) +

n

X

k=1

f [x₀, x₁, . . . , x_k]k!s k



h^k, (4)

師大

The binomial formula

s k



= s(s − 1) · · · (s − k + 1) k!

The forward difference notation 4

4f (x_i) = f (xi+1) − f (xi) and

4^kf (xi) = 4^k−1f (xi+1) − 4^k−1f (xi) = 4



4^k−1f (xi)

 , for i = 0, 1, . . . , n − 1.

With this notation,

f [x0, x1] = f (x1) − f (x0) x₁− x₀ = 1

h4f (x₀) f [x0, x1, x2] = f [x1, x2] − f [x0, x1]

x₂− x₀

=

1

h4f (x₁) −_h¹4f (x₀)

2h = 1

2h²4²f (x0),

師大

In general

f [x0, x1, . . . , x_k] = 1

k! h^k4^kf (x0).

The Newton forward-difference formula:

Pn(x) = f (x0) +

n

X

k=1

s k



4^kf (x0).

The interpolation nodes are arranged as xn, xn−1, . . . , x0: P_n(x) = f [x_n] + f [x_n, x_n−1](x − x_n)

+f [x_n, x_n−1, x_n−2](x − x_n)(x − x_n−1) + · · · +f [x_n, . . . , x₀](x − x_n)(x − x_n−1) · · · (x − x₁).

師大

If the nodes are equally spaced with h = xn− x₀

n , xi = xn− (n − i)h, x = xn+ sh, then

P_n(x) = f [x_n] +

n

X

k=1

f [x_n, x_n−1, . . . , x_n−k]sh(s + 1)h · · · (s + k − 1)h

= f (x_n) +

n

X

k=1

f [x_n, x_n−1, . . . , x_n−k](−1)^kk!−s k

 h^k

Binomial formula:

−s k



= −s(−s − 1) · · · (−s − k + 1)

k! = (−1)^ks(s + 1) · · · (s + k − 1)

k! .

師大

Backward-difference:

∇^kf (x_i) = ∇^k−1f (x_i) − ∇^k−1f (x_i−1) = ∇

∇^k−1f (x_i) , then

f [xn, xn−1] = 1

h∇f (x_n), f [xn, xn−1, xn−2] = 1

2h²∇²f (xn), and, in general,

f [xn, xn−1, . . . , x_n−k] = 1

k! h^k∇^kf (xn).

Newton backward-difference formula:

Pn(x) = f (x0) +

n

X

k=1

(−1)^k−s k



∇^kf (xk).

師大

Hermite Interpolation

Given n + 1 data points x0 < x1< · · · < xn, and

y₀⁽⁰⁾= f (x₀) y₁⁽⁰⁾ = f (x₁) · · · yn⁽⁰⁾= f (x_n) y⁽¹⁾₀ = f⁰(x0) y₁⁽¹⁾= f⁰(x1) · · · y⁽¹⁾n = f⁰(xn)

... ... ...

y₀^(m0) = f^(m0)(x₀) y^(m₁ ¹⁾= f^(m1)(x₁) · · · yn^(mn)= f^(mn)(x_n)

↓ ↓ ↓

m₀+ 1 conditions m₁+ 1 conditions · · · m_n+ 1 conditions for some function f ∈ C^m[a, b] where m = max{m₀, m₁, . . . , m_n}.

師大

Determine a polynomial P of degree at most N , where N =

n

X

i=0

m_i

!

+ n, (5)

to satisfy the following interpolation conditions:

P^(k)(xi) = y^(k)_i , k = 0, 1, . . . mi, i = 0, 1, . . . , n. (6) If n = 0, then P is the m₀th Taylor polynomial for f at x₀.

If m_i = 0 for each i, then P is the nth Lagrange polynomial interpolating f on x0, . . . , xn.

If m_i = 1 for each i, then P is called the Hermite polynomials.

師大

Theorem

If f ∈ C¹[a, b] and x₀, . . . , x_n∈ [a, b] are distinct, then the polynomial of least degree agreeing with f and f⁰ at x₀, . . . , x_n is unique and is given by

H2n+1(x) =

n

X

j=0

f (xj)Hn,j(x) +

n

X

j=0

f⁰(xj) bHn,j(x),

where

H_n,j(x) =1 − 2(x − x_j)L⁰_n,j(x_j) L²_n,j(x), Hb_n,j(x) = (x − x_j)L²_n,j(x), and

Ln,j(x) = (x − x₀) · · · (x − x_j−1)(x − x_j+1) · · · (x − x_n) (xj− x₀) · · · (xj− x_j−1)(xj − x_j+1) · · · (xj− x_n). Moreover, if f ∈ C²ⁿ⁺²[a, b], then ∃ ξ(x) ∈ [a, b] s.t.

f (x) = H2n+1(x) + (x − x0)²· · · (x − x_n)²

(2n + 2)! f⁽²ⁿ⁺²⁾(ξ(x)).

師大

Proof:

Since

Ln,j(xi) =

 0, if i 6= j, 1, if i = j, it implies that

H_n,j(x_i) =

 0, if i 6= j,

1, if i = j, Hb_n,j(x_i) = 0 ∀ i, j.

As a consequence, H_2n+1(x_i) =

n

X

j=0,j6=1

f (x_j) · 0 + f (x_i) · 1 +

n

X

j=0

f⁰(x_j) · 0 = f (x_i), ∀ i.

By definitions of H_n,j(x) and bH_n,j(x),

H_n,j⁰ (x) = −2L⁰_n,j(x_j)L²_n,j(x) + 21 − 2(x − x_j)L⁰_n,j(x_j) L_n,j(x)L⁰_n,j(x), Hb_n,j⁰ (x) = L²_n,j(x) + 2(x − xj)Ln,j(x)L⁰_n,j(x).

師大

It implies that if i 6= j,

H_n,j⁰ (xi) = −2L⁰_n,j(xj)L²_n,j(xi) + 21 − 2(x_i− x_j)L⁰_n,j(xj) L_n,j(xi)L⁰_n,j(xi)

= −2L⁰_n,j(xj) · 0 + 21 − 2(x_i− x_j)L⁰_n,j(xj) · 0 · L⁰_n,j(xi) = 0, Hb_n,j⁰ (x_i) = L²_n,j(x_i) + 2(x_i− x_j)L_n,j(x_i)L⁰_n,j(x_i)

= 0 + 2(xi− x_j) · 0 · L⁰_n,j(xi) = 0 and

H_n,i⁰ (xi) = −2L⁰_n,i(xi)L²_n,i(xi) + 21 − 2(x_i− x_i)L⁰_n,i(xi) L_n,i(xi)L⁰_n,i(xi)

= −2L⁰_n,i(xi) + 2L⁰_n,i(xi) = 0,

Hb_n,i⁰ (xi) = L²_n,i(xi) + 2(xi− x_i)Ln,i(xi)L⁰_n,i(xi) = 1.

Therefore, ∀ i, H_2n+1⁰ (x_i) =

n

X

j=0

f (x_j)H_n,j⁰ (x_i) +

n

X

j=0,j6=i

f⁰(x_j) bH_n,j⁰ (x_i) + f⁰(x_i) bH_n,i⁰ (x_i) = f⁰(x_i).

That is H2n+1 agrees with f and H_2n+1⁰ with f⁰ at x0, . . . , xn.

師大

Proof of uniqueness:

Since deg(P ) ≤ 2n + 1, write

P (x) = a0+ a1x + · · · + a2n+1x²ⁿ⁺¹.

2n + 2 coefficients, a₀, a₁, . . . , a_2n+1, to be determined and 2n + 2 conditions given

P (xi) = f (xi), P⁰(xi) = f⁰(xi), for i = 0, . . . , n.

⇒ 2n + 2 linear equations in 2n + 2 unknowns to solve

⇒ show that the coefficient matrix A of this system is nonsingular.

To prove A is nonsingular, it suffices to prove that Au = 0 has only the trivial solution u = 0.

Au = 0 iff

P (x_i) = 0, P⁰(x_i) = 0, for i = 0, . . . , n.

⇒ P is a multiple of the polynomial given by q(x) =

n

Y

i=0

(x − xi)².

師大

However, deg(q) = 2n + 2 whereas P has degree at most N . Therefore, P (x) = 0, i.e. u = 0.

That is, A is nonsingular, and the Hermite interpolation problem has a unique solution.

師大

Divided Difference Method for Hermite Interpolation

Denote the N + 1 conditions pairs

(x₀, f (x₀)), (x₀, f⁰(x₀)), (x₁, f (x₁)), (x₁, f⁰(x₁)), . . . , (x_n, f (x_n)), (x_n, f⁰(x_n)) by consecutive coordinate pairs

(z0, f (x0)), (z1, f⁰(x0)), (z2, f (x1)), (z3, f⁰(x1)), . . . , (zN −1, f (xn)), (zN, f⁰(xn)).

Note that z₀ ≤ z₁ ≤ · · · ≤ z_N. The unique Hermite interpolating polynomial in Newton’s form can then be written as

H_2n+1(x) = f [z₀] +

N

X

k=1

f [z₀, z₁, . . . , z_k](x − z₀)(x − z₁) · · · (x − z_k−1).

If zi 6= z_i+k, then

f [z_i, z_i+1, . . . , z_i+k] = f [z_i+1, z_i+2, . . . , z_i+k] − f [z_i, z_i+1, . . . , z_i+k−1]

zi+k− z_i .

師大

However the divided-difference formula has to be modified because there may be repetitions among z_i. Since z_2i= z_2i+1= x_i for each i, let

ξ2i< ξ2i+1

be distinct coordinates and ξj → z_j, for 2i ≤ j ≤ 2i + 1. Then, by applying Theorem 3.6,

f [z2i, z2i+1] = lim

ξj→zj

2i≤j≤2i+1

f [ξ2i, ξ2i+1] = f⁰(z2i) = f⁰(xi)

師大

z f (z)

z0= x0 f [z0] = f (x0)

f [z0, z1] = f⁰(x0)

z₁= x₀ f [z₁] = f (x₀) f [z₀, z₁, z₂] =^{f [z1,z}_z^{2]−f [z0,z1]}

2−z0

f [z₁, z₂] = ^{f [z}_z^{2]−f [z1]}

2−z1

z2= x1 f [z2] = f (x1) f [z1, z2, z3] =^{f [z2,z}_z^{3]−f [z1,z2]}

3−z1

f [z2, z3] = f⁰(x1)

z3= x1 f [z3] = f (x1) f [z2, z3, z4] =^{f [z3,z}_z^{4]−f [z2,z3]}

4−z2

f [z₃, z₄] = ^{f [z}_z^{4]−f [z3]}

4−z3

z4= x2 f [z4] = f (x2) f [z3, z4, z5] =^{f [z4,z}_z^{5]−f [z3,z4]}

5−z3

f [z₄, z₅] = f⁰(x₂) z₅= x₂ f [z₅] = f (x₂)

師大

Cubic spline interpolation

The previous sections concern the approximation of an arbitrary function on a closed interval by a polynomial. However, the oscillatory nature of high-degree polynomials restricts their use.

Piecewise polynomial interpolation: divide the interval into a collection of subintervals and construct different approximation on each subinterval.

The simplest piecewise polynomial approximation is piecewise linear interpolation.

A disadvantage of linear function approximation is that the

interpolating function is not smooth at each of the endpoints of the subintervals. It is often required that the approximating function is continuously differentiable.

An alternative procedure is to use a piecewise polynomial of Hermite type. However, to use Hermite piecewise polynomials for general interpolation, we need to know the derivatives of the function being approximated, which is frequently not available.

師大

The most common piecewise-polynomial use cubic polynomials between each successive pair of nodes and is called cubic spline interpolation.

Definition

Given a function f defined on [a, b], and a set of nodes

a = x0 < x1 < · · · < xn= b, a cubic spline interpolation S for f is a function that satisfies the following conditions:

(a) S is a cubic polynomial, denoted S_j(x), on [x_j, x_j+1] for each j = 0, 1, . . . , n − 1;

(b) S_j(x_j) = f (x_j) and S_j(x_j+1) = f (x_j+1) ∀ j = 0, 1, . . . , n − 1;

(c) Sj+1(xj+1) = Sj(xj+1)∀ j = 0, 1, . . . , n − 2;

(d) S_j+1⁰ (xj+1) = S_j⁰(xj+1)∀ j = 0, 1, . . . , n − 2;

(e) S_j+1⁰⁰ (xj+1) = S_j⁰⁰(xj+1)∀ j = 0, 1, . . . , n − 2;

(f) One of the following sets of boundary conditions is satisfied:

(i) S⁰⁰(x0) = S⁰⁰(xn) = 0 (free or natural boundary);

(ii) S⁰(x0) = f⁰(x0) and S⁰(xn) = f⁰(xn) (clamped boundary).

師大

Construction of cubic spline interpolation

Remark

In general, clamped boundary conditions lead to more accurate

approximations since they include more information about the function.

Write

S_j(x) = a_j+ b_j(x − x_j) + c_j(x − x_j)²+ d_j(x − x_j)³, ∀ j = 0, 1, . . . , n − 1.

That is

a_j = S_j(x_j) = f (x_j).

Since Sj+1(xj+1) = Sj(xj+1), it implies that

aj+1= aj+ bj(xj+1− x_j) + cj(xj+1− x_j)²+ dj(xj+1− x_j)³. Let

h_j = x_j+1− x_j, ∀ j = 0, 1, . . . , n − 1.

Then

aj+1= aj+ bjhj + cjh²_j + djh³_j. (7)

師大

Define b_n= S⁰(x_n) and observe that

S_j⁰(x) = b_j+ 2c_j(x − x_j) + 3d_j(x − x_j)² implies that

bj = S_j⁰(xj), ∀ j = 0, 1, . . . , n − 1.

Applying S_j+1⁰ (x_j+1) = S_j⁰(x_j+1) gives

bj+1= bj+ 2cjhj+ 3djh²_j, ∀ j = 0, 1, . . . , n − 1. (8) Defining c_n= S⁰⁰(x_n)/2 and applying S_j+1⁰⁰ (x_j+1) = S_j⁰⁰(x_j+1), we get

c_j+1 = c_j+ 3d_jh_j, ∀ j = 0, 1, . . . , n − 1.

Hence

dj = c_j+1− c_j 3hj

. (9)

師大

Substituting dj in (9) into (7) and (8), it obtains a_j+1= a_j+ b_jh_j +h²_j

3 (2c_j+ c_j+1) (10)

and

bj+1 = bj+ hj(cj+ cj+1). (11) Eq. (10) implies that

b_j = 1

h_j(a_j+1− a_j) − hj

3 (2c_j+ c_j+1), (12) and then

bj−1 = 1

h_j−1(aj− a_j−1) −hj−1

3 (2cj−1+ cj). (13)

師大

Substituting (12) and (13) into (11), we get hj−1cj−1+ 2(hj−1+ hj)cj+ hjcj+1 = 3

h_j(aj+1− a_j) − 3

h_j−1(aj− a_j−1).(14) for each j = 1, 2, . . . , n − 1.

Theorem

If f is defined at a = x0 < x1 < · · · < xn= b, then f has a unique natural spline interpolant S on the nodes x0, x1, . . . , xn that satisfies

S⁰⁰(a) = S⁰⁰(b) = 0.

Proof: The boundary conditions imply that cn = S⁰⁰(xn)/2 = 0,

0 = S⁰⁰(x₀) = 2c₀+ 6d₀(x₀− x₀) ⇒ c₀= 0.

師大

Substituting c0= cn= 0 into (14), it produces a linear system Ax = b, where A is an (n + 1) × (n + 1) matrix

A =



















1 0

h₀ 2(h₀+ h₁) h₁

h1 2(h1+ h2) h2

. .. . .. . ..

hn−2 2(hn−2+ hn−1) hn−1

0 1

















 ,

b =















0

3

h1(a₂− a₁) −_h³

0(a₁− a₀) ...

3

hn−1(a_n− a_n−1) −_h³

n−2(a_n−1− a_n−2) 0















, x =









 c₀ c₁ ... c_n









 .

The matrix A is strictly diagonally dominant, so A is nonsingular and the linear system has a unique solution for c0, c1, . . . , cn.

師大

Solving tridiagonal linear system

Let

A = LU ≡











`₁₁

`21 `22

. .. . ..

`<sub>n,n−1</sub> `_nn



























1 u₁₂ 1 u₂₃

. .. ...

. .. u_n−1,n 1















 .

Comparing the coefficients of matrix in both sides, we have a11 = `11,

ai,i−1 = `i,i−1, for i = 2, 3, . . . , n,

a_ii = `<sub>i,i−1</sub>u<sub>i−1,i</sub>+ `_ii, for i = 2, 3, . . . , n, a_i,i+1 = `_iiu_i,i+1, for i = 1, 2, . . . , n − 1.

師大

It implies that

`<sub>11</sub>= a<sub>11</sub>, u<sub>12</sub>= a<sub>12</sub>/`₁₁, and, for i = 2, . . . , n − 1,

`i,i−1 = ai,i−1,

`ii = aii− `_i,i−1ui−1,i, ui,i+1 = ai,i+1/`ii, and

`_n,n−1 = a_n,n−1,

`nn = ann− `_n,n−1un−1,n.

師大

Forward substitution

Solve Ly = b

⇒











`11y1 = b1,

`21y1+ `22y2 = b2, ...

`n,n−1yn−1+ `nnyn = bn.

⇒











y1 = b1/`11,

y₂ = (b₂− `<sub>21</sub>y<sub>1</sub>)/`₂₂, ...

yn = (bn− `<sub>n,n−1</sub>yn−1)/`nn.

⇒

 y₁ = b₁/`₁₁,

yi = (bi− `<sub>i,i−1</sub>yi−1)/`ii, for i = 2, . . . , n.

師大

Back substitution

Solve U x = y,

⇒











x1+ u12x2 = y1, ...

x_n−1+ u_n−1,nx_n = y_n−1, xn = yn,

⇒











xn = yn,

x_n−1 = y_n−1− u_n−1,nx_n, ...

x1 = y1− u₁₂x2,

⇒

 x_n = y_n,

xi = yi− u_i,i+1xi+1, for i = n − 1, . . . , 1.

師大

Theorem

If f is defined at a = x0 < x1< · · · < xn= b and differentiable at a and b, then f has a unique clamped spline interpolant S on the nodes

x₀, x₁, . . . , x_n that satisfies S⁰(a) = f⁰(a) and S⁰(b) = f⁰(b).

Proof: Since f⁰(a) = S⁰(a) = S⁰(x₀) = b₀, Eq. (12) with j = 0 implies f⁰(a) = 1

h0

(a1− a₀) −h0

3 (2c0+ c1).

Consequently,

2h₀c₀+ h₀c₁= 3 h0

(a₁− a₀) − 3f⁰(a).

Similarly,

f⁰(b) = b_n= b_n−1+ h_n−1(c_n−1+ c_n),

師大

so Eq. (12) with j = n − 1 implies that f⁰(b) = an− a_n−1

hn−1

−hn−1

3 (2cn−1+ cn) + hn−1(cn−1+ cn)

= an− a_n−1

h_n−1 +hn−1

3 (c_n−1+ 2c_n), and

h_n−1c_n−1+ 2h_n−1c_n= 3f⁰(b) − 3 hn−1

(a_n− a_n−1).

Eq. (14) together with the equations 2h0c0+ h0c1 = 3

h0

(a1− a₀) − 3f⁰(a) and

hn−1cn−1+ 2hn−1cn= 3f⁰(b) − 3

h_n−1(an− a_n−1)