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