<< . .

. 50
( : 95)



. . >>

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

<< . .

. 50
( : 95)



. . >>