t

t

t

yi+1 i+1

t

yi i

t

1

' hi E

t0

y0

x0 xi xi+1 xn

Figure 19.1 Cubic spline interpolation of n + 1 points in n intervals. The i-th

interval has a width of hi and runs from xi to xi+1 . First and second derivatives

are continuous at all nodes.

further coordinates in higher-dimensional spaces) are polynomial functions

of one (or more, in higher-dimensional spaces) parameter(s). To this cate-

gory belong B-splines and Bezier curves and surfaces. B-splines are treated

in Section 19.7; they are often used to construct smooth surfaces and multi-

dimensional interpolations, for example in the smooth-particle mesh-Ewald

(SPME) method for computing long-range interactions in periodic systems

(see Section 13.10.6 in Chapter 13 on page 373). For Bezier curves we re-

fer to the literature. For many more details, algorithms, programs in C,

and variations on this theme, the reader is referred to Engeln-M¨llges and

u

Uhlig (1996). Programs including higher-dimensional cases are also given

in Sp¨th (1973). A detailed textbook is de Boor (1978), while Numerical

a

Recipes (Press et al., 1992) is a practical reference for cubic splines.

In the following sections the spline methods are introduced, with emphasis

on the very useful cubic splines. A general Python program to compute one-

dimensional cubic splines is provided in Section 19.6. Section 19.7 treats

B-splines.

526 Splines for everything

19.2 Cubic splines through points

Consider n + 1 points xi , yi , i = 0, . . . , n (Fig. 19.1). We wish to construct a

function f (x) consisting of piece-wise functions fi (x ’ xi ), i = 0, . . . , n ’ 1,

de¬ned on the interval [0, hi = xi+1 ’ xi ], with the following properties:

(i) each fi (ξ) is a polynomial of the third degree;

fi (0) = yi , i = 0, . . . , n ’ 1;

(ii)

(iii) fn’1 (hn’1 ) = yn ;

fi (hi ) = fi+1 (0), i = 0, . . . , n ’ 2;

(iv)

fi (hi ) = fi+1 (0), i = 0, . . . , n ’ 2;

(v)

fi (hi ) = fi+1 (0), i = 0, . . . , n ’ 2.

(vi)

There are n intervals with 4n parameters to describe the n piecewise funct-

ions; the properties provide (n+1)+3(n’1) = 4n’2 equations. In order to

solve for the unknown parameters, two extra conditions are needed. These

are provided by speci¬cation of two further derivatives at the end nodes,

choosing from the ¬rst, second or third derivatives of either or both end

nodes.2 In the case of periodic cubic splines, with period xn ’ x0 , for which

yn = y0 , it is speci¬ed that the ¬rst and second derivative at the ¬rst point

are equal to the ¬rst and second derivative at the last point. If it is speci¬ed

that the second derivatives at both end points vanish, natural splines result.

It is clear that quadratic splines will be obtained when the function and its

¬rst derivative are continuous (one extra condition being required), and that

quartic and higher-order splines can be de¬ned as well.

The function in the ith interval is given by

(0 ¤ ξ ¤ hi ).

f (xi + ξ) = fi (ξ) = fi + gi ξ + pi ξ 2 + qi ξ 3 (19.3)

Here fi , gi , pi and qi are the four parameters to be determined from the

conditions given above. We see immediately that the values of the function

and its ¬rst derivative in point xi are

fi = f (xi ) = yi , (19.4)

and

gi = f (xi ) = fi , (19.5)

respectively, while 2pi is the second derivative at xi and 6qi is the third

derivative, which is constant throughout the i-th interval, and discontinuous

at the nodes.

2 It is wise to specify one condition at each end node; specifying two conditions at one end node

may lead to accumulating errors giving erratic oscillations.

19.2 Cubic splines through points 527

The continuity conditions lead to the following equations:

fi (hi ) = fi + gi hi + pi h2 + qi h3 = fi+1 , (19.6)

fi (hi ) = gi + 2pi hi + 3qi h2 = gi+1 = fi+1 , (19.7)

i

fi (hi ) = 2pi + 6qi hi = 2pi+1 = fi+1 , , (19.8)

valid for i = 0, . . . , n ’ 2.

Four parameters su¬ce for the reconstruction of f (x) in each interval. It

is convenient to use the function values and the ¬rst derivatives in the nodes

for reconstruction. Thus, for the i-th interval we obtain by solving pi and

qi from (19.6) and (19.7):

1 3(fi+1 ’ fi )

’ 2gi ’ gi+1 ,

pi = (19.9)

hi hi

1 ’2(fi+1 ’ fi )

qi = 2 + gi + gi+1 . (19.10)

hi

hi

Any value of f (x) within the i-th interval is easily found from (19.3) when

both the function and derivative values are given at the nodes of that inter-

val. So the interpolation is practically solved when the ¬rst derivatives are

known, and the task of constructing a cubic spline from a set of data points

is reduced to the task of ¬nding the ¬rst derivatives in all points (including

the end nodes).3

The continuity of the second derivatives is given by (19.8). Written in

terms of f and g, this condition gives n ’ 1 equations for n + 1 unknowns

g0 , . . . , gn :

1 1 1 1

gi + 2 + gi+1 + gi+2

hi hi hi+1 hi+1

fi+1 ’ fi fi+2 ’ fi+1

, i = 0, . . . , n ’ 2.

=3 + (19.11)

h2 h2

i i+1

These equations must be augmented by the two additional conditions, and

“ depending on the exact conditions “ will lead to a matrix equation with a

symmetric, tridiagonal matrix. For example, if the ¬rst derivatives g0 and

gn at the end points are known, they can be removed from the unknown

vector and lead to a matrix equation with g = [g1 , . . . gn’1 ] as the unknown

vector:

Ag = b, (19.12)

3 We could have chosen to solve for the second derivatives at all nodes, as is often done in the

literature, since the functions in each interval are also easily constructed from knowledge of

the values of the function and its second derivatives at the nodes. The computational e¬ort is

similar in both cases.

528 Splines for everything

⎛ ⎞⎛ ⎞ ⎛ ⎞

d1 s1 b1 ’ g0 /h0

g1

⎜ s1 ⎟⎜ ⎟⎜ ⎟

d2 s2

⎜ ⎟⎜ g2 b2

⎟⎜ ⎟

⎜ ⎟⎜ ⎟⎜ ⎟

.. .. ..