p=2 N0

1/2 1/2

(3)

p=3 N0

2/3

1/6 1/6 (4)

p=4 N0

11/24 11/24

(5)

1/24 1/24

p=5 N0

t

’1 0 1 2 3 4 5 6

(p)

Figure 19.9 The ¬rst ¬ve normalized cardinal B-spline functions N0 (t).

The convolution theorem implies that the Fourier transform of N (p) is the

product of the Fourier transforms of N (q) and N (p’q) . Since the Fourier

transform of the unit pulse N (1) is proportional to sin 1 ω/ω (ω is the conju-

2

gate variable in reciprocal space),12 the Fourier transform of N (p) is propor-

tional to the p-th power of that function, disregarding a phase factor due

to the shifted origin. This transform decays faster in reciprocal space for

higher-order splines.

A B-spline function of order p consists of piecewise polynomials of degree

p ’ 1 and is p ’ 2 times continuously di¬erentiable; thus N (4) (t) is twice

di¬erentiable, which makes it (and higher-order functions) suitable to ex-

press energy functions: the ¬rst derivative (needed for forces) is su¬ciently

smooth as it has a continuous derivative. For the derivatives the following

12 This applies to the central function; the Fourier transform of the function N (1) contains a

phase factor exp( 1 iω) because the origin is shifted with respect to the symmetry axis.

2

19.7 B-splines 551

relation holds:

(p)

dN0 (t) (p’1) (p’1)

(t) ’ N1

= N0 (t). (19.74)

dt

B-spline curves are not interpolating, as the control points are not equal

to the nodal values. Each nodal value is composed of nearby control points

according to (19.62). One may construct interpolating curves by deriving

the control points from a corresponding set of equations. For closed, i.e.,

periodic, B-spline curves, function values are periodic on [0, n) and there

are exactly n nodes. One therefore needs n known values to solve for the

coe¬cients ci . When the known values are the nodal function values, an

interpolating periodic spline is obtained:

n’1

(p)

f (t) = ci Ni+kn (t). (19.75)

i=0 k∈Z

The sum over k takes care of the periodicity, but only a few values are

included because of the compact support of the functions. For example, the

interpolating fourth-order B-spline for a sine function sin 1 πt on the periodic

2

interval [0, 4) (Fig. 19.10), given the nodal values (0, 1, 0, ’1) at the nodes

(0,1,2,3) is

3 3

f (t) = (N’1 + N3 ) ’ (N1 + N5 ). (19.76)

2 2

These coe¬cients are the solution of the following set of equations:

⎛ ⎞⎛ ⎞⎛ ⎞

0141 c0 0

1 ⎜ 1 0 1 4 ⎟ ⎜ c1 ⎟ ⎜ 1⎟

⎜ ⎟⎜ ⎟⎜ ⎟,

⎝ 4 1 0 1 ⎠ ⎝ c2 ⎠ = ⎝ (19.77)

0⎠

6

’1

1410 c3

and have the values (0, ’3/2, 0, 3/2). The supporting B-spline functions are

depicted in Fig. 19.10. Not surprising this B-spline is identical to the cubic

interpolating spline shown in Fig. 19.2 on page 531.

The recurrence relation (19.67) allows the following simple algorithm to

compute spline values (Griebel et al., 2003). The value of p must be an

integer not smaller than 2. While a recurrent program is elegant, a non-

recurrent version would be more e¬cient.

python program 19.5 Bspline

Computes the value of the cardinal B-spline function N (p) (t) of real t and

integer order p ¤ 2.

552 Splines for everything

sin (πt/2)

N’3(t) c1 = ’3/2

N’2(t) c2 = 0

N’1(t) c3 = 3/2

c0 = 0

N0(t)

c1 = ’3/2

N1(t)

c2 = 0

N2(t)

c3 = 3/2

N3(t)

’2 0 2 4 6 t

Figure 19.10 The sine function, periodic on the interval [0, 4), approximated by a

fourth-order B-spline.

01 def Bspline(t,p):

02 if (t<=0. or t>=p): return 0.

03 if (p==2): return 1.-abs(t-1.)

04 return (t*Bspline(t,p-1) + (p-t)*Bspline(t-1.,p-1))/(p-1.)

Comments

Recurrence is pursued until the order p = 2 is reached, which is the “hat-function” speci¬ed in

line 2.

For non-periodic (open) B-splines that are to be ¬tted on a function known

at n + 1 points t = 0, . . . , n, one needs coe¬cients for all the B-spline funct-

ions that share their support with the interval [0, n), which are n + p ’ 1

in number. This means that, in addition to the n + 1 nodal values, (p ’ 2)

extra conditions must be speci¬ed. These can be points outside the interval

or ¬rst or second derivative speci¬cations. Instead of function values at the

nodes, values at arbitrary (e.g., random) arguments may be speci¬ed, as

long as they provide a su¬cient number of independent data to solve for the

coe¬cients.

While in general a set of equations must be solved to ¬nd the spline

19.7 B-splines 553

coe¬cients, there is an exception for exponential functions. Functions of

the type f (t) = z t , among which the important class of exponentials f (t) =

exp(iωt), can be approximated by a B-spline series

∞

f (t) ≈ (p)

cn Nn (t). (19.78)

n=’∞

These are called Euler exponential splines (Schoenberg, 1973; Silliman,

1974). We wish to ¬nd the coe¬cients cn such that the spline series equals