0.8 1.6755 -0.64 -4.13 -2.08

1.2 0.5575 -2.79 -2.69 1.43 2.93

1.6 0.0038 -1.38 1.76 3.71 1.62 -0.81

2.0 0.7206 1.79 3.97 1.83 -1.17 -1.55 -0.36

TABLE 8.1. Divided di¬erences for the function f (x) = 1 + sin(3x) in the case

in which the evaluation of f at x = 0.2 is also available. The newly computed

values are denoted by italics

8.2 Newton Form of the Interpolating Polynomial 337

Notice that f [x0 , . . . , xn ] = 0 for any f ∈ Pn’1 . This property, how-

ever, is not always veri¬ed numerically, since the computation of divided

di¬erences might be highly a¬ected by rounding errors.

Example 8.4 Consider again the divided di¬erences for the function f (x) =

1 + sin(3x) on the interval (0, 0.0002). The function behaves like 1 + 3x in a

su¬ciently small neighbourhood of 0, so that we expect to ¬nd smaller numbers as

the order of divided di¬erences increases. However, the results obtained running

Program 66, and shown in Table 8.2 in exponential notation up to the ¬rst 4

signi¬cant ¬gures (although 16 digits have been employed in the calculations),

exhibit a substantially di¬erent pattern. The small rounding errors introduced in

the computation of divided di¬erences of low order have dramatically propagated

•

on the higher order divided di¬erences.

D2 fi D3 fi D4 fi D 5 fi

xi f (xi ) f [xi , xi’1 ]

0 1.0000

4.0e-5 1.0001 3.000

8.0e-5 1.0002 3.000 -5.39e-4

1.2e-4 1.0004 3.000 -1.08e-3 -4.50

1.6e-4 1.0005 3.000 -1.62e-3 -4.49 1.80e+1

’1.2e + 5

2.0e-4 1.0006 3.000 -2.15e-3 -4.49 -7.23

TABLE 8.2. Divided di¬erences for the function f (x) = 1+sin(3x) on the interval

(0,0.0002). Notice the completely wrong value in the last column (it should be

approximately equal to 0), due to the propagation of rounding errors throughout

the algorithm

8.2.2 The Interpolation Error Using Divided Di¬erences

Consider the nodes x0 , . . . , xn and let Πn f be the interpolating polynomial

of f on such nodes. Now let x be a node distinct from the previous ones;

letting xn+1 = x, we denote by Πn+1 f the interpolating polynomial of f

at the nodes xk , k = 0, . . . , n + 1. Using the Newton divided di¬erences

formula, we get

Πn+1 f (t) = Πn f (t) + (t ’ x0 ) . . . (t ’ xn )f [x0 , . . . , xn , t].

Since Πn+1 f (x) = f (x), we obtain the following formula for the interpola-

tion error at t = x

En (x) = f (x) ’ Πn f (x) = Πn+1 f (x) ’ Πn f (x)

= (x ’ x0 ) . . . (x ’ xn )f [x0 , . . . , xn , x] (8.20)

= ωn+1 (x)f [x0 , . . . , xn , x].

338 8. Polynomial Interpolation

Assuming f ∈ C (n+1) (Ix ) and comparing (8.20) with (8.7), yields

f (n+1) (ξ)

f [x0 , . . . , xn , x] = (8.21)

(n + 1)!

for a suitable ξ ∈ Ix . Since (8.21) resembles the remainder of the Tay-

lor series expansion of f , the Newton formula (8.17) for the interpolating

polynomial is often regarded as being a truncated expansion around x0

provided that |xn ’ x0 | is not too big.

8.3 Piecewise Lagrange Interpolation

In Section 8.1.1 we have outlined the fact that, for equally spaced inter-

polating nodes, uniform convergence of Πn f to f is not guaranteed as

n ’ ∞. On the other hand, using equally spaced nodes is clearly computa-

tionally convenient and, moreover, Lagrange interpolation of low degree is

su¬ciently accurate, provided su¬ciently small interpolation intervals are

considered.

Therefore, it is natural to introduce a partition Th of [a, b] into K subin-

tervals Ij = [xj , xj+1 ] of length hj , with h = max0¤j¤K’1 hj , such that

[a, b] = ∪K’1 Ij and then to employ Lagrange interpolation on each Ij

j=0

(i)

using n + 1 equally spaced nodes xj , 0 ¤ i ¤ n with a small n.

For k ≥ 1, we introduce on Th the piecewise polynomial space

Xh = v ∈ C 0 ([a, b]) : v|Ij ∈ Pk (Ij ) ∀Ij ∈ Th

k

(8.22)

which is the space of the continuous functions over [a, b] whose restric-

tions on each Ij are polynomials of degree ¤ k. Then, for any continuous

function f in [a, b], the piecewise interpolation polynomial Πk f coincides

h

on each Ij with the interpolating polynomial of f|Ij at the n + 1 nodes

(i)

xj , 0 ¤ i ¤ n . As a consequence, if f ∈ C k+1 ([a, b]), using (8.7) within

each interval we obtain the following error estimate

f ’ Πk f ¤ Chk+1 f (k+1) ∞. (8.23)

∞

h

Note that a small interpolation error can be obtained even for low k pro-

vided that h is su¬ciently “small”.

Example 8.5 Let us go back to the function of Runge™s counterexample. Now,

piecewise polynomials of degree k = 1 and k = 2 are employed. We check ex-

perimentally for the behavior of the error as h decreases. In Table 8.3 we show

the absolute errors measured in the maximum norm over the interval [’5, 5] and

the corresponding estimates of the convergence order p with respect to h. Except

when using an excessively small number of subintervals, the results con¬rm the

•

theoretical estimate (8.23), that is p = k + 1.

8.3 Piecewise Lagrange Interpolation 339

f ’ Πh ∞ f ’ Πh ∞

h p p

1 2

5 0.4153 0.0835

2.5 0.1787 1.216 0.0971 -0.217

1.25 0.0631 1.501 0.0477 1.024

0.625 0.0535 0.237 0.0082 2.537

0.3125 0.0206 1.374 0.0010 3.038

0.15625 0.0058 1.819 1.3828e-04 2.856

0.078125 0.0015 1.954 1.7715e-05 2.964

TABLE 8.3. Interpolation error for Lagrange piecewise interpolation of degree

k = 1 and k = 2, in the case of Runge™s function (8.12); p denotes the trend of

the exponent of h. Notice that, as h ’ 0, p ’ k + 1, as predicted by (8.23)

Besides estimate (8.23), convergence results in integral norms exist (see

[QV94], [EEHJ96]). For this purpose, we introduce the following space

±

b

L (a, b) = f : (a, b) ’ R, |f (x)| dx < +∞ ,

2 2

(8.24)

a

with

« 1/2

b

= |f (x)|2 dx

f . (8.25)

L2 (a,b)

a

Formula (8.25) de¬nes a norm for L2 (a, b). (We recall that norms and semi-

norms of functions can be de¬ned in a manner similar to what was done in

De¬nition 1.17 in the case of vectors). We warn the reader that the integral

of the function |f |2 in (8.24) has to be intended in the Lebesgue sense (see,

e.g., [Rud83]). In particular, f needs not be continuous everywhere.

Theorem 8.3 Let 0 ¤ m ¤ k + 1, with k ≥ 1 and assume that f (m) ∈

L2 (a, b) for 0 ¤ m ¤ k + 1; then there exists a positive constant C, inde-

pendent of h, such that

(f ’ Πk f )(m) ¤ Chk+1’m f (k+1) L2 (a,b) . (8.26)

L2 (a,b)

h

In particular, for k = 1, and m = 0 or m = 1, we obtain

f ’ Π1 f ¤ C1 h2 f L2 (a,b) ,

L2 (a,b)

h

(8.27)

(f ’ Π1 f ) ¤ C2 h f ,

L2 (a,b) L2 (a,b)

h

for two suitable positive constants C1 and C2 .

Proof. We only prove (8.27) and refer to [QV94], Chapter 3 for the proof of

(8.26) in the general case.

340 8. Polynomial Interpolation

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.

In the case k = 0, the interpolation nodes are collocated at the centers of

gravity of the triangles, while in the case k = 1 the nodes coincide with

the vertices of the triangles. This choice, that we are going to maintain

henceforth, is not the only one possible. The midpoints of the edges of the

triangles could be used as well, giving rise to a discontinuous piecewise

polynomial over „¦.

For k ≥ 0, the Lagrange piecewise interpolating polynomial of f , Πk f ∈

h

Pk („¦), is de¬ned as

c

N

Πk f (x, y) = f (zi )li (x, y). (8.37)

h

i=1

Notice that Π0 f is a piecewise constant function, while Π1 f is a linear

h h

function over each triangle, continuous at the vertices, and thus globally

continuous.

For any T ∈ Th , we shall denote by Πk f the restriction of the piecewise

T

interpolating polynomial of f over the element T . By de¬nition, Πk f ∈

T

Pk (T ); noticing that dk = dimPk (T ) = (k + 1)(k + 2)/2, we can therefore

write

dk ’1

(m)

∀T ∈ Th .

Πk f (x, y) = f (˜T )lm,T (x, y),

z (8.38)

T

m=0

(m)

In (8.38), we have denoted by zT , for m = 0, . . . , dk ’ 1, the piecewise

˜

interpolation nodes on T and by lm,T (x, y) the restriction to T of the La-

grange characteristic polynomial having index i in (8.37) which corresponds

(m)

˜

in the list of the “global” nodes zi to that of the “local” node zT .

’1

Keeping on with this notation, we have lj,T (x) = ˆj —¦ FT (x), where

l

ˆj = ˆj (ˆ ) is, for j = 0, . . . , dk ’ 1, the j-th Lagrange basis function for

l lx

ˆ ˆ

Pk (T ) generated on the reference element T . We notice that if k = 0 then

d0 = 1, that is, only one local interpolation node exists (coinciding with

the center of gravity of the triangle T ), while if k = 1 then d1 = 3, that is,

three local interpolation nodes exist, coinciding with the vertices of T . In

ˆ

Figure 8.8 we draw the local interpolation nodes on T for k = 0, 1 and 2.

As for the interpolation error estimate, denoting for any T ∈ Th by hT the

maximum length of the edges of T , the following result holds (see for the

proof, [CL91], Theorem 16.1, pp. 125-126 and [QV94], Remark 3.4.2, pp.

89-90)

f ’ Πk f ¤ Chk+1 f (k+1) k ≥ 0,

∞,T , (8.39)

∞,T

T T

where for every g ∈ C 0 (T ), g ∞,T = maxx∈T |g(x)|. In (8.39), C is a

positive constant independent of hT and f .

8.5 Extension to the Two-Dimensional Case 347

ˆ

FIGURE 8.8. Local interpolation nodes on T ; left, k = 0, center k = 1, right,

k=2

Let us assume that the triangulation Th is regular, i.e., there exists a

positive constant σ such that

hT

¤ σ,

max

T ∈Th ρT

where ∀T ∈ Th , ρT is the diameter of the inscribed circle to T , Then, it

is possible to derive from (8.39) the following interpolation error estimate

over the whole domain „¦

f ’ Πk f ¤ Chk+1 f (k+1) k ≥ 0, ∀f ∈ C k+1 („¦).

∞,„¦ ,

∞,„¦

h

(8.40)

The theory of piecewise interpolation is a basic tool of the ¬nite element

method, a computational technique that is widely used in the numerical

approximation of partial di¬erential equations (see Chapter 12 for the one-

dimensional case and [QV94] for a complete presentation of the method).

Example 8.7 We compare the convergence of the piecewise polynomial interpo-

2 2

lation of degree 0, 1 and 2, on the function f (x, y) = e’(x +y ) on „¦ = (’1, 1)2 .

We show in Table 8.4 the error Ek = f ’ Πk f ∞,„¦ , for k = 0, 1, 2, and the order

h

of convergence pk as a function of the mesh size h = 2/N for N = 2, . . . , 32.

Clearly, linear convergence is observed for interpolation of degree 0 while the

order of convergence is quadratic with respect to h for interpolation of degree 1

•

and cubic for interpolation of degree 2.

E0 E1 E2

h p0 p1 p2

1 0.4384 0.2387 0.016

1.6678 · 10’3

1

0.2931 0.5809 0.1037 1.2028 3.2639

2

2.8151 · 10’4

1

0.1579 0.8924 0.0298 1.7990 2.5667

4

3.5165 · 10’5

1

0.0795 0.9900 0.0077 1.9524 3.001

8

4.555 · 10’6

1

0.0399 0.9946 0.0019 2.0189 2.9486

16

TABLE 8.4. Convergence rates and orders for piecewise interpolations of degree

0, 1 and 2

348 8. Polynomial Interpolation

8.6 Approximation by Splines

In this section we address the matter of approximating a given function us-

ing splines, which allow for a piecewise interpolation with a global smooth-

ness.

De¬nition 8.1 Let x0 , . . . , xn , be n + 1 distinct nodes of [a, b], with a =

x0 < x1 < . . . < xn = b. The function sk (x) on the interval [a,b] is a spline

of degree k relative to the nodes xj if

sk|[xj ,xj+1 ] ∈ Pk , j = 0, 1, . . . , n ’ 1 (8.41)

sk ∈ C k’1 [a, b]. (8.42)

Denoting by Sk the space of splines sk on [a, b] relative to n + 1 distinct

nodes, then dim Sk = n + k. Obviously, any polynomial of degree k on [a, b]

is a spline; however, in the practice a spline is represented by a di¬erent

polynomial on each subinterval and for this reason there could be a discon-

tinuity in its k-th derivative at the internal nodes x1 , . . . , xn’1 . The nodes

for which this actually happens are called active nodes.

It is simple to check that conditions (8.41) and (8.42) do not su¬ce to

characterize a spline of degree k. Indeed, the restriction sk,j = sk|[xj ,xj+1 ]

can be represented as

k

sij (x ’ xj )i , if x ∈ [xj , xj+1 ] (8.43)

sk,j (x) =

i=0

so that (k + 1)n coe¬cients sij must be determined. On the other hand,

from (8.42) it follows that

(m) (m)

sk,j’1 (xj ) = sk,j (xj ), j = 1, . . . , n ’ 1, m = 0, ..., k ’ 1

which amounts to setting k(n ’ 1) conditions. As a consequence, the re-

maining degrees of freedom are (k + 1)n ’ k(n ’ 1) = k + n.

Even if the spline were interpolatory, that is, such that sk (xj ) = fj for

j = 0, . . . , n, where f0 , . . . , fn are given values, there would still be k ’ 1

unsaturated degrees of freedom. For this reason further constraints are

usually imposed, which lead to:

1. periodic splines, if

(m) (m)

sk (a) = sk (b), m = 0, 1, . . . , k ’ 1; (8.44)

8.6 Approximation by Splines 349

2. natural splines, if for k = 2l ’ 1, with l ≥ 2

(l+j) (l+j)

(b) = 0, j = 0, 1, . . . , l ’ 2. (8.45)

sk (a) = sk

From (8.43) it turns out that a spline can be conveniently represented using

k + n spline basis functions, such that (8.42) is automatically satis¬ed. The

simplest choice, which consists of employing a suitably enriched monomial

basis (see Exercise 10), is not satisfactory from the numerical standpoint,

since it is ill-conditioned. In Sections 8.6.1 and 8.6.2 possible examples of

spline basis functions will be provided: cardinal splines for the speci¬c case

k = 3 and B-splines for a generic k.

8.6.1 Interpolatory Cubic Splines

Interpolatory cubic splines are particularly signi¬cant since: i. they are

the splines of minimum degree that yield C 2 approximations; ii. they are

su¬ciently smooth in the presence of small curvatures.

Let us thus consider, in [a, b], n + 1 ordered nodes a = x0 < x1 < . . . <

xn = b and the corresponding evaluations fi , i = 0, . . . , n. Our aim is to

provide an e¬cient procedure for constructing the cubic spline interpolating

those values. Since the spline is of degree 3, its second-order derivative must

be continuous. Let us introduce the following notation

fi = s3 (xi ), mi = s3 (xi ), Mi = s3 (xi ), i = 0, . . . , n.

Since s3,i’1 ∈ P3 , s3,i’1 is linear and

xi ’ x x ’ xi’1

for x ∈ [xi’1 , xi ]

s3,i’1 (x) = Mi’1 + Mi (8.46)

hi hi

where hi = xi ’ xi’1 . Integrating (8.46) twice we get

(xi ’ x)3 (x ’ xi’1 )3

+ Ci’1 (x ’ xi’1 ) + Ci’1 ,

s3,i’1 (x) = Mi’1 + Mi

6hi 6hi

and the constants Ci’1 and Ci’1 are determined by imposing the end point

values s3 (xi’1 ) = fi’1 and s3 (xi ) = fi . This yields, for i = 1, . . . , n ’ 1

fi ’ fi’1

h2 hi

= fi’1 ’ Mi’1 i , Ci’1 = ’ (Mi ’ Mi’1 ).

Ci’1

6 hi 6

Let us now enforce the continuity of the ¬rst derivatives at xi ; we get

fi ’ fi’1

hi hi

s3 (x’ ) = Mi’1 + Mi +

i

6 3 hi

fi+1 ’ fi

hi+1 hi+1

=’ Mi ’ = s3 (x+ ),

Mi+1 + i

3 6 hi+1

350 8. Polynomial Interpolation

where s3 (x± ) = lim s3 (xi ± t). This leads to the following linear system

i t’0

(called M-continuity system)

i = 1, . . . , n ’ 1

µi Mi’1 + 2Mi + »i Mi+1 = di (8.47)

where we have set

hi hi+1

µi = , »i = ,

hi + hi+1 hi + hi+1

fi+1 ’ fi fi ’ fi’1

6

’ i = 1, . . . , n ’ 1.

di = ,

hi + hi+1 hi+1 hi

System (8.47) has n + 1 unknowns and n ’ 1 equations; thus, 2(= k ’ 1)

conditions are still lacking. In general, these conditions can be of the form

2M0 + »0 M1 = d0 , µn Mn’1 + 2Mn = dn ,

with 0 ¤ »0 , µn ¤ 1 and d0 , dn given values. For instance, in order to

obtain the natural splines (satisfying s3 (a) = s3 (b) = 0), we must set the

above coe¬cients equal to zero. A popular choice sets »0 = µn = 1 and

d0 = d1 , dn = dn’1 , which corresponds to prolongating the spline outside

the end points of the interval [a, b] and treating a and b as internal points.

This strategy produces a spline with a “smooth” behavior. In general, the

resulting linear system is tridiagonal of the form

®

® ®

2 »0 0 ... 0

M0 d0

.

µ1 2 M1 d1

.

»1 .

. .

.. .. ..

0 . = . (8.48)

. . .0 . .

. ° Mn’1 » ° dn’1 »

°. »n’1 »

. µn’1 2

Mn dn

0 ... 0 µn 2

and it can be e¬ciently solved using the Thomas algorithm (3.53).

A closure condition for system (8.48), which can be useful when the

derivatives f (a) and f (b) are not available, consists of enforcing the con-

tinuity of s3 (x) at x1 and xn’1 . Since the nodes x1 and xn’1 do not

actually contribute in constructing the cubic spline, it is called a not-a-

knot spline, with “active” knots {x0 , x2 , . . . , xn’2 , xn } and interpolating f

at all the nodes {x0 , x1 , x2 , . . . , xn’2 , xn’1 , xn }.

Remark 8.2 (Speci¬c software) Several packages exist for dealing with

interpolating splines. In the case of cubic splines, we mention the command

spline, which uses the not-a-knot condition introduced above, or, in gen-

eral, the spline toolbox of MATLAB [dB90] and the library FITPACK

[Die87a], [Die87b].

8.6 Approximation by Splines 351

A completely di¬erent approach for generating s3 consists of providing

a basis {•i } for the space S3 of cubic splines, whose dimension is equal to

n + 3. We consider here the case in which the n + 3 basis functions •i have

global support in the interval [a, b], referring to Section 8.6.2 for the case

of a basis with local support.

Functions •i , for i, j = 0, . . . , n, are de¬ned through the following inter-

polation constraints

•i (xj ) = δij , •i (x0 ) = •i (xn ) = 0,

and two suitable splines must be added, •n+1 and •n+2 . For instance, if

the spline must satisfy some assigned conditions on the derivative at the

end points, we ask that

•n+1 (xj ) = 0, j = 0, ..., n •n+1 (x0 ) = 1, •n+1 (xn ) = 0,

•n+2 (xj ) = 0, j = 0, ..., n •n+2 (x0 ) = 0, •n+2 (xn ) = 1.

By doing so, the spline takes the form

n

s3 (x) = fi •i (x) + f0 •n+1 (x) + fn •n+2 (x),

i=0

where f0 and fn are two given values. The resulting basis {•i , i = 0, ..., n + 2}

is called a cardinal spline basis and is frequently employed in the numerical

solution of di¬erential or integral equations. Figure 8.9 shows a generic car-

dinal spline, which is computed over a virtually unbounded interval where

the interpolation nodes xj are the integers. The spline changes sign in any

adjacent intervals [xj’1 , xj ] and [xj , xj+1 ] and rapidly decays to zero.

Restricting ourselves to the positive axis, it can be shown (see [SL89])

that the extremant of the function on the interval [xj , xj+1 ] is equal to

the extremant on the interval [xj+1 , xj+2 ] multiplied by a decaying factor

» ∈ (0, 1). In such a way, possible errors arising over an interval are rapidly

damped on the next one, thus ensuring the stability of the algorithm.

Let us summarize the main properties of interpolating cubic splines, re-

ferring to [Sch81] and [dB83] for the proofs and more general results.

Property 8.2 Let f ∈ C 2 ([a, b]), and let s3 be the natural cubic spline

interpolating f . Then

b b

[s3 (x)]2 dx ¤ [f (x)]2 dx, (8.49)

a a

where equality holds if and only if f = s3 .

The above result is known as the minimum norm property and has the

meaning of the minimum energy principle in mechanics. Property (8.49)

352 8. Polynomial Interpolation

0.8

0.6

0.4

0.2

0

’0.2

’4 ’3 ’2 ’1 0 1 2 3 4 5

FIGURE 8.9. Cardinal spline

still holds if conditions on the ¬rst derivative of the spline at the end points

are assigned instead of natural conditions (in such a case, the spline is called

constrained, see Exercise 11).

The cubic interpolating spline sf of a function f ∈ C 2 ([a, b]), with

sf (a) = f (a) and sf (b) = f (b), also satis¬es the following property

b b

[f (x) ’ sf (x)]2 dx ¤ [f (x) ’ s (x)]2 dx, ∀s ∈ S3 .

a a

As far as the error estimate is concerned, the following result holds.

Property 8.3 Let f ∈ C 4 ([a, b]) and ¬x a partition of [a, b] into subinter-

vals of width hi such that h = maxi hi and β = h/ mini hi . Let s3 be the

cubic spline interpolating f . Then

(r)

f (r) ’ s3 ¤ Cr h4’r f (4) ∞, r = 0, 1, 2, 3, (8.50)

∞

with C0 = 5/384, C1 = 1/24, C2 = 3/8 and C3 = (β + β ’1 )/2.

As a consequence, spline s3 and its ¬rst and second order derivatives

uniformly converge to f and to its derivatives, as h tends to zero. The third

order derivative converges as well, provided that β is uniformly bounded.

Example 8.8 Figure 8.10 shows the cubic spline approximating the function in

the Runge™s example, and its ¬rst, second and third order derivatives, on a grid

of 11 equally spaced nodes, while in Table 8.5 the error s3 ’ f ∞ is reported as

a function of h together with the computed order of convergence p. The results

clearly demonstrate that p tends to 4 (the theoretical order) as h tends to zero.

•

8.6 Approximation by Splines 353

h 1 0.5 0.25 0.125 0.0625

s3 ’ f 0.022 0.0032 2.7741e-4 1.5983e-5 9.6343e-7

∞

p “ 2.7881 3.5197 4.1175 4.0522

TABLE 8.5. Experimental interpolation error for Runge™s function using cubic

splines

1 0.8

0.9 (a)

0.6

(b)

0.8

0.4

0.7

0.2

0.6

0.5 0

0.4

’0.2

0.3

’0.4

0.2

’0.6

0.1

0 ’0.8

’5 ’4 ’3 ’2 ’1 0 1 2 3 4 5 ’5 ’4 ’3 ’2 ’1 0 1 2 3 4 5

5

1

4 (d)

(c)

0.5

3

2

0

1

0

’0.5

’1

’1

’2

’3

’1.5

’4

’5

’2

’5 ’4 ’3 ’2 ’1 0 1 2 3 4 5

’5 ’4 ’3 ’2 ’1 0 1 2 3 4 5

FIGURE 8.10. Interpolating spline (a) and its ¬rst (b), second (c) and third (d)

order derivatives (in solid line) for the function of Runge™s example (in dashed

line)

8.6.2 B-splines

Let us go back to splines of a generic degree k, and consider the B-spline

(or bell-spline) basis, referring to divided di¬erences introduced in Section

8.2.1.

De¬nition 8.2 The normalized B-spline Bi,k+1 of degree k relative to the

distinct nodes xi , . . . , xi+k+1 is de¬ned as

Bi,k+1 (x) = (xi+k+1 ’ xi )g[xi , . . . , xi+k+1 ], (8.51)

where

(t ’ x)k if x ¤ t,

g(t) = (t ’ x)k = (8.52)

+

0 otherwise.

354 8. Polynomial Interpolation

Substituting (8.18) into (8.51) yields the following explicit representation

k+1

(xj+i ’ x)k

+

Bi,k+1 (x) = (xi+k+1 ’ xi ) . (8.53)

k+1

j=0

(xi+j ’ xi+l )

l=0

l=j

From (8.53) it turns out that the active nodes of Bi,k+1 (x) are xi , . . . , xi+k+1

and that Bi,k+1 (x) is non null only within the interval [xi , xi+k+1 ].

Actually, it can be proved that it is the unique non null spline of min-

imum support relative to nodes xi , . . . , xi+k+1 [Sch67]. It can also be

(l) (l)

shown that Bi,k+1 (x) ≥ 0 [dB83] and |Bi,k+1 (xi )| = |Bi,k+1 (xi+k+1 )| for

l = 0, . . . , k ’1 [Sch81]. B-splines admit the following recursive formulation

([dB72], [Cox72])

1 if x ∈ [xi , xi+1 ],

Bi,1 (x) =

0 otherwise, (8.54)

x ’ xi xi+k+1 ’ x

Bi+1,k (x), k ≥ 1,

Bi,k+1 (x) = Bi,k (x) +

xi+k ’ xi xi+k+1 ’ xi+1

which is usually preferred to (8.53) when evaluating a B-spline at a given

point.

Remark 8.3 It is possible to de¬ne B-splines even in the case of partially

coincident nodes, by suitably extending the de¬nition of divided di¬erences.

This leads to a new recursive form of Newton divided di¬erences given by

(see for further details [Die93])

±

f [x1 , . . . , xn ] ’ f [x0 , . . . , xn’1 ] if x < x < . . . < x

0 1 n

xn ’ x0

f [x0 , . . . , xn ] =

f (n+1) (x0 )

if x0 = x1 = . . . = xn .

(n + 1)!

Assuming that m (with 1 < m < k + 2) of the k + 2 nodes xi , . . . , xi+k+1

are coincident and equal to », then (8.46) will contain a linear combination

k+1’j

of the functions (» ’ x)+ , for j = 1, . . . , m. As a consequence, the

B-spline can have continuous derivatives at » only up to order k ’ m and,

therefore, it is discontinuous if m = k + 1. It can be checked [Die93] that,

if xi’1 < xi = . . . = xi+k < xi+k+1 , then

±

k

xi+k+1 ’ x

if x ∈ [xi , xi+k+1 ],

xi+k+1 ’ xi

Bi,k+1 (x) =

0 otherwise,

8.6 Approximation by Splines 355

while for xi < xi+1 = . . . = xi+k+1 < xi+k+2

±

k

x ’ xi

if x ∈ [xi , xi+k+1 ],

xi+k+1 ’ xi

Bi,k+1 (x) =

0 otherwise.

Combining these formulae with the recursive relation (8.54) allows for con-

structing B-splines with coincident nodes.

Example 8.9 Let us examine the special case of cubic B-splines on equally

spaced nodes xi+1 = xi + h for i = 0, ..., n ’ 1. Equation (8.53) becomes

6h3 Bi,4 (x) =

±

(x ’ xi ) , if x ∈ [xi , xi+1 ],

3

3

h + 3h2 (x ’ x ) + 3h(x ’ x )2 ’ 3(x ’ x )3 , if x ∈ [x , x ],

i+1 i+1 i+1 i+1 i+2

h3 + 3h2 (xi+3 ’ x) + 3h(xi+3 ’ x)2 ’ 3(xi+3 ’ x)3 , if x ∈ [xi+2 , xi+3 ],

(xi+4 ’ x)3 ,

if x ∈ [xi+3 , xi+4 ],

0 otherwise.

In Figure 8.11 the graph of Bi,4 is shown in the case of distinct nodes and of

•

partially coincident nodes.

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

’2 ’1 0 1 2

FIGURE 8.11. B-spline with distinct nodes (in solid line) and with three coin-

cident nodes at the origin (in dashed line). Notice the discontinuity of the ¬rst

derivative

Given n + 1 distinct nodes xj , j = 0, . . . , n, n ’ k linearly independent

B-splines of degree k can be constructed, though 2k degrees of freedom are

356 8. Polynomial Interpolation

still available to generate a basis for Sk . One way of proceeding consists of

introducing 2k ¬ctitious nodes

x’k ¤ x’k+1 ¤ . . . ¤ x’1 ¤ x0 = a,

(8.55)

b = xn ¤ xn+1 ¤ . . . ¤ xn+k

which the B-splines Bi,k+1 , with i = ’k, . . . , ’1 and i = n ’ k, . . . , n ’ 1,

are associated with. By doing so, any spline sk ∈ Sk can be uniquely written

as

n’1

sk (x) = ci Bi,k+1 (x). (8.56)

i=’k

The real numbers ci are the B-spline coe¬cients of sk . Nodes (8.55) are

usually chosen as coincident or periodic.

1. Coincident: this choice is suitable for enforcing the values attained

by a spline at the end points of its de¬nition interval. In such a case,

indeed, thanks to Remark 8.3 about B-splines with coincident nodes,

we get

sk (a) = c’k , sk (b) = cn’1 . (8.57)

2. Periodic, that is

x’i = xn’i ’ b + a, xi+n = xi + b ’ a, i = 1, . . . , k.

This choice is useful if the periodicity conditions (8.44) have to be

imposed.

Remark 8.4 (Inserting nodes) Using B-splines instead of cardinal spli-

nes is advantageous when handling, with a reduced computational e¬ort, a

given con¬guration of nodes for which a spline sk is known. In particular,

assume that the coe¬cients ci of sk (in form (8.56)) are available over the

nodes x’k , x’k+1 , . . . , xn+k , and that we wish to add to these a new node

x.

The spline sk ∈ Sk , de¬ned over the new set of nodes, admits the follow-

˜

ing representation with respect to a new B-spline basis Bi,k+1

n’1

sk (x) = di Bi,k+1 (x).

i=’k

The new coe¬cients di can be computed starting from the known coe¬-

cients ci using the following algorithm [Boe80]:

8.7 Splines in Parametric Form 357

let x ∈ [xj , xj+1 ); then, construct a new set of nodes {yi } such that

yi = xi for i = ’k, . . . , j, yj+1 = x,

yi = xi’1 for i = j + 2, . . . , n + k + 1;

de¬ne

±

for i = ’k, . . . , j ’ k,

1

y

j+1 ’ yi

for i = j ’ k + 1, . . . , j,

ωi =

yi+k+1 ’ yi

0 for i = j + 1, . . . , n;

compute

di = ωi ci + (1 ’ ωi )ci i = ’k, ..., n ’ 1.

This algorithm has good stability properties and can be generalized to the

case where more than one node is inserted at the same time (see [Die93]).

8.7 Splines in Parametric Form

Using interpolating splines presents the following two drawbacks:

1. the resulting approximation is of good quality only if the function

f does not exhibit large derivatives (in particular, we require that

|f (x)| < 1 for every x). Otherwise, oscillating behaviors may arise

in the spline, as demonstrated by the example considered in Figure

8.12 which shows, in solid line, the cubic interpolating spline over the

following set of data (from [SL89])

xi 8.125 8.4 9 9.845 9.6 9.959 10.166 10.2

fi 0.0774 0.099 0.28 0.6 0.708 1.3 1.8 2.177

2. sk depends on the choice of the coordinate system. In fact, performing

a clockwise rotation of 36 degrees of the coordinate system in the

above example, would lead to the spline without spurious oscillations

reported in the boxed frame in Figure 8.12.

All the interpolation procedures considered so far depend on the cho-

sen Cartesian reference system, which is a negative feature if the

spline is used for a graphical representation of a given ¬gure (for in-

stance, an ellipse). Indeed, we would like such a representation to

be independent of the reference system, that is, to have a geometric

invariance property.

358 8. Polynomial Interpolation

2.5

’4.2

’4.4

2

’4.6

’4.8

1.5

’5

’5.2

6 7 8 9 10

1

0.5

0

8 8.5 9 9.5 10 10.5

FIGURE 8.12. Geometric noninvariance for an interpolating cubic spline s3 : the

set of data for s3 in the boxed frame is the same as in the main ¬gure, rotated

by 36 degrees. The rotation diminishes the slope of the interpolated curve and

eliminates any oscillation from s3 . Notice that resorting to a parametric spline

(dashed line) removes the oscillations in s3 without any rotation of the reference

system

A solution is provided by parametric splines, in which any component of the

curve, written in parametric form, is approximated by a spline function.

Consider a plane curve in parametric form P(t) = (x(t), y(t)), with t ∈

[0, T ], then take the set of the points in the plane of coordinates Pi =

(xi , yi ), for i = 0, . . . , n, and introduce a partition onto [0, T ]: 0 = t0 <

t 1 < . . . < tn = T .

Using the two sets of values {ti , xi } and {ti , yi } as interpolation data,

we obtain the two splines sk,x and sk,y , with respect to the independent

variable t, that interpolate x(t) and y(t), respectively. The parametric curve

Sk (t) = (sk,x (t), sk,y (t)) is called the parametric spline. Obviously, di¬erent

parameterizations of the interval [0, T ] yield di¬erent splines (see Figure

8.13).

A reasonable choice of the parameterization makes use of the length of

each segment Pi’1 Pi ,

(xi ’ xi’1 )2 + (yi ’ yi’1 )2 , i = 1, . . . , n.

li =

i

Setting t0 = 0 and ti = k=1 lk for i = 1, . . . , n, every ti represents the

cumulative length of the piecewise line that joins the points from P0 to

Pi . This function is called the cumulative length spline and approximates

satisfactorily even those curves with large curvature. Moreover, it can also

be proved (see [SL89]) that it is geometrically invariant.

Program 68 implements the construction of cumulative parametric cu-

bic splines in two dimensions (it can be easily generalized to the three-

8.7 Splines in Parametric Form 359

4

2

0

’2

’4

’2 0 2 4 6

FIGURE 8.13. Parametric splines for a spiral-like node distribution. The spline

of cumulative length is drawn in solid line

dimensional case). Composite parametric splines can be generated as well

by enforcing suitable continuity conditions (see [SL89]).

Program 68 - par spline : Parametric splines

function [xi,yi] = par spline (x, y)

t (1) = 0;

for i = 1:length (x)-1

t (i+1) = t (i) + sqrt ( (x(i+1)-x(i))ˆ2 + (y(i+1)-y(i))ˆ2 );

end

z = [t(1):(t(length(t))-t(1))/100:t(length(t))];

xi = spline (t,x,z);

yi = spline (t,y,z);

8.7.1 B´zier Curves and Parametric B-splines

e

The B´zier curves and parametric B-splines are widely employed in graph-

e

ical applications, where the nodes™ locations might be a¬ected by some

uncertainty.

Let P0 , P1 , . . . , Pn be n + 1 points ordered in the plane. The oriented

polygon formed by them is called the characteristic polygon or B´zier poly-

e

gon. Let us introduce the Bernstein polynomials over the interval [0, 1]

de¬ned as

n!

n

tk (1 ’ t)n’k = tk (1 ’ t)n’k ,

bn,k (t) =

k!(n ’ k)!

k

360 8. Polynomial Interpolation

for n = 0, 1, . . . and k = 0, . . . , n. They can be obtained by the following

recursive formula

bn,0 (t) = (1 ’ t)n

bn,k (t) = (1 ’ t)bn’1,k (t) + tbn’1,k’1 (t), k = 1, . . . , n, t ∈ [0, 1].

It is easily seen that bn,k ∈ Pn , for k = 0, . . . , n. Also, {bn,k , k = 0, . . . , n}

provides a basis for Pn . The B´zier curve is de¬ned as follows

e

n

0 ¤ t ¤ 1.

Bn (P0 , P1 , . . . , Pn , t) = Pk bn,k (t), (8.58)

k=0

This expression can be regarded as a weighted average of the points Pk ,

with weights bn,k (t).

The B´zier curves can also be obtained by a pure geometric approach

e

starting from the characteristic polygon. Indeed, for any ¬xed t ∈ [0, 1],

we de¬ne Pi,1 (t) = (1 ’ t)Pi + tPi+1 for i = 0, . . . , n ’ 1 and, for t

¬xed, the piecewise line that joins the new nodes Pi,1 (t) forms a polygon

of n ’ 1 edges. We can now repeat the procedure by generating the new

vertices Pi,2 (t) (i = 0, . . . , n ’ 2), and terminating as soon as the polygon