Errors in cubic spline interpolation

ARCHIEF

Journal of Engineering Mathematics, Vol. 3, No. 2, ApriI 1969

Wolters-Noordhoff Publishing Groningen

Printed in theNetherlands

Errors in Cubic Spline Interpolation

P. SONNEVELD

Deift University of Technology, Dept. of Mathematics, Delft, the Netherlands. (Received January 13, 1969)

SUMMARY

This paper deals with the problem of finding error bounds for cubic spline interpolation of functions of the class

C4 [a, b], and C [a, b], by examining a relationship between cubic spline interpolation and piecewise cubic Hermitian interpolation.

The method also gives an indication of what happens, in the case of almost uniform meshes, especially if the

approx-imated function is in the class C5[a, b].

Comparison is made with recent work carried out by K. E. Atkinson [3], in dealing with natural cubic spline inter-polation.

1. Introduction

Let h denote a subdivision of the interval [a, b] of the real x-axis

h={a=x0.czx1<...<x=b}.

(1.1)

A cubic spline, or spline on the mesh h is defined by the following relations yh(x)eC2[a, b]

(1.2)

for

x1 x x.1.

,

O i n-1.

The class of all splines on the mesh h will be denoted by Yh.

From the continuity requirements in (1.2) the following set of linear relations between the

polynomial coefficients can be derived a. b

q0

64Yo4Po2Pi

6'1Ynr4Pn2Pn-i

,

q-I n 2

q =

q-1 +-{p1p...1},

d. r

q+1q

0jn-1

i+ i c.

for ljn-1

where

1jn-1

(1.3)

Technische Hogeschool

₁₀₇

Delit

Llyi

h1

-The set (1.3) consists of 3n equations between 4n + 3 coefficients.

lt follows that if n + 3 independent extra relations are given, a spline is uniquely determined.

We can define several types of splines, interpolating a given functionf(x) in the n + i

mesh-points x, and satisfying two extra so-called end conditions.

Journal of Engineering Math., Vol. 3 (1969) 107-117

= j+ i

u= lAi.

j

j+1

108 P. Sonneveld

The first of these interpolating splines, denoted by Y (f; x), is defined by

y,,(f;x)eYh.

yh(f;x)=f(xl),

0in.

(1.4)

y(f;x)=f'(x),

¡=0 and

i=n.

All coefficients Pi can be determined by solving the system (1.3.a) together with the

condi-tions: p0=f'(x0), p=f'(x). The system (1.3.b,c) yields the values of the q, and finally (l.3.d) yields the values of the rj.

The spline interpolation y,,(f; x) will be referred to as the optimal interpolating spline, since

it has the best approximation property

Jba

y'(f; x)

f"(x)}2dx

J {z"(x) f"(x)}2dx

for all ze

_h

(1.5) Another type of spline interpolation, denoted by y (f; x), is defined by

y'(f; x)e Y,

y'(f; x.)

=f(x),

0 i n.

y"(f; x) =f"(x),

i = O and ¡

n.

This sphine interpolai ion an be computed by suh1ítutíng f" (x0), q =f" (xe) in (J .3.b), and

solving the Pj from the system (t .3.a,b); (I .3.c)and (I .3.d) serve to dctemniine q1, q,

_{..., q.}

and r0, r1, ..., r_ by substitution.

A third type of spline interpolation is the so-called natural sphine interpolation, defined by 9,,(f; x) Y,,

9,,(f;x,)=f(x1),

0ín.

_(1.7)

yh(f;xj)=O,

i=Oand i=n.

This spline has the minimum norm property

Jb

b

{j(f;

x)}2dxJ

{z"(f; x)}2dx (1.8)

for all interpolating z(f x)eC2[a, b]

Defining the splines 10(x) and 1(x) by

loeYh, mnEl'h,

10(x1) =1(x)=0,

0in.

l'(xo)

=

l(x)

= 1

l(x) = l''(x) =

0.

will show that the following relation between y,, (f; x) and j,, (f; x) holds

yh(f; x) = y'(f; x)f"(x0)I0(x)f"(x)1(x)

(1.10) The error Eh(f; x) of a spline interpolation is defined by

E,,(f; x)

=f(x)y,,(f;

x).

(1.11)

Now let

1h11= max h and IIg(x)tI =

max g(x)I,

a5x5b

then iffC4 [a, b], the error satisfies the following inequality

IE,,(f;x)II

K(f)

11h114

IE(f;

x)II K*(f)11h114. (1.12)

Journal of Engineering Math.,VoL 3 (1969) 107-117

(1.6)

The purpose of this paper is to find the best possible values for the

_{constants K(f) and}

K*(f).

In [2] and [3] an integral representation for the error is used in order to obtain error

estimates. Eh(f; x)

=JK(x t)fIv(t)dt

(1.13) E

_(f;

x) b K* (x, t)flv(t)dt

where K(x, t) andK*(x, _{t) are functions of the class C2 [a, b] with respect to x and t.}

Obviously the following relations hold Eh(f;x)I

JbIK(x,

t)Idt.I!fIvlI

(1.14)

E(f;

x)I

K(x, t)dt.!IfIvII.

Atkinson [3] derives the following bounds for the integral of IK*(x, t)t, incases of uniform meshes

.00156IhIt2(x_x)(x+i_x)<JIK*(x,

t)!dt<

<1.18HhII2(x-x1)(x11-x),

x.xxj1, 0in-1.

(1.15)

An improved version of this result will be given in section III.

2. Relation between CubicSpline Interpolation and Piecewise Cubic Hermitian Interpolation The classical piecewise cubic Hermitian interpolation can be defined by the following relations

y8(f; x)eC' [a,b]

yH(f;

x) = y+a1(xx1)+4b1(xx)2 +*c1(xx1)3

forxxx1, 0in-1

(2.1)

y'H(f; x) =f'(x1)

y(f; x1) =f(x)

_{0 i n.}

Therefore the piecewise cubic Hermitian interpolation is a piecewise cubic function, having

a continuous first derivative on [a, b]. Furthermore, iffe C4 [a, b], theerror can be expressed as follows

(xxj2(xx+

i)2

EH(f; x) =f(x) yH(f; x)

= 4! fIV()

for x1xx1, x1x1, Oin-1.

(2.2)

which is a classical result.

Therefore the order of the approximation is_O(11h114),which is in agreement with the errors

Eh_(f, x) and E (f; x) of the spline interpolations.

The above description shows that a similarity exists between cubic_{spline interpolation and}

piecewise cubic Hermitian interpolation.

Let us denote the difference of the two approximations by z(x) in order to examine this

similarity

z(x) _=YH(f;

_{x)yh(f; x) = Eh(f; x)EH(f; x)}

_(2.3)

and similarly:

110 P.Sonneveld

z(x) =

zQ1(x), i=O

= z'(x1)

z (x) =YH (f; x) - y (f; x)

=

E (f; x) - Eli (f; x).

(2.4)

Then it follows from (1.12) and (2.2)

I!zH = O(11h114), IIz*II = O(hII4), if

fEC' [a, b]

(2.5)

Obviously z(x) and z*(x) are piecewise cubics, satisfying the relations z(x)eC' [a, b], z*(x)ECI [a, b]

z'(x1)=

E(f; x)E,(f; x.)

=

E(f; x.),

for 0 i n,

z(x1) = 0, z*(x1)_=0,

for 0 i n.

(2.6)

z*'(x1)

_{= E"(f;}

x.)E(f;x)

= E'(f;

_x),

for 0 i n.

Now define a set of functions Q(x) as follows

xx1)

(x_x1+1)2

x1 x

xj+ ,

O i n-1,

h+1 Q.(x) =

_{Ix - x._}

xx1)

h.

₎

x1_1xx1, 1in.

(2.7)

xxi_i, xx1+1.

It follows that Q.(x) is a piecewise cubic function belonging to C' [a, b], and satisfying the

relations

Q1(xi)= O

Q(xj)=5ij={? ' Y

for0in, 0jn.

(2.8)

By differentiating (2.7) two times with respect to x, we get

Q'(x+)

= 2

i+ i

for 0in, 0jn-1

Q' (Xj -) = {2 _{+ }}

for 0in, 1j_n.

(2.9)

Now from (2.8) and (2.6) it follows that z(x) and z* (x) can be expressed in terms of the

func-tions Q(x)

z*(x) = _zr'Q1(x),

_zr=z*(xt).

_(2.10)

1=0

From (2.9) it is obvious that z (x) and z (x) have second derivatives which in general are

discontinuous in the mesh-points. Furthermore using the definitions (2.3) and (2.4) we find the following relations

5z"(x1)

= E'(f;

x1)-5E(f; x1)

= 5E',(f;

x.)

z*"(x1)= 5E4(f; x.)

for i i n-1

(2.11)

where 5g (x)=g (x + ) - g (x -).

Substituting (2.9) and (2.10) in (2.11), we find the following_{n - I relations between the z and}

5E(f; x.)

=

E(f;x1)

(2.12)

A1z1+2z'+p1z1

=

E(f;x1)

where 2 and a, have the meaning as defined in (1.3). Furthermore we find from (2.3)

z'(x0)= z'(x) =0 _(2.13)

and from 2.4)

= E"(f;

x0)E(f; x0)

= E(f;

x0)

(2 14)

= E"(f;

x) -

E(f;

x)

= -

E(f; x)

from which it follows by substitution of (2.9) and (2.10)

2z'+z'= +4h1E,(f;x0)

2z'+z.1

= --hE(f;

x).

The system (2.12.a) together with (2.13), yields _{a uniquely solvable system of equations for} thez if 5E,(J; xj are known, and so does the _{system (2.12.b) together with (2.15) for the zí'.}

In order to obtain an expression for 5E',,(f; x.) we examine Taylor's expansion forf(x) in the

neighbourhood of x= x.

f(x) =f(x1) + (xx1)f'(x1) +(xx)2f"(x1)

_{+*(xx1)3f"'(x)}

x(

t)

6 Iv(t)dt. _(2.16) =S.(x)+E1(x).

The piecewise cubic Hermitian interpolation off(x) can be written as

yH(f; x)=yH(Sl;x)+yH(El; x), _(2.17)

and since S1(x) is its own piecewise cubic Hermitian interpolation _{this can be reduced}

yH(f;x)=S.(x)+ y11(E; x).

(2.18)

Hence it follows forEH(f; x):

E11(f; x)

_{= SI(x)+El(x)SI(x)yH(EI; x)}

=E. (x) - YE, (E1; x)

₌

EH (E1; x).

Now since E.(x1)=E(x1)=0, the following relations can easily be verified

Yii(X)_Ei(Xi){}+

i+ i

+E(x1+1)Q11(x),

forx1xx11

_(2.20)

yH(Ei;x)=EI(xl_l){IQQ1H}+

+E(x1__1)Q1 (x), for x1_

x x

Differentiating (2.20) two times with respect_{to x and substituting x=x, we find by making}

use of (2.9)

Journal of EnqineeringMath.,Vol. 3 (1969) 107-117

(2.15)

112 _{P. Sonneveld}

with

R

=

vJQï(t)fI%'(t)dt,

Oj n

v=phj+1,

0jn-1

y,,= 2,,h,, =h,,

The systems (2.25.a) and (2.25.b) have the same structureas the defining. equations (1.3),

and can also be derived directly from the defining equations. This is done by substituting

p

=

y(f; x.)

= f(x) z

in (1.3.a). Applying a Taylor series expansion with integral remainder

to the resulting right-hand members, yields again the systems (2.25).

Although this procedure is much easier to deal with, the present derivation of the equations (2.25) yields some more information, and seems to be more consistent with the development of an error analysis. One of the consequences of the system (2.25) is the possibility of deriving the kernels K(x, t) and K*(x, t) as defined in (1.13) and all their properties from the system (2.25) together with (2.3) and (2.20). This, however, is beyond the purpose of this paper.

3. Error Bounds

For the purpose of obtaining bounds on z(x) and z*(x)_{we make use of the fact that a uniform} upper bound for the right-hand members of(2.25)is also a uniform upper bound for the solu-tion. For the sake of simplicity this will be shown for the_{system (2.25.a) only.}

Journal of Engineering Math., VoL 3 (1969) 107-117

Ï fr

y(E1;x1+) =

2

{3E(x1+1)E(x+1)h+1}

h.+ 2 (2.21) 2

y(E;x1) =

Substituting the integral expression for Ej(x), given by (2.16), we can reduce (2.21)

to the

following form

y(E1; x+)

= h2

(x_t)fIv(t)dt

y(E1; xe)

=J

X' h2j

(x1t)f"(t)dt.

xj,. _{(x1 - t)2} (x

_-

t)2

i+1 (2.22)

Now if we substitute these results in (2.19) and bear in mind that E'(x1) =0, we finally arrive

at the following result

xi +i

E(f;x+)=

Q1(t)f'(t)dt, 0in-1

E(f; x.)

=

J

Q.(t)f"'(t)dt,

I i n.

Xi - (2.23)

from which it follows

b

c5E(f; x.)

=J

Q(t)fR'(t)dx =JQ(t)fw(t)dx fori

i n-1

. (2.24)

With the aid of (2.23) and (2.24), we can now reduce the equations (2.12), (2.13) and (2.14) to the following form

J2jz_i+2z+/Ajz+j

=

R,

1jn1

tz'o=z=0.

0jn

(2.25)

The system reads

z'0 = z;, = 0.

1jn-1

_(3.1)

Let z = max

_{i jn-1}

IzI, then the following inequality holds:

j2zj = 1R1A1z_j zI IRI +

IzI

since 2+u,=l. It follows

zI

IRI

Now let M = max IR), then

i

jn-i

IzM,

1jn-1.

In the same way we have

z7tM*,

Ojn

where M* = max IR) _(3.5)

In (2.25) R is given by an integral, which can be rearranged in the following form

'1

hh

¡

(1t)2t{h2

fIV(x+ht)_h2fW(x_ht)}dt

j+1.'

R

2

h+h+1J0

for 1jn-1

R0 =

_{J(1 t)2tf"(x0+h1t)dt}

_(3.6)

From (3.9) we can construct a bound on z(x) by using the relation (2.7) and (2.10)

z(x)

(xx)(x1+jx)

_{{z(x1+ 1x) z+1 (xx1)}}

for x.xx,1,0in-1

from which it follows

Iz(x)I ₌

(xx)(x+x) 11h113

24

for x1xx1, 0in-1

Journal of Engineering Math.,Vol. 3(1969)107-117

(3.11)

R = j-h,

I

Now since (1t)2t

jo

rl

(1t)2tf"(xht)dt.

is a non-negative function it follows

<

h.hJ+i(h+h+l)11fIvf1

Ijn-1

IR) ROI

hIjf"lI.

o

-

1

,3 dV

"n = -:

n

j

(3.7)

Using the mesh norm, (3.7) can be simplified by IR) _{I1h1I3 IIf"lI}

0jn.

Applying (3.4) we get

(3.8)

114

Finally using the relations (2.3) and (2.2)

Eh(f;x) =

EH(f;x)+z(x)

(x_xi)2(xi1_x)2fIv()

EH(f;x)

₂₄

we arrive at the result

11h1131 (x - xt)(x,+

ix)

(x x)(xj+

x) +

11f WII IEh(f;x)l 24

i

_h11

for x,xx1, Oinl.

The same procedure applied toz*(x) leads to the same result for E(f; x):

IE(f; x)I

(xx)(xi x)

(xx)(x+ x)

Ij

_{} lIfIl.}

(3.14)

This last result is somewhat more specific than the result Atkinson obtained in [3] for

_{E (f; x).}

In fact his bounds for the integral in (1.15) can be improved upon in the following way: Let t be any positive number. Then for every fixed xE[x; x.+1]

_{a functionfexists with}

If

!VII =1, so that:

K*(x, t)ldtt<

E'(f; x)f

K*(x,_t)ldt

Ja Ja

Itfollows from (315) and (3.14):

11h113)

ÇIK*(x,t)Idt

<

(xx)(x+1x)

{(X_X)(X1+1_X) +

h1

-Ia 24

Furthermore for any y > i a function g(x) exists, satisfying the relations:

gIv(x)=(_1y{1

(x_xi)_(xi+i_x)}

j+1 h+1

with

Then Ig"I

=1, and it follows from (3.6): R.

=(-1)'lRJ.

I

hJhJ+l(h+hì+l)(1

) Rl = h+h1+1

(1e),

IRAI =

(1),

e ROI

-Introducing the mesh ratio fi as follows:

II I fi= max

1jn

h we can write

R1=(-1YIR)

m IRM

11h IP

M=

24 1h113 j1 12

24 )

(v+2)(v+3) 12 (v+2)(v+3) M

Now we can express

z'

as follows

052M

and applying the same reasoning as in the beginning of this chapter_{we get} Journal of Engineering Math., VoI. 3(1969)107-117

P. Sonneveld (3.13) (3.15) (3.16) (3.17)

ljn-1

(3.18) ¡ a (3.19) (3.20) (3.21)

I5)

=ômax(MIRI)--Mm.

(3.22) Hence the z' satisfy the relations:

7* _{= 1)jz7'}

-i

Ojn.

(3.23)

mIzM

Hence the error E (g; x) can be expressed as follows

Ç(xx1)(x+j x) m+ O.(Mmfl

(3.24)

24 +

h11

J

with 0 O 1. So E(g; x)I has the following lower bound:

Ç(xx)(x+x)

m )

+ IE'(g;x). _(3.25)

(xx)(x+1x)<

24

Since

IE' (g; x)I <I K (x, t)p dtjg"Ij₌

t

K (x, t) dt

Ja Ja

and using (3.20), we find if y tends to infinity:

(xx)(x+x)Ç

1h112 b

<t

lK(x, t)ldt.

(3.26) 24

1(xx)(x1+1x)

+ f33 _{J =}

Combining (3.16) and (3.26) we arrive at:

(xx)(x+1x) I uhu2

i

'

<

_IK*(x,t)tdt

24

(xx)(x1+1x) +

f33 J = Ja

(xx)(x11--x) Ç

- 24 1(x_x)(x1+1_x)+ßuuhIu2} (3.27) If the mesh is uniform, then ß=1, and we have exactly:

Jb

(xx)(x+1x)

IK*(x,t)Idt

=

24 {(xx)(x1+1x)+llh2}. (3.28)

4. Improvement of the Bounds on Eh(f; x) if fcC5{a, b]

The error of optimal cubic spline interpolation can be more specifically bounded iffe C5 [a, b].

In that case an integration by parts can be applied to the integrals in (3.6)

hh+

₁₁₍₁

_t)2t{h+ Pv

(x+hj+it)_hPv (x_ht)}dt

=

RJ_h+,

_.'

hh

_j+t

h2_j+i

- 2

_h+h+1

11(1t)3

+ _(4.1) 12 h3

(h1 _hi)hilfIv+h+ 1+

_{hh+ fV(fl}

-

24 with x._1=<

_x+1, I=<Jn-1

Introducing the local mesh norm

IhII = max{h1, h.1}

Journal of Engineering Math., Vol. 3 (1969) 107-117 (4.2)

116 _{P. Sonneveld} and the local mesh ratio:

(11h11 IIhlI

ß=max

L ' L ( Vii _t1i+1 then it follows 1. z.Iz. z. ._

(

1. 3 i+1 i' i+1 iJ ₁ ¡

li \

l/

Now if the mesh is uniform, ß. equals I and hence we may regard i

-

_i/ß

as a local measure for nonuniformity

If the local nonuniformity is bounded by:

t<e,

for 1in-1

then

E1(1E)

ift

t.(le1)t(1t) ift.

Define the monotonic nondecreasing function 9(t) as

Il,

(c)=t4t(1_e)

Oc

then it follows

Ihh+1(h+1-h1)I

The expression (4.1) for R. can now be bounded:

IR) 11h113 ilf'll +

llhilll

fyi!)

1f the mesh is nearly uniform, the first term of (4.10) will be very small; quantitatively we can write

R) = o(t_{11h11 3l1f"II) + O (11h114 II! VII) (4.11)}

If t=O(1lhI12)

IR) 11hlI 11f VII + O(II/2lI 11f IVII) (4.12)

The result (4.12) has a bearing on parametric spline interpolation of smooth curves, where only nearly uniform meshes can be constructed if uniform meshes are required.

Applying (3.4) we get:

lzI 11h11 311flvii + 1IhI1411fVII . (4.13)

and according to (2.6) it follows:

lE(f; x)I

9(t)11h113Ilf"Il

+!IhII'lIf "II

(4.14)

This result is an extension of the statement

jE(f; x)l = O(11h114) (4.15) for nearly uniform meshes.

Finally by applying (3.10) and (3.12) we arrive at the following bound for E(f; x):

Journal of Engineering Math., Vol. 3 (1969) 107-117

(4.3) (4.5) (4.6) (4.7) (4.8) (4.9) (4.10)

E(f;

x)I

(x x)(x

- x)IIhIP

IIf:I!

[i

+ p9(e)

+ IL0

}

XXX11

(4.16)

which bound is much more specific than that given by (3.14). For uniform meshes the expres-sion (4.16) reads:

Eh(f;x)!

(xxJ(x+1x)JIhI2

IIfviiJIhI}

xjxx.

(4.17)

Further refinements of the error analysis can be made by examining the equations (2.25)

thoroughly; this, however, is beyond the purpose of this paper.

Extensions of the theory are possible in cases of periodic splines, which occur in approxima-tion of periodic funcapproxima-tions and in parametric spline interpolaapproxima-tion of closed smooth curves.

REFERENCES

.1. H. Ahlberg, E. N. Nilson and J. L. Walsh; The Theoryof splinesand their applications,Academic Press, New

York and London, 1969.

G. Birkhoff and C. de Boor; Error bounds for spline interpolations,_{JournalofMath. and Mech.,13,5. pp.827-835,}

1964.

K. E. Atkinson; On the order of convergence of natural cubic spline interpolation, SIAM Journal_{of Num. An.,}

5,1, pp. 89-101, 1968.