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 = x0 < x1 < · · · < xn = b which divide the interval J = [x0 , xn ] into n cells, the i-th being Ji = [xi , xi+1 ] of length hi = xi+1 − xi > 0. Implicitly we are given a function f (x) defined on J. We seek n polynomials {pi (x) : i = 0, . . . , n − 1} which will define a function “near” f (x) by defining the “nearby” function as equal to pi (x) on the cell [xi , xi+1 ]. Thus there is no “single analytic definition” 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 pi (x) will have the 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 define 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 satisfies two interpolation conditions, and in addition to those there may be “end conditions” at x0 = a and xn = b and conditions that link the “pieces” at the nodes. When we need a name for the function defined on all of J by using the polynomial pieces, we call it s(x); the piece defined on Ji will be called pi (x). When we need a single estimate for the lengths of the cells, we use h = maxi {|hi | : 0 ≤ i < n}. We shall tacitly assume that the function f (x) has as much differentiability 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 differentiable at least, and sometimes “twice” will be replaced by “four-” or “fivetimes.” People who know what these things are will realize that in some situations it would suffice that the functions have distributional derivatives that belong to L2 (J) (buzzword = Sobol¨ev space). {Remark: The convention used above to number intervals Ji and increments hi numbers the intervals 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 coding purposes, but the reason for this feeling is probably that the early versions of Fortran were unwilling 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 fight the textbook on this matter.} 2. Lower-degree and Uncoupled Examples. Here are some simple cases of piecewise polynomial approximation. If all the pi (x)’s are linear, then just the requirement that pi (x) interpolate the values of f (x) at the endpoints of Ji uniquely determines the polynomial pi (x): there is only one linear interpolant between the two values at xi and xi+1 . The function s(x) is thus uniquely determined. An estimate of its deviation from f (x) can be made on an interval-by-interval basis: in Ji we can write f (x) = pi (x) +

f 00 (ξx ) (x − xi )(x − xi+1 ) . 2!

h2 The maximum value of |(x − xi )(x − xi+1 )| is i , attained at the midpoint of Ji . We thus have a global, 4 uniform error estimate M 2 h2 · |f (x) − s(x)| ≤ 2 4 00 where M2 = max{|f (x)| : x ∈ J}. The object we have constructed is called the piecewise-linear interpolant of f (x), for obvious reasons. The next case would be to take all the pi (x)’s quadratic, and a natural first hope would be that we might be able to match up the first derivative of s(x) with that of f (x) at the nodes. That hope does not 1

642:573:01 Fall 2000

long survive the following count. If all the pi (x)’s are quadratic, then there are 3n coefficients required in order completely to determine s(x). The requirement that pi (x) take the same values as f (x) at each node imposes 2 conditions on each pi (x); there are n intervals, so that requires the coefficients to satisfy 2n linear equations. If we simply want to have s0 (x) continuous at each of the interior nodes x1 , . . . , xn−1 , we impose an additional (n − 1) equations p0i (xi+1 ) = p0i+1 (xi+1 ),

i = 0, . . . , n − 2

that the coefficients must satisfy: we have imposed 2n + (n − 1) = 3n − 1 linear conditions on the coefficients, and we would expect that “there is only one degree of freedom remaining.” Certainly an attempt to make s0 (xi ) = f 0 (xi ) at each node would add (n + 1) equations, bringing the count to 4n, and it does not seem 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 [α, β] and some number γ, there is a unique quadratic polynomial q(x) that interpolates f (x) at the endpoints and has q 0 (α) = γ; in Newton form that polynomial is q(x) = f (α) + γ · (x − α) +

f [α, β] − γ · (x − α)2 β−α

as the reader may see by constructing a divided-difference table. One can now compute q 0 (β) = 2 · f [α, β] − γ directly. So if in piecewise-quadratic interpolation we specify a value for p00 (x0 ) = s0 (x0 ), then p00 (x1 ) is determined, and if s0 (x) is to be continuous at x1 then p01 (x1 ) is determined . . . and so on down the line; specifying the derivative at x0 has taken away the one degree of freedom left in the choice of s(x). {Take [α, β] successively to equal [x0 , x1 ], [x1 , x2 ], . . . .} The coefficients of the pi (x)’s are “coupled” at the nodes 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 order to determine the entire cubic spline uniquely move from left to right along the interval as one computes the successive pi (x)’s in Newton form (which incidentally happens 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 (x0 ) you can easily convince yourself that it may be a worse approximation! {Suppose f (x) is a constant function: if s0 (x0 ) 6= 0 then the piecewise-quadratic will oscillate around the constant, while the piecewise-linear would equal the constant. Of course that was a loaded example—but suppose f (x) wiggles a little in J0 but is constant in [x1 , xn ]. Then even matching s0 (x0 ) and f 0 (x0 ) may not help.} If we are to have s0 (x) behave nicely with respect to f 0 (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 pi (x) are to be cubics, then it will require 4n coefficients to determine all the polynomials. The requirement that the piecewise function they define be an interpolant of f (x) imposes 2(n − 1) + 2 = 2n conditions, all of which can be expressed linearly in the coefficients and which look as if they should be linearly independent. In piecewise-cubic Hermite interpolation we require that the first derivative of the i-th interpolating polynomial interpolate the derivative of f (x) at the nodes. This imposes 2n additional conditions and determines the pi (x)’s uniquely. Indeed, finding the pi (x)’s then becomes a purely “local” problem: one finds the cubic Hermite interpolant of f (x) on [xi , xi+1 ] and that is pi (x); no interaction between cells takes place. The whole situation is exactly analogous to piecewise-linear interpolation; all that happens is that the first derivative is interpolated along with the function. To see what happens on a single interval, look at Atkinson’s treatment of Hermite interpolation on a single interval [a, b], which is the example on pp. 162–163. 2

642:573:01 Fall 2000

Note that here we do have a good error estimate: applying Atkinson’s formula (3.6.15) piece-by-piece gives us M4 · h4 |f (x) − s(x)| ≤ 384 where M4 = max{|f (4) (x)| : x0 ≤ x ≤ xn }. The error → 0 as h → 0 twice as fast as the error for piecewise1 ≈ 10−2.6 in there to help make the error linear interpolation did, and it is nice to have the factor 384 smaller too. In general, however, one will have p00i−1 (xi ) 6= p00i (xi ), i = 1, . . . , n − 1; the first derivative of the interpolant is continuous, but there will usually be jumps in its second derivative (which is not well-defined at the nodes). The definition 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 p00i−1 (xi ) = p00i (xi )

i = 1, . . . , n − 1

so as to get a continuous second derivative. That already gives 2n + (n − 1) = 3n − 1 conditions on the coefficients, so there is not enough freedom remaining to force the first derivatives to agree with f 0 (x) at the nodes. What one can do is to force the first derivatives of the cubic pieces to agree with each other at the nodes: i = 1, . . . , n − 1 . p0i−1 (xi ) = p0i (xi ) This gives n − 1 additional conditions, bringing the count up to (3n − 1) + (n − 1) = 4n − 2. While the common values {bi : i = 1, . . . , n − 1} of the first derivatives will probably not satisfy bi = f 0 (xi ), we may hope that the difference will be fairly small if f (x) is sufficiently smooth and if the {hi }n−1 i=0 are sufficiently small. One expects that 2 additional conditions are needed if all the pieces of the piecewise cubic interpolant are to be determined uniquely. While there is no canonical choice of these conditions, three popular choices are (1) The Complete Cubic Spline Interpolant, produced by choosing p00 (x0 ) = f 0 (x0 )

and p0n−1 (xn ) = f 0 (xn ) ;

“the first derivative is correct at the endpoints.” (2) The Natural Cubic Spline Interpolant, produced by choosing the {bi } so that p000 (x0 ) = 0 and p00n−1 (xn ) = 0 ; “no curvature is present at the endpoints.” It can be argued that this is not at all natural, since f 00 (x) need not be zero at the endpoints, and indeed the error estimates for this interpolant are not as small as those for the complete spline. (While it is not obvious that this condition can be produced by a suitable choice of the {bi }, we shall see later that it can be.) (3) The Periodic Cubic Spline Interpolant, produced by choosing the {bi } so that p00 (x0 ) = p0n (xn ) and p000 (x0 ) = p00n (xn ) ; it is most natural to use these conditions with a function for which f (x0 ) = f (xn ), so that the function f (x) “wants to have a periodic extension” to the real line with period cells that are the translates of J. {Of course a function can have periodic derivatives without itself being periodic; the interested reader may want to contemplate the situation further.} Atkinson (pp. 166–176) deals with the natural cubic spline interpolant, (2) above. These notes will deal with the complete cubic spline interpolant, (1) above; error estimates are fairly straightforward to make and make best possible, and the natural spline can then be treated as a perturbed version of the complete spline. Periodic splines will largely be left to the reader. 3

connect the values and derivatives of the spline with the coefficients of pi (x) when the latter is written in 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. {Why do something nonconstructive in a course in numerical analysis? Because it illustrates the idea of finding solutions by solving variational problems, which is an important idea in applied mathematics generally, and one which recurs in the finite-element method for numerical solution of boundary value problems of differential equations.} The basic lemma is sometimes called Holladay’s identity.} Lemma: In the general setup of these notes, let f (x) be twice continuously differentiable. Suppose s(x) is a cubic spline interpolant of f (x) that satisfies either the natural-interpolant endpoint conditions s00 (x0 ) = 0 and s00 (xn ) = 0 , the complete-interpolant endpoint conditions s0 (x0 ) = f 0 (x0 ) and s0 (xn ) = f 0 (xn ) , or the periodic endpoint conditions (assuming also that f (x0 ) = f (xn )) s0 (x0 ) = s0 (xn ) and s00 (x0 ) = s00 (xn ) . If g(x) is any C2 -function that interpolates f (x) at the nodes {xi } (including f (x) itself), then Z

Non-Digression Before the Proof. Having spent some time Hilbert space, we should recognize Z xin n |g 00 (x)|2 dx, really involves an “inner that anything involving the integral of the square of a function, like x0

product,” in this case having the form

Z

xn

hf, gi =

f 00 (x)g 00 (x) dx .

x0

This is not (strictly speaking) an inner product, but rather a “quadratic form” defined 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 derivatives, this inner product can give a zero norm to nonzero functions. All Z xn

Because s (x) is continuous on [a, b]—it “takes the same value from the left at each node xi as it does from the right”—if we add these equations for i = 0 to n−1 the contributions from the first term on the last set-off line will “telescope.” The contributions from the second term on that line will be zero, because σ(xi ) = 0 at all nodes. Therefore Z b n−1 X s(x)σ(x) dx = [s00 (xi+1 )σ 0 (xi+1 ) − s00 (xi )σ 0 (xi )] hs, σi = a

i=0 00

= s (b)σ 0 (b) − s00 (a)σ 0 (a) . (1)

Strictly speaking, this is a “sesquilinear form” because it is conjugate-linear, rather than linear, in the argument g. However, we

shall only deal with real-valued functions, so the complex-conjugation will really never have an effect. (2) “Homogeneous boundary conditions”—requiring functions to take the value zero (possibly in some generalized sense) at the boundary of a region—traditionally bear the name Dirichlet.

or the periodic endpoint conditions hold (depending on the case considered), in which not all of the numbers {fi , bi , ci , di } equal zero. Now suppose the latter case could actually occur—that there were a solution of these homogeneous equations. The resulting function s(x) defined by the pi (x)’s would indeed be C2 , 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 satisfies the defining conditions to be such a spline interpolant of itself. By uniqueness of the spline interpolant, s(x) = 0 identically. But then each pi (x) = 0 identically, so all its coefficients are zero, so all the {fi , bi , ci , di } are zero contrary to their choice! We conclude that the equations which must be solved to find the coefficients of the polynomials making up a spline interpolant are determined by 4n equations in 4n unknowns where the equations have no nontrivial homogeneous solutions; the equations are therefore uniquely solvable for any choice of r. h. s.’s, and so they can always be solved to yield the coefficients 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 finding the coefficients 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 {bi } is actually most easily done by writing the polynomials in Newton form and solving the following Hermite interpolation problem (the bi are not yet known and are to be determined) pi (xi ) = fi , p0i (xi ) = bi , pi (xi+1 ) = fi+1 , p0i (xi+1 ) = bi+1 . Newton interpolation then gives us pi (x) as pi (x) = fi + pi [xi , xi ](x − xi ) + pi [xi , xi , xi+1 ](x − xi )2 + pi [xi , xi , xi+1 , xi+1 ](x − xi )2 (x − xi+1 ) and the divided differences are explicitly computed in terms of the values of fi , fi+1 , bi , bi+1 and f [xi , xi+1 ] in the first divided-difference table of Table 1. To get this into Taylor form is easy: rewrite the single factor (x − xi+1 ) in the form x − xi+1 = x − xi + xi − xi+1 = (x − xi ) − hi and the Newton form turns into the Taylor form pi (x) = fi + pi [xi , xi ](x − xi ) + {pi [xi , xi , xi+1 ] − hi · pi [xi , xi , xi+1 , xi+1 ]}(x − xi )2 + pi [xi , xi , xi+1 , xi+1 ](x − xi )3 . These calculations find pi (x) “looking right from xi .” On the other hand, “looking left from xi ” almost the same calculations would have done the interpolation problem to find pi−1 (x), because cubic Hermite 7

which is to say that we shall be finding the coefficients in the complete cubic spline interpolant of f (x). The factors of 2 are present only for the sake of symmetry of notation. The equations determining the bi ’s then take the form 2b0 µ1 b0 + 2b1 + λ1 b2 µ2 b1 + 2b2 + λ2 b3 ···

···

···

..

.

···

···

=

2f 0 (x0 )

···

···

=

3 · {µ1 f [x1 , x0 ] + λ1 f [x1 , x2 ]}

···

···

3 · {µ2 f [x2 , x1 ] + λ2 f [x2 , x3 ]}

···

···

= .. .

2bn

=

2f 0 (xn )

···

and are (n + 1) linear equations in the (n + 1) unknowns (b0 , . . . , bn ). 6. The R. H. Sides of the equations we have just derived are fairly easy to compute and to estimate in terms of f 0 (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 f 0 (x). But we don’t yet know that the equations can be solved at all, or that (b0 , . . . , bn )T will not exhibit unpleasant behavior. It is therefore a marvelous piece of luck that these equations belong to one of the most tractable classes of linear systems: those whose matrix of coefficients 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 § only. Let us agree that the norm of a vector X = (x0 , . . . , xn )T will be the `∞ -norm, i.e., will be given by kXk∞ = max{|xi | : i = 0, . . . , n} . (4)

In fact, if we further “normalized” these equations by dividing the top and bottom ones by 2 and the others by 3, we would see

that what these equations do is balance weighted averages of the bi ’s on the l. h. sides against weighted averages of values (at nearby points) of f 0 (x) on the r. h. sides. (5) See the separate notes on systems of this kind for a more complete discussion of these matrices than is given here.

This lemma is basically a very special case of the general Gerˇ sgorin theorem on locating the eigenvalues of a matrix—and

estimating the norm of its inverse—in terms of the diagonal elements of the matrix. (7) This is a “discrete maximum principle,” much like the maximum principle of potential theory.

10

642:573:01 Fall 2000

7. Approximation by an Interpolating Spline. Even when the interpoland f (x) is only C2 , it turns out that the spline interpolants approximate f (x) and their derivatives approximate f 0 (x) quite well, uniformly on J. For smoother interpolands, the approximations are good out to f (3) (x). For the C2 case, the estimate employs the “minimum-energy” characterization of the cubic spline interpolants, and the estimates are not of the order of magnitude in h that one will learn to expect for smoother functions. The basic tool is the Schwarz inequality for integrals, which asserts that for two (Lebesguemeasurable) real-valued functions f and g on an interval [α, β], the inequality s s Z β Z β Z β 2 |f (x)g(x)| dx ≤ f (x) dx · g(x)2 dx α

In the last line, we evaluated the integral of 1 over [ξ, x] and replaced the integral of the square-of-thedifference of 2nd derivatives over [ξ, x] by the integral over the whole [x0 , xn ], which could only have made the integral bigger. But now we can estimate both factors in the last line: it is obvious that |x − ξ|1/2 ≤ h1/2 , while if we replace g 00 by f 00 in Holladay’s identity on pp. 4–5 above we see that the nonnegativity of the integral of s002 then implies Z b Z b |f 00 (t)|2 dt ≥ |s00 (t) − f 00 (t)|2 dt . a

It would suffice for this that f (x) belong to the Sobol¨ ev space W 2,2 (J), which would imply that it had a continuous first derivative

but only that it had a second derivative in the “weak” or “distributional” sense that belonged to L2 (J). The fact that one can assume only this rather shadowy existence for the second derivative is the reason for the low order of approximation.

11

642:573:01 Fall 2000

Proof. Let x ∈ J be given, and let [xi , xi+1 ] be a cell to which it belongs. Then since f (t) − s(t) = 0 for t = xi and t = xi+1 , we can compute f (x) − s(x) {up to its sign} by integrating f 0 (t) − s0 (t) from the nearer endpoint of [xi , xi+1 ]—call it x∗ —to x: Z x 0 0 (f (t) − s (t)) dt ≤ max{|x − x∗ |} · max{|f 0 (t) − s0 (t)| : t ∈ J} |f (x) − s(x)| = x∗ s s Z b Z b 1 1 |f 00 (t)|2 dt = h3/2 · |f 00 (t)|2 dt ≤ h · h1/2 · 2 2 a a as advertised. These are “modern” estimates—estimates of this type will recur at the end of the second second-semester course, when one develops the finite-element method of solving boundary-value problems in differential equations. However, they’re not very sharp: the exponents in h3/2 and h1/2 really “should” be 4 and 3 respectively, at least when f (x) is a sufficiently smooth function. Developing estimates of that order requires a great deal of classical technique, and we shall demonstrate that technique in the next §. Meanwhile, even the coarse estimates available from the “energy integral” show that as h → 0 we can expect cubic spline interpolants of C2 functions to converge to the interpoland—and, even better, the derivatives converge also {although “one power of h more slowly”}. That fact at least reassures us that the cubic spline interpolant is a reasonable subject of investigation. 8. Classical Error Estimates for Cubic Splines. We can now start to estimate the difference between the given function f (x) and its complete cubic spline interpolant in the classical way, with bounds involving the maximum absolute values of the higher derivatives of f (x). This is time-consuming but each step is elementary. We shall carry out the steps in the following order: (A) Working locally (for x in a single cell [xi , xi+1 ]) we can rewrite the difference between f (x) and pi (x) in the form f (x) − pi (x) = f (x) − H(x) + H(x) − pi (x) where H(x) is the cubic Hermite interpolant of f (x) (and f 0 (x)) on [xi , xi+1 ], use the triangle inequality |f (x) − pi (x)| ≤ |f (x) − H(x)| + |H(x) − pi (x)| , and estimate the two terms separately. The first term is the error term from Hermite interpolation, for which bounds are known (Atkinson p. 162, (3.6.15)). (B) That leaves the second error term ei (x) = pi (x) − H(x) which is a polynomial of degree ≤ 3 satisfying ei (xi ) = 0 = ei (xi+1 ), e0i (xi ) = bi − f 0 (xi ), e0i (xi+1 ) = bi+1 − f 0 (xi+1 ) . If we set

then ei can easily be estimated in terms of max{|i |, |i+1 |}. It is clear from these definitions that 0 = 0, and the only consistent definition of n is to set it = 0 also. (C) The (n + 1)-tuple (0, 1 , . . . , n−1 , 0)T satisfies a system of (n + 1) linear equations of the form AX = Y with the same matrix of coefficients A as the system that determined (b0 , b1 , . . . , bn )T . The r. h. s. of that system can be written out explicitly and its elements estimated. Because AX = Y implies kXk∞ ≤ kY k∞ , the estimate thus obtained also controls max{|0 |, . . . , |n |}; using the result of part (B) above, we then make that estimate control the error cubics ei (x) (i = 0, . . . , n − 1), thus finishing the estimation process for the complete cubic spline interpolant. 12

642:573:01 Fall 2000

(D) The error in the natural cubic spline interpolant, and the errors in the cubic splines produced by other possible choices of endpoint conditions, are most easily handled by estimating the errors between those interpolants and the complete interpolant. Hypotheses will be introduced naturally along the way: usually they will require that f (x) have 4 or 5 continuous derivatives on J. For future reference we recall the standard notation Mk = max{|f (k) (x)| : x0 ≤ x ≤ xn } for all functions f (x) for which this expression makes sense. Recall also that h = max{h0 , . . . , hn } . 9. Estimate (A). Not much to do here: Atkinson’s (3.6.13) gives |f (x) − H(x)| ≤

h4 · M4 384

which is sharp (take f (x) = x4 ). 10. Estimate (B). The error cubic ei (x) is a polynomial whose degree is at most 3, and it has a root at each end of [xi , xi+1 ], so it must have the form ei (x) = (x − xi )(x − xi+1 )[a(x − xi ) + b] . Its derivative thus has the form e0i (x) = (x − xi+1 )[a(x − xi ) + b] + (x − xi )[a(x − xi ) + b] + a(x − xi )(x − xi+1 ) , and evaluation, first at xi and then at xi+1 , shows that the linear function in the square brackets has to interpolate the values i+1 i at x = xi and at x = xi+1 . −hi hi Therefore this linear factor is

i+1 (x − xi ) − i (xi+1 − x) , and so h2i

ei (x) = (x − xi )(x − xi+1 )

i+1 (x − xi ) − i (xi+1 − x) . h2i

The product |(x−xi )(x−xi+1 )| assumes its maximum at the midpoint of the cell [xi , xi+1 ], and the maximum value is (h/2)2 . The linear factor assumes its maximum value at an endpoint, so its maximum value can be max{|i |, |i+1 |} . This gives the estimate estimated by hi |ei (x)| ≤

since the derivative of the spline matches the derivative of f (x) at the endpoints. Because of what we know about the matrix of coefficients of this system of equations, we know that if we find a bound for the absolute values of all the r. h. sides of these equations, the same bound will dominate the absolute values of all the i ’s. The first thing on the agenda, therefore, is to find a sharp bound for those r. h. sides. One method of attack is to rewrite the differences on the r. h. sides as integrals. This method recurs again and again in error analysis, and so it will pay us to write it up in as general a way as possible. Thus let F (t) be a 5-times-continuously differentiable function defined on an open interval −η < t < η centered on 0, let |h| < η, and consider 3 · F [0, h] − F 0 (h) − 2 · F 0 (0) =

1 h

Z

h

{[F 0 (t) − F 0 (h)] + 2 · [F 0 (t) − F 0 (0)]} dt .

(11.2)

0

Integrate the term [F 0 (t) − F 0 (h)] dt by parts, differentiating the first factor and using as the integral of dt the function t, which takes the value zero at the endpoint 0. Similarly integrate the term 2[F 0 (t) − F 0 (0)] dt by parts, differentiating the first factor but this time using as the integral of dt the function (t − h) which zeros out at the endpoint h. The result will be that the integral is rewritten as −1 h

It follows (let g(s) = s(hi−1 + s)2 in the first integral, for example) that each of the the two integrals in the h4 hi−1 h4i hi · M4 respectively. So · i−1 · M4 and · curly brackets can be estimated in absolute value by hi−1 12 hi 12 the entire r. h. side is now estimated in absolute value by the expression hi−1 hi 1 · · (h2i−1 + h2i ) · M4 . 24 hi−1 + hi

and that estimates the error of the complete cubic spline interpolant, as desired. In order to make a worst-case estimate of the integrals above by the mean-value theorem, we had to act as if it were possible for f (4) (x) to equal +1 on one side of xi and −1 on the other side. While in fact f (4) (x) cannot do this if it is continuous, it is possible to find (draw pictures!) functions whose behavior is arbitrarily close to that, so no improvement in the estimate of the {|i |} is possible. The function x4 · sgn(x) has that impossible 4-th derivative, and its spline interpolant (Fig. 1) demonstrates that the slope of the interpolant can differ a good deal from that of the curve at a node. In the plot shown in Fig. 1 we have x0 = −1, x1 = 0, x2 = 1. This function is a limiting case of 4-times continuously differentiable functions, so they can be “arbitrarily close to this bad.” 12. Refinements. It can be shown that the constant 5/384 of the preceding § is the best possible constant valid for arbitrary 4-times differentiable f (x) and arbitrary nodes whose maximum spacing is h > 0. However, some improvements are possible in the case of equal increments (all hi = h). In this case, and assuming that the function F (x) is 5-times continuously differentiable, we can replace F (4) (s) in the formula (11.7) above by Z s

Evaluating the integral of the first term ds, interchanging the order of integration of the iterated integral, and computing the inner integral in the new iterated integral will all combine to give ( ) Z h h4 (4) (h − t)3 −1 −h2 00 F (0) + F (0) + (3t + h) F (5) (t) dt . (12.3) h 2 12 12 0 (h − t)3 (3t + h) ; this time it is nonnegative on the 12 interval of t between 0 and h, regardless of the sign of h.

Here again one should note the sign of the “kernel”

The relation (12.3) can now be used in the calculations in (11.8) above, in place of the relation (11.7) that we previously used. It is useful in two special cases. In the special case where hi = hi−1 = h of (11.8) above, it is easy to check that not only do the two terms containing f 00 (xi ) cancel, but also the two terms containing f (4) (xi ). As a result, the r. h. s. of the i-th equation for the {i } takes the form (Z ) Z h 0 (h + t)3 (h − t)3 −1 (5) (5) · (h − 3t) F (xi − t) dt + (h + 3t) F (xi + t) dt . (12.4) 4h 12 12 −h 0 h5 on each of [−h, 0] and [0, h], so the entire expression (12.4) can be estimated 30 5 h4 1 2h · · M5 = · M5 . This looks as though we have gained a whole “power-of-h” in absolute value by 4h 30 60 on our former error estimate (11.12)! That is too good to be true: all we have improved is our estimate of h4 h M5 , and although we get to multiply by in estimating max{|ei (x)| : x0 ≤ x ≤ max{|i | : 0 ≤ i ≤ n} ≤ 60 4 xn }, we find that, as in (11.12), we still have to write The integral of the “kernel” is

|f (x) − pi (x)| ≤ |f (x) − H(x)| + |ei (x)| ≤ 16

h4 h5 · M4 + · M5 ; 384 240

(12.5) 642:573:01 Fall 2000

the first term is still of order h4 . Nonetheless, the second term still becomes small in comparison to the first, h4 · M4 · (1 + O(h))—which is so for small values of h > 0 we can manage to get our error bounded by 384 almost 5 times as optimistic as the best-generally-possible error estimate. But we shouldn’t get our hopes up too high—even in the case of f (x) = x4 , whose 5-th derivative is identically zero and for which ei (x) vanishes identically on equally spaced nodes. The Hermite-interpolant part of the error will still produce the deviations shown in Fig. 2 {where f (x) = 1 − x4 , x0 = −1, x1 = 0, and x2 = 1} to remind us that the piecewise cubic will only resemble a curve that is not a cubic in case the pieces are fairly short. The “bulging” present even though the 1st and 2nd derivatives agree at the origin shows what to expect. On the pessimistic side, we can see that for very disparate choices of hi and hi+1 we have to expect considerable inaccuracy. Indeed, for 5-times-continuously differentiable f (x) the expression (11.8) is now −h4i−1 hi hi−1 h4i −1 (4) 5 + · f (xi ) + O(h ) . · · (12.6) 2(hi−1 + hi ) hi−1 12 hi 12 The absolute value of the coefficient of f (4) (xi ) in the leading term of this expression is 1 hi hi−1 · (h2i − h2i−1 ) hi hi−1 |hi − hi−1 | = . 24 hi−1 + hi 24

(12.7)

To estimate this, do calculus: maximize xy(y − x) on the triangle 0 ≤ x, 0 ≤ y ≤ h, y ≥ x. The gradient of this function is (y 2 − 2xy, 2xy − x2 ) which has no zeros interior to the triangle, so the maximum is attained on the boundary—obviously not on the line x = 0 (!)—so therefore on the line y = h. Maximizing hx(h − x) gives x = h/2, so the worst case of our estimate occurs for hi = h/2, hi−1 = h (or vice versa), and is h3 /96. It follows that at an xi which divides the interval [hi−1 , hi+1 ] in the ratio 2:1, one must expect an i of this size. Consequently, if the nodes {xi : 0 ≤ i ≤ n} divide the interval J in ratios alternately 1:2 and 2:1, one cannot expect a bound on the error cubics ei (x) which is much smaller than the numbers h4 h h3 · · |F (4) (xi )| = · |F (4) (xi )|, and thus an error estimate on the order of 4 96 384 |f (x) − pi (x)| ≤

2h4 · M4 384

(12.8)

would be about the best one could really plan on. Fig. 3 shows what can happen {with f (x) = 1 − x4 , x0 = −1, x1 = −1/2, x2 = 1/2, x3 = 1}. The bulge is quite unexpectedly large, and worst at a point where f (x) is very flat, so smoothness of the function f (x) doesn’t help at all. 13. Cubic spline interpolants of f (x) have the further property that their derivatives manage to approximate the derivatives of f (x), at least up to order 3. The work of proving this splits into two parts as before: estimate the error for the appropriate derivative of the Hermite interpolant on [xi , xi+1 ], and then estimate the size of the appropriate derivative of ei (x). For the first derivative of the Hermite interpolant on [xi , xi+1 ], we have the exact error term f [xi , xi , xi+1 , xi+1 , x](x − xi )2 (x − xi+1 )2

Because we are seeking a bound valid for all f (x) and the geometry is symmetrical, we can restrict our xi + xi+1 , i.e., the case in which x lies in the left half of [xi , xi+1 ]. There are two attention to the case x ≤ 2 2xi + xi+1 or not, and it is routine calculus to verify that the largest value cases according to whether x ≤ 3 of the estimate occurs at the midpoint of [xi , xi+1 ], giving the value hi hi M4 3 M4 h2i · · + ·h . = 4! 4 2 2 96

(13.6)

{Any interested reader seeking to verify this bound should move the origin of coordinates to the center of [xi , xi+1 ] before doing the calculus.} This is not the sharpest-possible estimate, but the constants are adequate and the order of magnitude certainly correct. For the first derivative of ei (x) on [xi , xi+1 ], we have h2i ·

one can apply the same method as above, viz., get a bound of the correct order for the error of the k-th derivative of the Hermite interpolant H(x) on [xi , xi+1 ] and get a bound of the correct order for the k-th derivative of ei (x) (the second and third derivatives of ei (x) will involve 1/h and 1/h2 respectively). The 18

642:573:01 Fall 2000

details are similar to those above, and are left as an exercise for any interested reader. I do not know what the best constants are, but K3 is probably irrational. It should perhaps be pointed out that since the {i } are exactly the differences {bi − fi0 } between the derivative of the spline and the derivative of the function at the i-th node, the error p0i (xi ) − f 0 (xi ) in the approximation of the derivative at the nodes is on the order of h4 when f (x) is 5-times continuously differentiable and the nodes are equally spaced. {This follows from our investigation of that case above.} However, at points other than the nodes the order of approximation is the same as that of the derivative of the piecewise Hermite interpolant, which is only h3 . Still, this fact is occasionally useful in applications where a good approximate derivative at some point is desired: one just makes sure that that point is a node (and so perhaps one has to make one more measurement of an empirical f (x)). 14. Error in the Natural Spline Interpolant. This is the spline interpolant with zero curvature at the ends defined on p. 3 above. All we need do is look at the Taylor form of pi (x) at xi (bottom of p. 7 above) and refer to Table 1 for the divided differences; we see that for i = 0 we have b0 + b1 − 2f [x0 , x1 ] 1 00 p (x0 ) = {f [x0 , x1 ] − b0 } − (14.1) 2! 0 h0 so if we are to have p000 (x0 ) = c0 (no harm in adding a little generality), that requirement replaces the 0-th equation in the system on p. 9 above by 2b0 + b1 = 3f [x0 , x1 ] − c0

The resulting matrix is just as strictly diagonally dominant as the one for the complete cubic spline, so it has the same solvability and “norm-reducing” properties that the previous one did. The natural cubic spline interpolant thus exists and is uniquely determined by f (x) (as we already knew from the minimum-energy approach to splines). There are some problems with the error estimates, however. We must now define 0 and n by adjoining the equations 2(f00 + 0 ) + (f10 + 1 ) = 3f [x0 , x1 ] − c0

for some sample point ξn ∈ [xn−1 , xn ]. However, it is now evident that unless the forced second derivatives c0 and cn at the endpoints are chosen equal to the actual values of the second derivatives at the endpoints, we’ll lose a factor of h2 in our uniform estimate of the functions |ei (x)|, i = 0, . . . , n: the r. h. s. of the system of equations AX = Y (X = vector of i ’s—Y = r. h. sides of the equations above and the corresponding equations above) that we must solve is only bounded by a multiple of h instead of a multiple of h3 , and therefore we can expect a uniform error estimate that is only on the order of h2 (recall the factor of h/4 that comes from the estimate in §10 above of |ei (x)| in terms of max{|i | : i = 0, . . . , n}). The pessimistic appearance of this estimate is deceptive: the only reason that it looks so bad is that it is required to hold uniformly for all points x ∈ [x0 , xn ]. But it is unreasonable to expect good approximation near a endpoint at which one has wilfully forced a derivative of the spline to take a value that f (x) thinks is wrong! Indeed, if c0 and cn are chosen equal to f 00 (x0 ) and f 00 (xn ) in the estimate above, then the O(h)terms drop out and the approximation becomes O(h4 ) as before. So in fact, the disturbance in the accuracy of the approximation turns out to be an end effect : the difference between a spline interpolant that does funny things at an endpoint and the complete cubic spline interpolant decays more rapidly than any power of h {or of 1/(the distance from x to the endpoint)} as h → 0 or as x moves away from the endpoint, and so the order-of-magnitude estimates for the errors in other interpolants become indistinguishable from the ones we have developed for the complete cubic spline interpolant, once one gets away from the offending endpoint. 15. End Conditions other than the Complete End Conditions. The plan for investigating the approximation error of cubic spline interpolants with end conditions different from those of the complete spline parallels that used for the complete spline. We begin by setting up some notation. Let {pi (x) : 0 ≤ i < n} and {bi : 0 ≤ i ≤ n} be the polynomials and slopes that define the complete cubic spline interpolant of the 4- or 5-times continuously differentiable function f (x) on [x0 , xn ]. Let another cubic spline interpolant be defined by the cubic polynomials {qi (x) : 0 ≤ i < n} on the respective intervals [xi , xi+1 ], with slopes {mi : 0 ≤ i ≤ n} at the nodes. Then their difference is a piecewise cubic function defined by {zi (x) : 0 ≤ i < n}, with qi (x) = pi (x) + zi (x),

This format is certainly inclusive enough to include the case of the natural cubic spline interpolant in the preceding §. Then it is again clear that the matrix of coefficients is “strongly diagonally dominant” Pas in the discussion above, so the equations will have a unique solution. {Indeed, it would suffice to have |αi | < 2 P and |βi | < 2 to insure the solvability of the equations, and in fact all we need in what follows is that the equations be [uniquely] solvable for all choices of the numbers k0 and kn .} 20

642:573:01 Fall 2000

In this case, however, we shall not control the {ζi } and the |qi (x)| by controlling the size of all the {ζi } uniformly. Rather, we shall control their size away from the endpoints by appealing to the equations µi ζi−1 + 2ζi + λi ζi+1 = 0,

i = 1, ..., n − 1 .

(15.5)

For the sake of simplicity and because this is the most important case, consider the case where all hi = h, so all µi = λi = 1/2. Then the equations all have the form 1ζi−1 + 4ζi + 1ζi+1 = 0,

i = 1, ..., n − 1 .

(15.6)

This is(9) a 2-nd order linear homogeneous difference equation (with constant coefficients) for the sequence {ζi : 1 ≤ i ≤ n − 1}. These are the discrete analogues of 2-nd order differential equations with constant coefficients. All such equations have solutions of the form ζi = (const.) · ri

(It is handy that r1 = r2−1 , and we shall use that fact without further reference). It is elementary to prove that any sequence {ζi } constructed from r1 and r2 by taking arbitrary constants A and B and setting ζi = A · r2−i + B · r2i

(15.10)

satisfies the difference equation, and that every solution is of that form. {For readers who want to check that assertion: determine A and B so that the values at i = 1 and i = 2 come out as the given sequence does, and by recursion all the values must come out the same.} It is easy to check that the equations governing the ζi ’s are of the form 2ζ0

+ .. .

α1 ζ1 1 ζi−1 2

.. . β0 ζ0

+

β1 ζ1

··· .. .

+ .. .

+ 2ζi

+

.. .

.. .

.. .

+

···

+ .. .

···

αn ζn

1 ζi+1 2

+ βn−1 ζn−1

+2ζn

= .. .

K0 .. .

=

0

.. .

.. .

(15.11)

= Kn

where the numbers K0 and Kn are two constants, different in general from the previous k0 and kn . Consider the system of equations with the same coefficient matrix .. .

··· .. .

.. .

+

2ηi

+

.. .

.. .

.. .

2η0 .. . 1 ηi−1 2 .. .

1 ηi+1 2

2ηn (9)

21

= .. .

2C0 .. .

=

0

.. .

.. .

=

2Cn

(15.12)

For a thorough treatment of linear difference equations, please see the separate set of notes on this subject.

642:573:01 Fall 2000

and solve these equations in the case where the r. h. s. is the “1-st standard basis vector” (1, 0, ..., 0)T (so C0 = 1, Cn = 0). For these equations we must have ηn = 0, η0 = 1, and for 0 ≤ i ≤ n we must have ηi = A · r2−i + B · r2i

(15.13)

for some constants A and B. Since these constants must satisfy 1 = η0 = A + B 0 = ηn = A · r2−n + B · r2n ,

(15.14)

it is elementary to solve for the coefficients and find that ηi = The estimate |ηi | ≤ √

1 · [r2i − r22n−i ] . 1 − r22n

(15.15)

√ 1 · 2 · |r2 |i = ( 3 + 1) · |r2 |i 3−1

(15.16)

is just as elementary, since 2n − i ≥ i holds for all 0 ≤ i ≤ n and |r2 | < 1. This estimate shows that the |ηi | decay exponentially as i increases, and—of crucial importance—the estimate of their size is independent of n. By symmetry or by repeating the argument, the solution of those equations with the r. h. s. (0, 0, ..., 0, 1)T is subject to a similar estimate √ (15.17) |ηi | ≤ ( 3 + 1) · |r2 |n−i . It follows {the reader should be sure to understand the details!} that the solutions of the equations as given are subject to the estimate √ |ηi | ≤ ( 3 + 1) · {|C0 | · |r2 |i + |Cn | · |r2 |n−i } ,

(15.18)

the r. h. s. of which is fairly small in the middle of the range of indices. Since the equations that originally determined the {ζi } can be rewritten in the form = .. .

2ζ0 .. . 1 ζi−1 2 .. .

.. .

.. .

.. .

+

2ζi

+

.. .

.. .

.. .

1 ζi+1 2

2ζn

K0 −

P i

=

0

.. .

.. . P

= Kn −

αi ζi

.. . (15.19)

i

βi ζi

we see that since it is known that max{|ζi |} ≤ max{K0 , Kn }, the r. h. sides of the 0-th and n-th equations are bounded in absolute value by 2 · max{K0 , Kn }, and so finally we see that the |ζi | are subject to the estimate √ (15.20) |ζi | ≤ ( 3 + 1) · max{K0 , Kn } · [|r2 |i + |r2 |n−i ] . This setting is general enough to cover the case of the natural cubic spline interpolant: in this case the i = 0 equations are 2b0 = 2f 0 (x0 ) for the complete spline 2m0 + m1 = 3f [x0 , x1 ] for the natural spline 22

where 1 is the familiar “error in the derivative of the complete cubic spline interpolant at the node x1 ”. We know from the results of §13 above that the absolute value of the r .h. s. of this equation is actually of the form h (15.23) − f 00 (x0 ) + O(h3 ) 2 so it’s certainly bounded. A similar equation and estimate are available for n at the n-th node. Consequently the errors {ζi : 0 ≤ i ≤ n} are subject to an estimate |ζi | ≤ K · h · [|r2 |i + |r2 |n−i ]

(15.24)

where K > 0 is some constant that we could estimate (in terms of M2 and M4 , say) if it were important to do so. The factor [|r2 |i + |r2 |n−i ] is “small in the middle of the range 1 ≤ i ≤ n − 1” and we can exploit this to make the differences {zi (x)} between the splines small, in two ways. First, suppose we want error estimates that are valid uniformly in an interval J(δ) = [x0 + δ, xn − δ], in which we “agree to stay at least a fixed distance δ > 0 from the endpoints.” Then in order for the cell [xi , xi+1 ] to be contained in J(δ)—as a moment’s geometrical thought will show—we must have δ/h ≤ i and (δ/h) + 1 ≤ n − i. This gives [|r2 |i + |r2 |n−i ] ≤ (const.) · e−kδ/h

whenever the index i belongs to an interval [xi , xi+1 ] in which x can lie. As a result, using the estimate developed in §10 above for the cubic on [xi , xi+1 ] with endpoint values 0, 0 and endpoint derivatives i , i+1 , but replacing the ’s by ζ’s, we see that we have an estimate |zi (x)| ≤ (Const.) ·

h2 −kδ/h ·e 4

(15.27)

h2 —look at e−kδ/h , which → 0 faster than any 4 power of h as h → 0! {Use L’Hˆopital’s rule on ln(ha e−kδ/h ) for a ≥ 0—explicit bounds are also easy to compute.} Thus the difference between the natural cubic spline interpolant of f (x) and its complete cubic spline interpolant “damps out like e−1/h ” uniformly in any fixed interval J(δ) contained in the open interval (x0 , xn ). The reader can easily see that this result can be applied to any of the cubic spline interpolants of more general type for which we set up the estimation machinery above.

valid uniformly for x ∈ J(δ). Never mind the

Another way to look at the estimate we just derived is the following: take the estimate |zi (x)| ≤ (Const.) ·

h2 −kδ/h ·e , 4

(15.28)

let a small error tolerance η > 0 be given, and try to make the r. h. s. of the estimate be < h4 η by taking δ > 2h sufficiently large. {Why h4 ? because that is the order of the error in the complete cubic spline, and since we are estimating error as measured from that interpolant, there is no hope of improving the order of 23

Because the expression in square brackets → ∞ as h → 0, we can make the r. h. s. of this inequality be > 2h by taking h sufficiently small; and because − ln(4η/(Const.)) is large and positive for small η > 0, once h is “sufficiently small” for such an it is also ”sufficiently small” for all smaller η’s. On the other hand, since h| ln h| → 0 as h → 0, we can make the whole r. h. s. small by taking h > 0 sufficiently small. Thus there is a function of h—the r. h. s. of this inequality—which is of the order O(h ln h), such that for sufficiently small h, if x ∈ [x0 , xn ] stays as far from x0 and xn as specified by that inequality, then the estimate M4 + η · h4 |f (x) − s(x)| ≤ 384

(15.31)

holds, where η > 0 is any preassigned positive number. This result is valid for a large class of cubic spline interpolants including the natural one.(10) Estimates of this type can be extended to the derivatives of cubic spline interpolants with little trouble by imitating the methods used for the natural spline above. The details can easily be filled in by the interested reader. Attached figures with numbers > 3 show how the effects of the most wrong-headed choices of b0 and bn damp out quite quickly when the hi are chosen equal and are only moderately small.

(10)

This observation was made in the first ediction of Atkinson’s book—cf. the inequality (3.84), p. 147—but was dropped from the