Cubic spline algorithm

Abstract: Here is a summary of end conditions for cubic spline interpolation, as described in Matlab Spline Toolbox. The reference cited in that toolbox is A Practical Guide to Splines, (Applied Math. Sciences Vol. 27, Springer Verlag, New York(1978),xxiv + 392p).

Introduction

A Cubic spline is piecewise cubic polynomial { Pn }0£ n <N which can be fitted to a sequence { (xn,yn) }0£ n <N of data points. With knots at xn it exactly interpolates the data points. The piecewise portions are defined so that at the knots the polynomial function and its first two derivatives are continuous. Additionnal end conditions are required to well define the cubic spline. Here are the equations corresponting to the previous assumptions :

Pn(xn) = yn
with 0£ n < N,
Pn(xn+1) = Pn+1(xn+1)
P'n(xn+1) = P'n+1(xn+1)
P''n(xn+1) = P''n+1(xn+1)
ü
ý
þ
with 0£ n < N-1.

1   Matlab Cubic Spline Formulation

1.1   Equations

There are several ways to define the local polynomial function. The following formulation has the advantage to use only three paramters an, bn and cn. With xn £ x < xn+1, we define
    Pn(x) = yn + an(x-xn) + bn(x-xn)2 + cn(x-xn)3     (1)
    P'n(x) = an + 2bn(x-xn) + 3cn(x-xn)2     (2)
    P''n(x) = 2bn + 6cn(x-xn)     (3)
    P'''n(x) = 6cn     (4)

1.2   Discretisation

First of all, we define the step hn for all 0£ n < N-1 : hn = xn+1 - xn

The idea in Matlab discretisation is to solve a linear system in {y'n} and define a relation between paramters (an, bn, cn) and y'n. Using the previous equations (1-4), we write the equation for polynomila function and it first derivative at knots xn and xn+1.

    P'n(xn) : an = y'n     (5)
    P'n(xn+1) : y'n + 2bnhn + 3cnhn2 = y'n+1     (6)
    Pn(xn) : yn = yn     (7)
    Pn(xn+1) : yn + y'nhn + bnhn2 + cnhn3 = yn+1     (8)

Substracting equation (7) from equation (8), we obtain that
y'n + bn hn + cn hn2 =
yn+1-yn
hn
and a first formula for bn
bn =
yn+1-yn
hn
- y'n
hn
- cn hn.     (9)

Then substituing equation (9) in equation (6)
2 ì
í
î
yn+1-yn
hn
- y'n - cn hn2 ü
ý
þ		 + 3 cn hn2 = y'n+1 - y'n

which leads to
cn =
1
hn2
ì
í
î		 y'n+1 + y'n - 2
yn+1-yn
hn
ü
ý
þ		 .     (10)

Finally, we define the equation satisfied by y'n using the continuity of the second derivative
P''n(xn) = P''n-1(xn).     (11)

We expand the left hand side
P''n(xn) = 2 bn     (12)
=
2
hn
ì
í
î
yn+1-yn
hn
- y'n - ì
í
î		 y'n+1 + y'n - 2
yn+1-yn
hn
ü
ý
þ		 ü
ý
þ
    (13)
=
2
hn
ì
í
î		 3
yn+1-yn
hn
- y'n+1 - 2 y'n ü
ý
þ
    (14)

and the right hand side
P''n-1(xn) = 2 bn-1 + 6 cn-1 hn-1     (15)
=
2
hn-1
ì
í
î		 -3
yn-yn-1
hn-1
+ 2 y'n + y'n-1 ü
ý
þ
    (16)

Using the continuity equation (11) we finally obtain the core equation of the linear system for 0 < n < N-1 :
hn y'n-1 + 2 (hn+hn-1)y'n + hn-1y'n+1 = 3
yn+1-yn
hn
+ 3
yn-yn-1
hn-1
    (17)

To fully define the linear system, we have to add end contions. In the next subsection we present the common end conditions.

1.3   End Conditions

Several conditions are possible :

The last step is the discretization of the end conditions.

This document was translated from LATEX by HEVEA.