SPLINE INTERPOLANT - John Travis

Consider the cubic polynomials Sk(x) = ak + bk(x-xk-1) + ck(x-xk-1)2+ dk(x-xk-1)3 , where xk-1< x < xk, for k=1,2,...,n
 

We want to piece these cubics together to get:

-- Interpolation...for k=1, 2, ...,n

(a) Sk(xk-1) = yk-1

(b) Sk(xk) = yk

-- Derivative Continuity...for k=1,2,...,n-1
(c) Sk'(xk) = Sk+1'(xk)
-- Second Derivative Continuity...for k=1,2,...,n-1
(d) Sk''(xk) = Sk+1''(xk)
Apply these conditions to the cubic spline formulas.
(a) yields (a') ak = yk-1, for k=1,2,..., n

(b) yields (b') ak + bkDxk+ ckDxk2 + dkDxk3= yk, for k=1,2,...,n

(c) yields (c') bk + 2ckDxk+ 3dkDxk2 = bk+1, for k=1, 2,...,n-1

(d) yields (d') 2ck + 6dkDxk= 2ck+1, for k=1, 2,...,n-1

Solving (d') yields dk= ( ck+1 - ck) / (3Dxk), for k=1, 2,...,n-1
 

If we define cn+1 to be S''(xn)/2, then we can consider the above as valid for k=1,2,...,n. Plugging this into (b') yields

yk-1+ bkDxk+ ckDxk2+ [( ck+1 - ck) / (3Dxk)]Dxk3= yk, or

yk-1+ bkDxk+ ( ck+1 + 2ck)Dxk2/3 = yk, for k=1,2,...,n, or by solving for bk,

bk= [ yk - yk-1- ( ck+1 + 2ck)Dxk2/3 ]/Dxk, or

(b'') bk= (yk - yk-1)/Dxk- ( ck+1 + 2ck)Dxk/3

Relabeling by replacing k with k+1 yields for k=0,1,...,n-1
(b''') bk+1= (yk+1 - yk)/Dxk+1- ( ck+2 + 2ck+1)Dxk+1/3
Similarly, plugging into (c') yields for k=1, 2,...,n-1
bk+ 2ckDxk+ 3[( ck+1 - ck) / (3Dxk)]Dxk2= bk+1, or

bk+ ( ck+1 + ck)Dxk= bk+1, for k=1,2,...,n-1.

So, using (b'') and (b''') here yields for k=1,2,...,n-1,
(yk- yk-1 )/Dxk- ( ck+1 + 2ck)Dxk/3 + ( ck+1 + ck)Dxk= (yk+1 - yk)/Dxk+1- ( ck+2 + 2ck+1)Dxk+1/3, or

( ck+2+ 2ck+1 )Dxk+1- ( ck+1 + 2ck)Dxk+ 3( ck+1 + ck)Dxk= 3(yk+1 - yk)/Dxk+1- 3(yk - yk-1)/Dxk, or

ckDxk+ 2ck+1(Dxk+1+ Dxk) + ck+2Dxk+1= 3(yk+1 - yk)/Dxk+1- 3(yk - yk-1)/Dxk.

Writing out a few of these yields (noting that cn+1= S''(xn)/2)

k=1        c1Dx1+ 2c2(Dx2+ Dx1) + c3Dx2= 3(y2 - y1)/Dx2- 3(y1 - y0)/Dx1.

k=2        c2Dx2+ 2c3(Dx3+ Dx2) + c4Dx3= 3(y3 - y2)/Dx3- 3(y2 - y1)/Dx2,

. . .

k=n-1     cn-1Dxn-1+ 2cn(Dxn+ Dxn-1) + cn+1Dxn= 3(yn - yn-1)/Dxn- 3(yn-1 - yn-2)/Dxn-1.
 

Collecting this in a matrix yields n-1 equations in the n+1 unknowns. Adding two additional equations (usually endpoint conditions) yields a system that can be solved for the unknown coefficients ck. Using the bold equations above, we can then go back and determine all of the other coefficients.