We now consider an example: the construction of the distribution function

of 50 samples from a normal distribution. Figure 19.5 shows the cumulative

distribution of the data. In Fig. 19.6 the cumulative distribution is put into

bins of width 0.5, and a ¬tting spline is constructed with such a weight that

19.5 Splines for tabulation 539

1

w = 0.00035

χ2/n = 1.06

0.8 curv = 0.101

0.6

0.4

0.2

0

’3 ’2 ’1 0 1 2 3

Figure 19.6 A ¬tting spline approximating the cumulative distribution.

0.3

0.2

0.1

0

’3 ’2 ’1 0 1 2 3

Figure 19.7 Distribution constructed (by a spline) to the ¬rst derivative of the

cumulative ¬t.

χ2 /n = 1. Finally, Fig. 19.7 shows a spline ¬t through the derivative values

of Fig. 19.6. The distribution is still some way from being Gaussian, but it

is the smoothest (cumulative) distribution compatible with the data.

19.5 Splines for tabulation

Splines can be used to recover forces and energies from tabulated values.

Assume we have a centrosymmetric potential V (r) between a pair of par-

540 Splines for everything

ticles, de¬ned for 0 < r ¤ rc , with rc a cuto¬ distance beyond which the

potential is identically zero. The particles have cartesian coordinates r i and

r j ; we use the convention r = r i ’ r j and r = |r|. The force F = F i on

particle i, which is equal and opposite to the force acting on particle j, is

given by

dV r

F (r) = ’ . (19.50)

dr r

We wish to tabulate values of V such that:

• V and dV /dr are recovered simultaneously from the table;

• the derivative of the force is continuous “ this is accomplished by using

cubic spline interpolation for V ;

• the interpolation should be fast “ this requires storage of not only the

values of the potential, but (at least) also of its derivative, allowing local

cubic interpolation;

• a search procedure to ¬nd the table index for an arbitrary argument should

be avoided “ the table index should follow directly from the argument;

• the computationally expensive square root evaluation should be avoided

“ this means that table indices should be derived from r2 rather than r.

These requirements lead to the following procedures.

Construction of the table

(i) Decide on the number of table entries n + 1 (the table index runs

from 0 to n).

(ii) Decide on the cuto¬ radius rc .

(iii) De¬ne a scaled variable x, proportional to r2 , running from 0 to n.

nr2

x= 2 . (19.51)

rc

The table index i runs from 0 to n.

(iv) Construct a table with function values fi and derivatives gi . The

values should be the required potential and its derivative. De¬ning

i

r i = rc , (19.52)

n

the values are

fi = V (ri ), (19.53)

2

dV rc dV

gi = = (ri ). (19.54)

dx 2nri dr

19.5 Splines for tabulation 541

The values gi can also be determined from a spline ¬t to fi . In

that case the second derivatives are strictly continuous; if the gi ™s

are determined from analytical derivatives, the interpolated values

approach the analytical function more accurately, but there will be

slight discontinuities in the second derivatives at the nodes.

For fast recovery from memory fi and gi should be contiguously

stored, for example in a two-dimensional array, rather than in sepa-

rate arrays.8

Recovery from the table

We wish to obtain V (r) and F (r) for a given value of r2 , where

r2 = (xi ’ xj )2 + (yi ’ yj )2 + (zi ’ zj )2 . (19.55)

(i) Determine x = nr2 /rc , i = int (x) and ξ = x ’ i.

2

(ii) Fetch fi , gi , fi+1 , gi+1 from the table.

(iii) Determine

p = 3fi+1 ’ 3fi ’ 2gi ’ gi+1 , (19.56)

q = ’2fi+1 + 2fi + gi + gi+1 . (19.57)

(iv) Determine V and dV /dx:

V = fi + ξgi + ξ 2 p + ξ 3 q, (19.58)

dV

= gi + 2ξp + 3ξ 2 q. (19.59)

dx

(v) Reconstruct the force:

2n dV

F =’ (19.60)

r.

2

rc dx

We note that the determination of p and q takes six ¬‚oating point operations

(¬‚ops),9 while the actual interpolation takes ten ¬‚ops. One may choose to

store also pi and qi , thus saving 6 out of 16 ¬‚ops at the expense of doubling

the required memory. The total number of ¬‚ops is comparable to that

required to compute Lennard“Jones potential and forces; for almost any

other kind of interaction tables are more e¬cient than analytical evaluation.

Using r2 instead of r causes a rather inappropriate mapping of r onto the

table index, the spacing in r being much ¬ner for r close to rc than for values

close to the potential minimum, where high accuracy is required. A more

8 This requirement depends on the architecture of the computer used; it is certainly true for

vector processors that will fetch several words from memory at once. Processors that will hold

all tables in cache memory are less demanding in this respect.

For example, as follows: u = fi+1 ’ fi ; q = ’2u + gi + gi+1 ; p = ’q + u ’ gi .

9

542 Splines for everything

x/n

1

0.8

b

0.6