De¬ne e = f ’ Π1 f . Since e(xj ) = 0 for all j = 0, . . . , K, Rolle™s theorem

h

infers the existence of ξj ∈ (xj , xj+1 ), for j = 0, . . . , K ’ 1 such that e (ξj ) = 0.

Since Π1 f is a linear function on each Ij , for x ∈ Ij we obtain

h

x x

e (x) = e (s)ds = f (s)ds,

ξj ξj

whence

xj+1

|e (x)| ¤ |f (s)|ds, for x ∈ [xj , xj+1 ]. (8.28)

xj

We recall the Cauchy-Schwarz inequality

1/2 1/2

β β β

u(x)v(x)dx ¤ 2 2

u (x)dx v (x)dx (8.29)

± ± ±

which holds if u, v ∈ L2 (±, β). If we apply this inequality to (8.28) we obtain

« 1/2 « 1/2

xj+1 xj+1

¬ · ¬ ·

|e (x)| ¤ |f (s)|2 ds

12 dx

xj xj

« 1/2 (8.30)

xj+1

¬ ·

¤ h1/2 |f (s)|2 ds .

xj

To ¬nd a bound for |e(x)|, we notice that

x

e(x) = e (s)ds,

xj

so that, applying (8.30), we get

1/2

xj+1 xj+1

|e(x)| ¤ |e (s)|ds ¤ h |f (s)| ds

3/2 2

. (8.31)

xj xj

Then

xj+1 xj+1 xj+1 xj+1

|e (x)|2 dx ¤ h2 |f (s)|2 ds |e(x)|2 dx ¤ h4 |f (s)|2 ds,

and

xj xj xj xj

from which, summing over the index j from 0 to K ’ 1 and taking the square

root of both sides, we obtain

1/2 1/2

b b

|e (x)| dx ¤h |f (x)| dx

2 2

,

a a

and

1/2 1/2

b b

|e(x)| dx ¤h |f (x)| dx

2 2 2

,

a a

3

which is the desired estimate (8.27), with C1 = C2 = 1.

8.4 Hermite-Birko¬ Interpolation 341

8.4 Hermite-Birko¬ Interpolation

Lagrange polynomial interpolation can be generalized to the case in which

also the values of the derivatives of a function f are available at some (or

all) of the nodes xi .

Let us then suppose that (xi , f (k) (xi )) are given data, with i = 0, . . . , n,

n

k = 0, . . . , mi and mi ∈ N. Letting N = i=0 (mi + 1), it can be proved

(see [Dav63]) that, if the nodes {xi } are distinct, there exists a unique

polynomial HN ’1 ∈ PN ’1 , called the Hermite interpolation polynomial,

such that

(k) (k)

HN ’1 (xi ) = yi , i = 0, . . . , n k = 0, . . . , mi ,

of the form

n mi

(k)

HN ’1 (x) = yi Lik (x) (8.32)

i=0 k=0

(k)

where yi = f (k) (xi ), i = 0, . . . , n, k = 0, . . . , mi .

The functions Lik ∈ PN ’1 are called the Hermite characteristic polynomials

and are de¬ned through the relations

1 if i = j and k = p,

dp

(Lik )(xj ) =

dxp 0 otherwise.

De¬ning the polynomials

n mk +1

(x ’ xi )j x ’ xk

lij (x) = , i = 0, . . . , n, j = 0, . . . , mi ,

xi ’ xk

j! k=0

k=i

and letting Limi (x) = limi (x) for i = 0, . . . , n, we have the following recur-

sive formula for the polynomials Lij

mi

(k)

Lij (x) = lij (x) ’ j = mi ’ 1, mi ’ 2, . . . , 0.

lij (xi )Lik (x)

k=j+1

As for the interpolation error, the following estimate holds

f (N ) (ξ)

f (x) ’ HN ’1 (x) = „¦N (x) ∀x ∈ R

N!

where ξ ∈ I(x; x0 , . . . , xn ) and „¦N is the polynomial of degree N de¬ned

by

„¦N (x) = (x ’ x0 )m0 +1 (x ’ x1 )m1 +1 . . . (x ’ xn )mn +1 . (8.33)

342 8. Polynomial Interpolation

Example 8.6 (osculatory interpolation) Let us set mi = 1 for i = 0, . . . , n.

In this case N = 2n + 2 and the interpolating Hermite polynomial is called the

osculating polynomial, and it is given by

n

(1)

HN ’1 (x) = yi Ai (x) + yi Bi (x)

i=0

where Ai (x) = (1 ’ 2(x ’ xi )li (xi ))li (x)2 and Bi (x) = (x ’ xi )li (x)2 , for i =

0, . . . , n, with

n

1

li (xi ) = , i = 0, . . . , n.

xi ’ xk

k=0,k=i

As a comparison, we use Programs 65 and 67 to compute the Lagrange and

Hermite interpolating polynomials of the function f (x) = sin(4πx) on the interval

[0, 1] taking four equally spaced nodes (n = 3). Figure 8.4 shows the superposed

graphs of the function f (dashed line) and of the two polynomials Πn f (dotted

•

line) and HN ’1 (solid line).

1.5

1

0.5

0

’0.5

’1

’1.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

FIGURE 8.4. Lagrange and Hermite interpolation for the function

f (x) = sin(4πx) on the interval [0, 1]

Program 67 computes the values of the osculating polynomial at the ab-

scissae contained in the vector z. The input vectors x, y and dy contain the

interpolation nodes and the corresponding function evaluations of f and

f , respectively.

Program 67 - hermpol : Osculating polynomial

function [herm] = hermite(x,y,dy,z)

n = max(size(x)); m = max(size(z)); herm = [];

for j = 1:m

xx = z(j); hxv = 0;

for i = 1:n,

den = 1; num = 1; xn = x(i); derLi = 0;

for k = 1:n,

8.5 Extension to the Two-Dimensional Case 343

if k ˜= i, num = num*(xx-x(k)); arg = xn-x(k);

den = den*arg; derLi = derLi+1/arg;

end

end

Lix2 = (num/den)ˆ2; p = (1-2*(xx-xn)*derLi)*Lix2;

q = (xx-xn)*Lix2; hxv = hxv+(y(i)*p+dy(i)*q);

end

herm = [herm, hxv];

end

8.5 Extension to the Two-Dimensional Case

In this section we brie¬‚y address the extension of the previous concepts to

the two-dimensional case, referring to [SL89], [CHQZ88], [QV94] for more

details. We denote by „¦ a bounded domain in R2 and by x = (x, y) the

coordinate vector of a point in „¦.

8.5.1 Polynomial Interpolation

A particularly simple situation occurs when „¦ = [a, b] — [c, d], i.e., the

interpolation domain „¦ is the tensor product of two intervals. In such a

case, introducing the nodes a = x0 < x1 < . . . < xn = b and c = y0 <

y1 < . . . < ym = d, the interpolating polynomial Πn,m f can be written as

n m

Πn,m f (x, y) = i=0 j=0 ±ij li (x)lj (y), where li ∈ Pn , i = 0, . . . , n, and

lj ∈ Pm , j = 0, . . . , m, are the characteristic one-dimensional Lagrange

polynomials with respect to the x and y variables respectively, and where

±ij = f (xi , yj ).

The drawbacks of one-dimensional Lagrange interpolation are inherited

by the two-dimensional case, as con¬rmed by the example in Figure 8.5.

Remark 8.1 (The general case) If „¦ is not a rectangular domain or if

the interpolation nodes are not uniformly distributed over a Cartesian grid,

the interpolation problem is di¬cult to solve, and, generally speaking, it is

preferable to resort to a least-squares solution (see Section 10.7). We also

point out that in d dimensions (with d ≥ 2) the problem of ¬nding an

interpolating polynomial of degree n with respect to each space variable on

n + 1 distinct nodes might be ill-posed.

Consider, for example, a polynomial of degree 1 with respect to x and y

of the form p(x, y) = a3 xy +a2 x+a1 y +a0 to interpolate a function f at the

nodes (’1, 0), (0, ’1), (1, 0) and (0, 1). Although the nodes are distinct, the

problem (which is nonlinear) does not in general admit a unique solution;

actually, imposing the interpolation constraints, we end up with a system

that is satis¬ed by any value of the coe¬cient a3 .

344 8. Polynomial Interpolation

8

0.5

0.4 6

0.3

4

0.2

2

0.1

0

0

’2

’0.1

5

5

5

5

0

0

0

0

’5

’5 ’5

’5

FIGURE 8.5. Runge™s counterexample extended to the two-dimensional case:

interpolating polynomial on a 6 — 6 nodes grid (left) and on a 11 — 11 nodes grid

(right). Notice the change in the vertical scale between the two plots

8.5.2 Piecewise Polynomial Interpolation

In the multidimensional case, the higher ¬‚exibility of piecewise interpola-

tion allows for easy handling of domains of complex shape. Let us suppose

that „¦ is a polygon in R2 . Then, „¦ can be partitioned into K nonover-

lapping triangles (or elements) T , which de¬ne the so called triangulation

of the domain which will be denoted by Th . Clearly, „¦ = T . Suppose

T ∈Th

that the maximum length of the edges of the triangles is less than a positive

number h. As shown in Figure 8.6 (left), not any arbitrary triangulation is

allowed. Precisely, the admissible ones are those for which any pair of non

disjoint triangles may have a vertex or an edge in common.

T1

T1

T2

aT

T2

y y 3

FT

1

aT T

1

T1

aT

T 2

T1

T2 x x

T2

0 1

FIGURE 8.6. The left side picture shows admissible (above) and non admissible

(below) triangulations while the right side picture shows the a¬ne map from the

ˆ

reference triangle T to the generic element T ∈ Th

Any element T ∈ Th , of area equal to |T |, is the image through the a¬ne

ˆ

map x = FT (ˆ ) = BT x + bT of the reference triangle T , of vertices (0,0),

x

8.5 Extension to the Two-Dimensional Case 345

ˆ

(1,0) and (0,1) in the x = (ˆ, y ) plane (see Figure 8.6, right), where the

xˆ

invertible matrix BT and the right-hand side bT are given respectively by

x2 ’ x1 x3 ’ x1

, bT = (x1 , y1 )T ,

BT = (8.34)

y2 ’ y1 y3 ’ y1

(l)

while the coordinates of the vertices of T are denoted by aT = (xl , yl )T

for l = 1, 2, 3.

li (x,y)

li (x,y)

1

1

zi

zi

li (x) li (x)

1

1

zi

zi

FIGURE 8.7. Characteristic piecewise Lagrange polynomial, in one and two space

dimensions. Left, k = 0; right, k = 1

The a¬ne map (8.34) is of remarkable importance in practical computa-

tions, since, once a basis has been generated for representing the piecewise

ˆ

polynomial interpolant on T , it is possible, applying the change of coor-

dinates x = FT (ˆ ), to reconstruct the polynomial on each element T of

x

Th . We are thus interested in devising local basis functions, which can be

fully described over each triangle without needing any information from

adjacent triangles.

For this purpose, let us introduce on Th the set Z of the piecewise interpo-

lation nodes zi = (xi , yi )T , for i = 1, . . . , N , and denote by Pk („¦), k ≥ 0,

the space of algebraic polynomials of degree ¤ k in the space variables x, y

±

k

Pk („¦) = p(x, y) = aij x y , x, y ∈ „¦ .

ij

(8.35)

i,j=0

i+j¤k

Finally, for k ≥ 0, let Pc („¦) be the space of piecewise polynomials of degree

k

¤ k, such that, for any p ∈ Pc („¦), p|T ∈ Pk (T ) for any T ∈ Th . An ele-

k

mentary basis for Pc („¦) consists of the Lagrange characteristic polynomials

k

li = li (x, y), such that li ∈ Pc („¦) and

k

li (zj ) = δij , i, j = 1, . . . , N, (8.36)

346 8. Polynomial Interpolation

where δij is the Kronecker symbol. We show in Figure 8.7 the functions li for

k = 0, 1, together with their corresponding one-dimensional counterparts.