<<

. 14
( 26)



>>

0.4 1.9320 1.83 -2.46
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

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

<<

. 14
( 26)



>>