Piecewise polynomial interpolation, particularly with cubic splines

24 Pages·2000·0.246 MB·English

by  Walsh B.

PIECEWISE POLYNOMIAL INTERPOLATION,
PARTICULARLY WITH CUBIC SPLINES
1. Setup. Everything we do with piecewise-polynomial approximation will be set in the following
more-or-less standard framework. We are given (n+1) nodes, usually numbered from left to right,
a=x <x <(cid:1)(cid:1)(cid:1)<x =b
0 1 n
which divide the interval J =[x ;x ] into n cells, the i-thbeing J =[x ;x ] oflength h =x −x >0.
0 n i i i+1 i i+1 i
Implicitlywearegivenafunctionf(x)de(cid:12)nedonJ. Weseeknpolynomialsfp (x):i=0; :::; n−1gwhich
i
will de(cid:12)ne a function \near" f(x) by de(cid:12)ning the \nearby" function as equal to p (x) on the cell [x ;x ].
i i i+1
Thus there is no \single analytic de(cid:12)nition" for the approximating function|one has to know the cell in
which x lies in order to evaluate the approximating function at x. All the polynomials p (x) will have the
i
same degree, however|indeed, in the most important case they will all be cubics|and the \nearness" will
include the requirement that the piecewise polynomial function they de(cid:12)ne interpolate f(x) on the nodes,
so that
p0(x0)=f(x0)=f0; p0(x1)=p1(x1)=f(x1)=f1; :::; pn−1(xn−1)=f(xn−1)=fn−1;
and moreover, pn−1(xn)=f(xn)=fn :
Thus each polynomial satis(cid:12)es two interpolation conditions, and in addition to those there may be \end
conditions" atx =a andx =b and conditions that link the \pieces" atthe nodes. When we need a name
0 n
forthefunctionde(cid:12)nedonallofJ byusingthepolynomialpieces,wecallits(x);thepiecede(cid:12)nedonJ will
i
becalledp (x). Whenweneedasingleestimateforthelengthsofthecells,weuseh=max fjh j:0(cid:20)i<ng.
i i i
We shall tacitly assume that the function f(x) has as much di(cid:11)erentiability as it needs to make our
theorems make sense. In particular, functions to be interpolated by cubic splines will usually be assumed
to be twice-continuously di(cid:11)erentiable at least, and sometimes \twice" will be replaced by \four-" or \(cid:12)ve-
times." People who know what these things are will realize that in some situations it would su(cid:14)ce that the
functions have distributional derivatives that belong to L2(J) (buzzword = Sobol¨ev space).
fRemark: The convention used above to number intervals J and increments h numbers the intervals
i i
by their left-hand endpoints. This convention conforms to Atkinson’s book, but you should note that Stoer
& Bulirsch [for example] number by r. h. endpoints. I feel that numbering by r. h. endpoints is better for
codingpurposes,butthe reasonforthis feelingisprobablythattheearlyversionsofFortranwereunwilling
to admit zero as the value of an index. Possibly C programming would be easier with the intervals indexed
from 0 to n−1. In all cases, though, I see no reason gratuitously to (cid:12)ght the textbook on this matter.g
2. Lower-degree and Uncoupled Examples. Here are some simple cases of piecewise polynomial
approximation. If all the p (x)’s are linear, then just the requirement that p (x) interpolate the values of
i i
f(x) at the endpoints of J uniquely determines the polynomial p (x): there is only one linear interpolant
i i
between the two values at x and x . The function s(x) is thus uniquely determined. An estimate of its
i i+1
deviation from f(x) can be made on an interval-by-intervalbasis: in J we can write
i
f00((cid:24) )
f(x)=p (x)+ x (x−x )(x−x ):
i i i+1
2!
h2
The maximum value of j(x−x )(x−x )j is i, attained at the midpoint of J . We thus have a global,
i i+1 i
4
uniform error estimate
M h2
jf(x)−s(x)j(cid:20) 2 (cid:1)
2 4
where M = maxfjf00(x)j : x 2 Jg. The object we have constructed is called the piecewise-linear inter-
2
polant of f(x), for obvious reasons.
The next case would be to take all the p (x)’s quadratic, and a natural (cid:12)rst hope would be that we
i
might be able to match up the (cid:12)rst derivative of s(x) with that of f(x) at the nodes. That hope does not
1 642:573:01Fall2000
long survive the following count. If all the p (x)’s are quadratic, then there are 3n coe(cid:14)cients required in
i
order completely to determine s(x). The requirement that p (x) take the same values as f(x) at each node
i
imposes 2 conditions oneachp (x); there aren intervals,so that requiresthe coe(cid:14)cients to satisfy 2n linear
i
equations. If wesimply wantto haves0(x) continuous ateachofthe interiornodesx1; :::; xn−1,weimpose
an additional (n−1) equations
p0(x )=p0 (x ); i=0; :::; n−2
i i+1 i+1 i+1
thatthecoe(cid:14)cientsmustsatisfy: wehaveimposed2n+(n−1)=3n−1linearconditionsonthecoe(cid:14)cients,
and we would expect that \there is only one degree of freedom remaining." Certainly an attempt to make
s0(x ) = f0(x ) at each node would add (n+1) equations, bringing the count to 4n, and it does not seem
i i
reasonable to expect such a heavily overdetermined a set of equations in 3n unknowns to have a solution at
all.
In fact, it is easy to see what the remaining degree of freedom is. Given an interval [(cid:11);(cid:12)] and some
number γ, there is a unique quadratic polynomial q(x) that interpolates f(x) at the endpoints and has
q0((cid:11))=γ; in Newton form that polynomial is
f[(cid:11);(cid:12)]−γ
q(x)=f((cid:11))+γ(cid:1)(x−(cid:11))+ (cid:1)(x−(cid:11))2
(cid:12)−(cid:11)
as the reader may see by constructing a divided-di(cid:11)erence table. One can now compute
q0((cid:12))=2(cid:1)f[(cid:11);(cid:12)]−γ
directly. So if in piecewise-quadratic interpolation we specify a value for p0(x ) = s0(x ), then p0(x ) is
0 0 0 0 1
determined, and if s0(x) is to be continuous at x then p0(x ) is determined ::: and so on down the line;
1 1 1
specifying the derivative at x has taken away the one degree of freedom left in the choice of s(x). fTake
0
[(cid:11);(cid:12)] successively to equal [x ;x ], [x ;x ], ::: .g The coe(cid:14)cients of the p (x)’s are \coupled" at the nodes
0 1 1 2 i
by the requirement that s0(x) be continuous, but they’re not coupled in a very complicated way: one can
see the single choice that one makes in orderto determine the entire cubic spline uniquely movefrom left to
right along the interval as one computes the successive p (x)’s in Newton form (which incidentally happens
i
to be the same as Taylor form in this case).
Although a piecewise quadratic can be made C1 and it is more complicated than a piecewise linear
function, it is not clear that it approximates f(x) any better; in fact, by thinking of an extremely inept
choice of s0(x ) you can easily convince yourself that it may be a worse approximation! fSuppose f(x) is
0
a constant function: if s0(x )6=0 then the piecewise-quadratic will oscillate around the constant, while the
0
piecewise-linear would equal the constant. Of course that was a loaded example|but suppose f(x) wiggles
a little in J but is constant in [x ;x ]. Then even matching s0(x ) and f0(x ) may not help.g If we are to
0 1 n 0 0
have s0(x) behave nicely with respect to f0(x)|in addition to being continuous|it appears that we shall
need more freedom than the one degree that piecewise-quadratic functions would furnish. So on we go to
piecewise-cubic interpolation.
If all the p (x) are to be cubics, then it will require 4n coe(cid:14)cients to determine all the polynomials.
i
The requirementthatthe piecewisefunction they de(cid:12)ne be aninterpolantoff(x)imposes 2(n−1)+2=2n
conditions, all of which can be expressed linearly in the coe(cid:14)cients and which look as if they should be
linearly independent. In piecewise-cubic Hermite interpolation we require that the (cid:12)rst derivative of
thei-thinterpolatingpolynomialinterpolatethe derivativeoff(x)atthenodes. Thisimposes2nadditional
conditions and determines the p (x)’s uniquely. Indeed, (cid:12)nding the p (x)’s then becomes a purely \local"
i i
problem: one (cid:12)nds the cubic Hermite interpolant of f(x) on [x ;x ] and that is p (x); no interaction
i i+1 i
between cells takes place.
The whole situation is exactly analogous to piecewise-linear interpolation; all that happens is that the
(cid:12)rst derivative is interpolated along with the function. To see what happens on a single interval, look at
Atkinson’streatmentofHermiteinterpolationonasingleinterval[a;b],whichistheexampleonpp.162{163.
2 642:573:01Fall2000
Note that here we do have a good error estimate: applying Atkinson’s formula (3.6.15) piece-by-piece gives
us
M
jf(x)−s(x)j(cid:20) 4 (cid:1)h4
384
where M =maxfjf(4)(x)j:x (cid:20)x(cid:20)x g. The error !0 as h!0 twice as fast as the error for piecewise-
4 0 n
1
linear interpolation did, and it is nice to have the factor (cid:25) 10−2:6 in there to help make the error
384
smallertoo. Ingeneral,however,one willhavep00 (x )6=p00(x ), i=1; :::; n−1;the (cid:12)rstderivative ofthe
i−1 i i i
interpolant is continuous, but there will usually be jumps in its second derivative (which is not well-de(cid:12)ned
at the nodes).
The de(cid:12)nition of cubic splines is made (at least partly) in order to avoid those discontinuities: the
idea is that we will force the piecewise cubic to be an interpolant of f(x) but will force the equalities
p00 (x )=p00(x ) i=1; :::; n−1
i−1 i i i
so as to get a continuous second derivative. That already gives 2n+(n−1) = 3n−1 conditions on the
coe(cid:14)cients,so thereis notenoughfreedomremainingto force the (cid:12)rstderivativesto agreewithf0(x) atthe
nodes. What one can do is to force the (cid:12)rst derivatives of the cubic pieces to agree with each other at the
nodes:
p0 (x )=p0(x ) i=1; :::; n−1:
i−1 i i i
This gives n−1 additional conditions, bringing the count up to (3n−1)+(n−1) = 4n−2. While the
common values fb : i = 1; :::; n−1g of the (cid:12)rst derivatives will probably not satisfy b = f0(x ), we may
i i i
hope that the di(cid:11)erence will be fairly small if f(x) is su(cid:14)ciently smooth and if the fh gn−1 are su(cid:14)ciently
i i=0
small. Oneexpectsthat2additionalconditionsareneededifallthepiecesofthepiecewisecubicinterpolant
are to be determined uniquely. While there is no canonicalchoice of these conditions, three popular choices
are
(1) The Complete Cubic Spline Interpolant, produced by choosing
0 0 0 0
p (x )=f (x ) and p (x )=f (x );
0 0 0 n−1 n n
\the (cid:12)rst derivative is correct at the endpoints."
(2) The Natural Cubic Spline Interpolant, produced by choosing the fb g so that
i
00 00
p (x )=0 and p (x )=0;
0 0 n−1 n
\no curvature is present at the endpoints." It can be argued that this is not at all natural, since f00(x)
neednotbezeroattheendpoints,andindeedtheerrorestimatesforthisinterpolantarenotassmallas
thoseforthe completespline. (While itisnotobviousthatthis conditioncanbe producedbyasuitable
choice of the fb g, we shall see later that it can be.)
i
(3) The Periodic Cubic Spline Interpolant, produced by choosing the fb g so that
i
0 0 00 00
p (x )=p (x ) and p (x )=p (x );
0 0 n n 0 0 n n
it is most naturalto use these conditions with a function for which f(x )=f(x ), so that the function
0 n
f(x) \wants to have a periodic extension" to the real line with period cells that are the translates of J.
fOf course a function can have periodic derivatives without itself being periodic; the interested reader
may want to contemplate the situation further.g
Atkinson(pp. 166{176)dealswith the naturalcubic spline interpolant,(2)above. Thesenotes willdeal
withthe completecubicsplineinterpolant,(1)above;errorestimates arefairlystraightforwardto makeand
makebestpossible,andthe naturalsplinecanthenbetreatedasaperturbedversionofthe completespline.
Periodic splines will largely be left to the reader.
3 642:573:01Fall2000
3. Notational Di(cid:11)erences Between these notes and Atkinson pp. 166{176.
Atkinson denotes the function value at x by y ; these notes use f .
i i i
Atkinson writes, but never uses, the form
s(x)=ai+bix+cix2+dix3 for xi−1 (cid:20)x(cid:20)xi ;:
We shall implicitly use, but rarely write, the form
p (x)=f +b (x−x )+c (x−x )2+d (x−x )3
i i i i i i i i
for the i-th polynomial written in Taylor form with base point equal to x , the left endpoint of the i-th cell.
i
Note that the equations
Atkinson j These notes
y =s(x ) = f
i i i
s0(x ) = b
i i
M =s00(x ) = 2!(cid:1)c
i i i
s(3)(x ) = 3!(cid:1)d
i i
connect the values and derivatives of the spline with the coe(cid:14)cients of p (x) when the latter is written in
i
Taylor form.
4. Cubic Spline Interpolants as the Solution of a Variational Problem. It is a remarkable
fact that the existence and uniqueness of cubic spline interpolants, as well as some of their approximation
properties, can be established in a nonconstructive way. fWhy do something nonconstructive in a course in
numericalanalysis? Becauseitillustratestheideaof(cid:12)ndingsolutionsbysolvingvariationalproblems,which
is an important idea in applied mathematics generally, and one which recurs in the (cid:12)nite-element method
for numerical solution of boundary value problems of di(cid:11)erential equations.g The basic lemma is sometimes
called Holladay’s identity.g
Lemma: Inthegeneralsetupofthesenotes,letf(x)betwicecontinuouslydi(cid:11)erentiable. Supposes(x)
is a cubic spline interpolant of f(x) that satis(cid:12)es either the natural-interpolant endpoint conditions
00 00
s (x )=0 and s (x )=0;
0 n
the complete-interpolant endpoint conditions
0 0 0 0
s(x )=f (x ) and s(x )=f (x );
0 0 n n
or the periodic endpoint conditions (assuming also that f(x )=f(x ))
0 n
0 0 00 00
s(x )=s(x ) and s (x )=s (x ):
0 n 0 n
If g(x) is any C2-function that interpolates f(x) at the nodes fx g (including f(x) itself), then
i
Z Z Z
xnjg00(x)j2dx= xnjs00(x)j2dx+ xnjs00−g00j2dx
x0 x0 x0
holdswhens(x)satis(cid:12)esthenatural-interpolantendpointconditions. Thesameistrueforthecompletespline
interpolants(x)ifg(x)alsosatis(cid:12)esthe\completeendpointconditions"g0(x )=f0(x )andg0(x )=f0(x ),
0 0 n n
and is true for the periodic spline interpolant if g0(x )=g0(x ).
0 n
4 642:573:01Fall2000
Non-Digression Before the Proof. Having spent some timZe in Hilbert space, we should recognize
xn
thatanythinginvolvingtheintegralofthesquareofafunction,like jg00(x)j2dx,reallyinvolvesan\inner
x0
product," in this case having the form
Z
hf;gi= xnf00(x)g00(x)dx:
x0
This is not (strictly speaking) an inner product, but rather a \quadratic form" de(cid:12)ned on pairs of functions
f; g|which for our present purposes(1) we can take as belonging to C2([a;b]). Unlike the familiar inner
product that does not involve derivaZtives,this inner product can give a zero norm to nonzero functions. All
one can deduce from knowing that xnjg00(x)j2dx=0 is that g00(x)(cid:17) 0, which does not imply that f (cid:17)0
but rather that f00(x) (cid:17) 0 and thusx0that f0(x) (cid:17) A and f(x) (cid:17) Ax+B for some constants A; B 2 R.
However, everything else we know about inner products will work pretty well with this one. In particular,
we can see what really makes \Holladay’s identity" go:
Sub-lemma: Let (cid:27)(x) 2 C2([a;b]) satisfy the \Dirichlet(2) node conditions" (cid:27)(x ) = 0 at all of the
i
nodesfx :0(cid:20)i(cid:20)ng. Lets(x)be anycubicsplinewithnodesatthefx g|thatis,afunctions2C2([a;b])
i i
whose restriction to each subinterval J =[x ;x ] is a cubic polynomial. Then hs;(cid:27)i=0 under any of the
i i i+1
following three additional hypotheses:
(1) s(x) satis(cid:12)es the \natural end conditions" s00(x ) = 0 = s00(x ) (and no end conditions, other than
0 n
taking the value zero, need be imposed on (cid:27)(x));
(2) No end conditions are imposed on s(x), but (cid:27)(x) satis(cid:12)es the \Neumann complete end conditions"
(cid:27)0(x )=0=(cid:27)0(x );
0 n
(3) s(x)satis(cid:12)esthe\periodicendcondition"s00(x )=s00(x )and(cid:27)(x)satis(cid:12)esthe\periodicendcondition"
0 n
(cid:27)0(x )=(cid:27)0(x ).
0 n Z
Proof. Consider the part of the integral hs;(cid:27)i= xns00(x)(cid:27)00(x)dx over the interval [x ;x ], and
i i+1
x0
integrate by parts twice, letting \u" always be the derivative of s and \dv" always be the derivative of (cid:27).
One obtains
#
Z xi+1 xi+1 Z xi+1
s00(x)(cid:27)00(x)dx=s00(x)(cid:27)0(x) − s(3)(x)(cid:27)0(x)dx
xi =s00(x)(cid:27)0(x)#xxii+1 −(xsi(3)(x)(cid:27)(x)#xi+1 −Z xi+1s(4)(x)(cid:27)(x)dx)
#xi # xi xi
xi+1 xi+1
=s00(x)(cid:27)0(x) −s(3)(x)(cid:27)(x) because s(4)(x)(cid:17)0:
xi xi
Because s00(x) is continuous on [a;b]|it \takes the same value from the left at each node x as it does from
i
theright"|ifweaddtheseequationsfori=0ton−1thecontributionsfromthe(cid:12)rsttermonthelastset-o(cid:11)
line will\telescope." The contributionsfromthe secondtermonthat line will be zero,because (cid:27)(x )=0 at
i
all nodes. Therefore
Z
b nX−1
hs;(cid:27)i= s(x)(cid:27)(x)dx = [s00(x )(cid:27)0(x )−s00(x )(cid:27)0(x )]
i+1 i+1 i i
a i=0
=s00(b)(cid:27)0(b)−s00(a)(cid:27)0(a):
(1)
Strictlyspeaking,thisisa\sesquilinearform"becauseitisconjugate-linear,ratherthanlinear,intheargumentg. However,we
shallonlydealwithreal-valuedfunctions,sothecomplex-conjugationwillreallyneverhaveane(cid:11)ect.
(2)
\Homogeneous boundary conditions"|requiring functions to take the value zero (possibly in some generalized sense) at the
boundaryofaregion|traditionallybearthenameDirichlet.
5 642:573:01Fall2000
In case (1), clearly hs;(cid:27)i=s00(b)(cid:27)0(b)−s00(a)(cid:27)0(a)=0−0=0, and similarly in case (2) where (cid:27)0(a)=0=
(cid:27)0(b). While it maynotbe the casethatboth terms ofs00(b)(cid:27)0(b)−s00(a)(cid:27)0(a) arezeroin the third case,the
terms will cancel if s00(a)=s00(b) and (cid:27)0(a)=(cid:27)0(b).
Proof of the lemma. Write g = s−(s−g), let (cid:27)(x) denote s(x)−g(x), and observe that if hs;(cid:27)i = 0
then we can write
Z
b
jg00(x)j2dx=hg;gi=hs;si−2hs;(cid:27)i+h(cid:27);(cid:27)i=hs;si+h(cid:27);(cid:27)i
a Z Z
b b
= js00(x)j2dx+ js00(x)−g00(x)j2dx
a a
which is the conclusion of the lemma. If s00(x) is a natural cubic spline interpolant of f(x), then since
(cid:27)(x ) = 0 at all the nodes, we are in case (1) of the sublemma and hs;(cid:27)i = 0. If s(x) is a complete cubic
i
spline interpolant of f(x), then (cid:27)(x) = s(x) −g(x) satis(cid:12)es (cid:27)0(x) = s0(x) −f0(x) = 0 for x = x and
0
x = x , so we are in case (2) of the sublemma and hs;(cid:27)i = 0. Finally, if s(x) satis(cid:12)es the end conditions
n
s0(x )=s0(x ) ands00(x )=s00(x )(as wellas the interpolationconditions(x )=f(x )=f(x )=s(x )),
0 n 0 n 0 0 n n
then (cid:27)0(x) = s0(x)−g0(x) takes equal values at x and x also, so we are in case (3) of the sublemma and
0 n
hs;(cid:27)i=0. Thus the conclusion of the lemma holds in all cases.
Z
b
Corollary: The natural cubic spline interpolantof f(x) (if it exists) minimizes jg00(x)j2dx over the
a
class of all C2 interpolants g(x) of f(x) that equal f(x) at all nodes. The complete cubic spline interpolant
of f(x) (if it exists) minimizes the same integralover the class of all C2 interpolants of f(x) that equal f(x)
at all nodes and satisfy the endpoint conditions g0(x )=f0(x ) and g0(x )=f0(x ). Same for the periodic
0 0 n n
interpolant. In all these cases, the spline is the unique minimizing function.(3)
Proof. In all these cases, the appropriate case of Holladay’s identity says
Z Z Z
b b b
jg00(x)j2dx= js00(x)j2dx+ js00−g00j2dx:
a a a
Z
b
Evidentlyther.h.s.is(cid:21) js00j2dx,andequalitywillholdifandonlyifthecontinuousfunctions00(x)−g00(x)
a
is identically zero on [a;b]. So the integral of jg00j2 must be strictly largerunless the C2 function s(x)−g(x)
satis(cid:12)es the di(cid:11)erential equation y00 = 0 on [a;b]. The solutions of that equation are the linear functions
y(x)=Ax+B. Ineither case,since thefunction s(x)−g(x)=Ax+B takesthe value0bothatx=x and
0
x=x we must have A=0=B and s(x)−g(x)=0 identically, i.e., g(x) is identical to s(x). We already
n
saw minimization, and that’s uniqueness.
Wecanusethisuniquenessresulttoproveexistence,althoughinanon-constructiveway. Theargument
goes as follows. For a given function f(x) there are 4n−2 equations in the coe(cid:14)cients of the polynomials
p (x) that must be satis(cid:12)ed in order for them, considered piecewise, to determine a cubic spline interpolant
i
of f(x):
f =f(x ); i=0; :::; n−1
i i
p (x )=f(x ); i=0; :::; n−1
i i+1 i+1
p0(x )−p0 (x )=0; i=0; :::; n−2
i i+1 i+1 i+1
p00(x )−p00 (x )=0; i=0; :::; n−2
i i+1 i+1 i+1
and (cid:12)nally either the periodic endpoint conditions or
00 00
p (x )=0 and p (x )=0 or
0 0 n−1 n
0 0 0 0
p (x )=f (x ) and p (x )=f (x ):
0 0 0 n−1 n n
R
(3) The integral bjg00(x)j2dx hasaquasi-physical interpretationas\the energyrequiredtobendalinearlyelasticstrip of metal
a
intothecurvethatinterpolatesthefunctionvalues." Thispropositionisconsequentlysometimesreferredtoasa\minimum-energy"
characterizationofcubicsplineinterpolants.
6 642:573:01Fall2000
Theinterpolationconditionsareinhomogeneousequationsontheff ; b ; c ; d g;alltheremainingconditions
i i i i
arehomogeneous(r.h.s.zero)exceptforthe complete-splineendpointconditions. Thetotalequationcount
is 4n, and there are 4n unknown coe(cid:14)cients. By a well-known theorem on systems of linear equations,
either every inhomogeneous case of these equations has a unique solution or else there exists a solution of
the system
f =0; i=0; :::; n−1
i
p (x )=0; i=0; :::; n−1
i i+1
p0(x )−p0 (x )=0; i=0; :::; n−2
i i+1 i+1 i+1
p00(x )−p00 (x )=0; i=0; :::; n−2
i i+1 i+1 i+1
and (cid:12)nally either
00 00
p (x )=0 and p (x )=0 or
0 0 n−1 n
0 0
p (x )=0 and p (x )=0
0 0 n−1 n
orthe periodicendpointconditionshold(depending onthecaseconsidered),inwhichnotallofthe numbers
ff ; b ; c ; d g equal zero. Now suppose the latter case could actually occur|that there were a solution
i i i i
of these homogeneous equations. The resulting function s(x) de(cid:12)ned by the p (x)’s would indeed be C2,
i
and it would be a complete, natural or periodic spline interpolant of the function f(x) which is identically
zero. But f(x) identically zero also satis(cid:12)es the de(cid:12)ning conditions to be such a spline interpolant of itself.
By uniqueness of the spline interpolant, s(x) = 0 identically. But then each p (x) = 0 identically, so all
i
its coe(cid:14)cients are zero, so all the ff ; b ; c ; d g are zero contrary to their choice! We conclude that the
i i i i
equationswhichmustbesolvedto(cid:12)ndthecoe(cid:14)cientsofthepolynomialsmakingupasplineinterpolantare
determinedby 4n equationsin 4nunknowns wherethe equationshaveno nontrivialhomogeneoussolutions;
the equations are therefore uniquely solvable for any choice of r. h. s.’s, and so they can always be solvedto
yield the coe(cid:14)cients of the polynomials in the cubic spline interpolant of any given f(x). This establishes
the existence of cubic spline interpolants (we already had uniqueness).
We shall give an algorithm for (cid:12)nding the coe(cid:14)cients later on. Meanwhile, it is nice to know that one
is talking about things that actually exist when one proves theorems about cubic spline interpolants.
5. Setting up the equations. Determining the fb g is actually most easily done by writing the
i
polynomials in Newton form and solving the following Hermite interpolation problem (the b are not yet
i
known and are to be determined)
0 0
p (x )=f ; p (x )=b ; p (x )=f ; p (x )=b :
i i i i i i i i+1 i+1 i i+1 i+1
Newton interpolation then gives us p (x) as
i
p (x)=f +p [x ;x ](x−x )+p [x ;x ;x ](x−x )2+p [x ;x ;x ;x ](x−x )2(x−x )
i i i i i i i i i i+1 i i i i i+1 i+1 i i+1
and the divided di(cid:11)erences are explicitly computed in terms ofthe values off , f , b , b and f[x ;x ]
i i+1 i i+1 i i+1
in the (cid:12)rst divided-di(cid:11)erence table of Table 1.
To get this into Taylor form is easy: rewrite the single factor (x−x ) in the form
i+1
x−x =x−x +x −x =(x−x )−h
i+1 i i i+1 i i
and the Newton form turns into the Taylor form
p (x)=f +p [x ;x ](x−x )
i i i i i i
+fp [x ;x ;x ]−h (cid:1)p [x ;x ;x ;x ]g(x−x )2
i i i i+1 i i i i i+1 i+1 i
+p [x ;x ;x ;x ](x−x )3 :
i i i i+1 i+1 i
These calculations(cid:12)nd p (x) \lookingrightfrom x ." Onthe other hand, \lookingleft fromx " almost
i i i
the same calculations would have done the interpolation problem to (cid:12)nd pi−1(x), because cubic Hermite
7 642:573:01Fall2000
Table 1
Looking right from x :
i
x f(x ) p [(cid:1);(cid:1)] p [(cid:1);(cid:1);(cid:1)] p [(cid:1);(cid:1);(cid:1);(cid:1)]
k i i i
f[x ;x ]−b b +b −2f[x ;x ]
x f b i i+1 i i i+1 i i+1
i i i h h2
i i
b −f[x ;x ]
i+1 i i+1
x f f[x ;x ]
i i i i+1
h
i
x f b
i+1 i+1 i+1
x f
i+1 i+1
Looking left from x :
i
x f(xk) pi−1[(cid:1);(cid:1)] pi−1[(cid:1);(cid:1);(cid:1)] pi−1[(cid:1);(cid:1);(cid:1);(cid:1)]
x f b f[xi−1;xi]−bi bi−1+bi−2f[xi−1;xi]
i i i −hi−1 h2i−1
xi fi f[xi−1;xi] bi−1−−fh[ix−i1−1;xi]
xi−1 fi−1 bi−1
xi−1 fi−1
interpolation imposes two conditions at each node and thus the problem of determining pi−1(x) at xi is
symmetric to that of determining pi(x) there. In Newton form, pi−1(x) must be
pi−1(x)=fi+pi−1[xi;xi](x−xi)
+pi−1[xi;xi;xi−1](x−xi)2
+pi−1[xi;xi;xi−1;xi−1](x−xi)2(x−xi−1):
Again the Newton form can be put into Taylor form with base-point x by rewriting
i
x−xi−1 =x−xi+xi−xi−1 =(x−xi)+hi−1
and so the Taylor form of pi−1(x) at xi is
pi−1(x)=fi+pi−1[xi;xi](x−xi)
+fpi−1[xi;xi;xi−1]+hi−1(cid:1)pi−1[xi;xi;xi−1;xi−1]g(x−xi)2
+pi−1[xi;xi;xi−1;xi−1](x−xi)3 :
Since these polynomials are in Taylor form with base point x , one can read their derivatives at x right
i i
o(cid:11) their Taylor coe(cid:14)cients. That the 0-th and 1st Taylor coe(cid:14)cients are equal simply expresses the facts
that both polynomials take the value f and that both have derivative b at that point. The coe(cid:14)cients of
i i
(x−x )2, however,arethe secondderivativesof the polynomialsup to a factor of2!. Hence a necessaryand
i
su(cid:14)cient condition for the polynomials to have equal second derivatives at the break point x is that
i
pi[xi;xi;xi+1]−hi(cid:1)pi[xi;xi;xi+1;xi+1]=pi−1[xi;xi;xi−1]+hi−1(cid:1)pi−1[xi;xi;xi−1;xi−1]
which Table 1 tells us is
f[xi;xhi+i1]−bi −hi(cid:1) bi+bi+1−h22i(cid:1)f[xi;xi+1] = f[xi;−xhi−i−11]−bi +hi−1(cid:1) bi−1+bi−h2i2−(cid:1)1f[xi;xi−1]
8 642:573:01Fall2000
in terms of known quantities. With both sides multiplied by hi(cid:1)hi−1 this is
hi−1(cid:1)ff[xi;xi+1]−big−hi−1(cid:1)fbi+bi+1−2(cid:1)f[xi;xi+1]g=hi(cid:1)fbi−f[xi;xi−1]g+hi(cid:1)fbi−1+bi−2(cid:1)f[xi;xi−1]g
and this, with the unknown b ’s kept on the r. h. s. and the known quantities moved to the l. h. s., becomes
i
an equation
3(cid:1)hi(cid:1)f[xi;xi−1]+hi−1(cid:1)f[xi;xi+1]=hi(cid:1)bi−1+(hi−1+hi)(cid:1)2bi+hi−1(cid:1)bi+1
that is to be satis(cid:12)ed at each of the interior nodes, i = 1; :::; n−1. One more piece of beauti(cid:12)cation and
we’ll have something that’s worth trying to solve: divide the i-th equation by (hi−1+hi) and put
(cid:22) = hi ; (cid:21) = hi−1 :
i i
hi−1+hi hi−1+hi
Theseare\convexcombinationcoe(cid:14)cients"|pairsofnonnegativenumberssatisfying(cid:22) +(cid:21) =1|andwith
i i
the unknown b ’s on the l. h. s. where unknowns belong, the equations determining the b ’s take the form
i i
(cid:22)ibi−1+2bi+(cid:21)ibi+1 =3(cid:1)f(cid:22)if[xi;xi−1]+(cid:21)if[xi;xi+1]g
for i=1; :::; n−1. We need two more equations, and we take those to be
2b =2b (=2f0(x ))
0 0 0
2b =2b (=2f0(x ))
n n n
whichis tosaythatwe shallbe (cid:12)nding the coe(cid:14)cientsinthe complete cubic spline interpolantoff(x).
The factors of 2 are present only for the sake of symmetry of notation.
The equations determining the b ’s then take the form
i
2b (cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1) = 2f0(x )
0 0
(cid:22) b +2b +(cid:21) b (cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1) = 3(cid:1)f(cid:22) f[x ;x ]+(cid:21) f[x ;x ]g
1 0 1 1 2 1 1 0 1 1 2
(cid:22) b +2b +(cid:21) b (cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1) = 3(cid:1)f(cid:22) f[x ;x ]+(cid:21) f[x ;x ]g
2 1 2 2 3 2 2 1 2 2 3
(cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1) ... (cid:1)(cid:1)(cid:1) (cid:1)(cid:1)(cid:1) ... (cid:1)(cid:1)(cid:1)
2b = 2f0(x )
n n
and are (n+1) linear equations in the (n+1) unknowns (b ; :::; b ).
0 n
6. The R. H. Sides of the equations we have just derived are fairly easy to compute and to estimate
in terms of f0(x) (when it exists); indeed, it should be obvious to the reader that each can only be about 3
times(4) as large as any value of f0(x). But we don’t yet know that the equations can be solved at all, or
that (b ; :::; b )T will not exhibit unpleasant behavior. It is therefore a marvelous piece of luck that these
0 n
equations belong to one of the most tractable classes of linear systems: those whose matrix of coe(cid:14)cients
is strictly diagonally dominant.(5) Since we don’t need the most general properties of these systems,
the following lemma is enough to give us what we need now. To accomodate ancient habits we’ll write the
system in the form AX =Y, where X and Y are vectors, for the duration of this x only. Let us agree that
the norm of a vector X =(x ; :::; x )T will be the ‘1-norm, i.e., will be given by
0 n
kXk1 =maxfjxij:i=0; :::; ng:
(4)
Infact,ifwefurther\normalized"theseequationsbydividingthetopandbottomonesby2andtheothersby3,wewouldsee
thatwhattheseequationsdoisbalanceweightedaveragesofthebi’sonthel.h.sidesagainstweightedaveragesofvalues(atnearby
points)off0(x)onther.h.sides.
(5)
Seetheseparatenotesonsystemsofthiskindforamorecompletediscussionofthesematricesthanisgivenhere.
9 642:573:01Fall2000
Lemma: Let AX = Y be an (n+1)(cid:2)(n+1) system of linear equations whose matrix of coe(cid:14)cients
has the \tri-diagonal" form
2 3
2 (cid:21) 0 0 (cid:1)(cid:1)(cid:1) 0 0
0
6 7
66(cid:22)1 2 (cid:21)1 0 (cid:1)(cid:1)(cid:1) 0 0 77
6 7
6 7
6 0 (cid:22) 2 (cid:21) (cid:1)(cid:1)(cid:1) 0 0 7
6 2 2 7
6 7
6 0 0 (cid:22) 2 (cid:1)(cid:1)(cid:1) 0 0 7
6 3 7
666 0 ... ... ... ... 0 0 777
6 7
6 7
4 0 0 0 0 (cid:1)(cid:1)(cid:1) 2 (cid:21)n−15
0 0 0 0 (cid:1)(cid:1)(cid:1) (cid:22) 2
n
where(6) theo(cid:11)-diagonalcoe(cid:14)cientssatisfyj(cid:22) j+j(cid:21) j(cid:20)1. ThenforanygivenpairofvectorsX; Y satisfying
i i
the equations, we have kXk1 (cid:20) jykj, where 0 (cid:20) k (cid:20) n is any index(7) at which the maximum de(cid:12)ning
kXk1 =maxfjxij:i=0; :::; ng is attained. In particular,
kXk1 (cid:20)kYk1 :
Proof. Let X andY be suchvectorsandAbe sucha matrix,andlet the index k be as described. Then
(the appropriate terms being missing in the cases k=0 and k =n) we have
(cid:22)kxk−1+2xk+(cid:21)kxk+1 =yk
2xk =yk−(cid:22)kxk−1−(cid:21)kxk+1 :
Taking absolute values and using the triangle inequality, we get
2jxkj(cid:20)jykj+j(cid:22)kjjxk−1j+j(cid:21)kjjxk+1j(cid:20)jykj+j(cid:22)kjjxkj+j(cid:21)kjjxkj
=jy j+(j(cid:22) j+j(cid:21) j)jx j(cid:20)jy j+jx j:
k k k k k k
Finally, subtracting jx j from both sides gives us
k
kXk1=jxkj(cid:20)jykj(cid:20)kYk1
as desired.
Corollary: AX = Y has a unique solution X for every given Y (and that solution must satisfy the
norm inequality kXk1 (cid:20)kYk1).
Proof. Any solution of AX = 0 satis(cid:12)es kXk1 (cid:20) k0k1, so X = 0. Thus the associated homogeneous
(n+1)(cid:2)(n+1) system AX = 0 has only the trivial (all-zero) solution. It is a standard theorem of linear
algebra that this property implies that AX =Y is solvable (and uniquely so) for each given Y.
fRemark: Thatcorollarywasprovedunconstructively,butthesematriceshaveacomputationalproperty
that enables one to show that the system AX = Y can always be solved by Gaussian elimination without
pivoting. And the Gaussian elimination is trivial: in each column there is only one sub-diagonal element.
Thus machine solution of such \tridiagonaldiagonally dominant" systems is very fast compared to systems
in general. See the separate notes onXdiagonally dominant matrices, in which it is shown that the \strict
diagonaldominance"conditionja j> ja jispreservedastheGaussianeliminationalgorithmprogresses.
ii ij
j6=i
Strict diagonal dominance implies invertibility, by the Ger(cid:20)sgorintheorem.g
(6)
This lemma is basically a very special case of the general Ger(cid:20)sgorin theorem on locating the eigenvalues of a matrix|and
estimatingthenormofitsinverse|intermsofthediagonalelementsofthematrix.
(7)
Thisisa\discretemaximumprinciple,"muchlikethemaximumprincipleofpotentialtheory.
10 642:573:01Fall2000