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