Cubic Spline¶

We consider here a cubic spline passing through the points $$(x_i,y_i)$$ with $$a=x_1<\ldots<x_n=b$$, that is, a class function $${\mathcal C}^2$$ on $$[a, b]$$ and each restriction at the interval $$[x_{i-1},x_i]$$, $$1\leq i\leq n$$, is a polynomial of degree less than 3. We will note $$S$$ such a spline. His equation is given by

$S_i(x) = Ay_i + By_{i+1} + Cy''_i+ D y''_{i+1}, \qquad x_{i}\leq x\leq x_{i+1},$

where

$A = \frac{x_{i+1}-x}{x_{i+1} - x_i} \qquad \text{et} \qquad B = \frac{x-x_i}{x_{i+1} - x_i},$
$C = \frac{1}{6}\left(A^3-A\right)\left(x_{i+1}-x_i\right)^2 \qquad \text{et} \qquad D = \frac{1}{6}\left(B^3-B\right)\left(x_{i+1}-x_i\right)^2.$

If we derive this equation twice with respect to $$x$$, we get

$\frac{d^2S(x)}{d x} = Ay''_i + By''_{i+1}.$

Since $$A = 1$$ in $$x_i$$ and $$A = 0$$ in $$x_ {i + 1}$$ and conversely for $$B$$, we can see that the second derivative is continuous at the interface of the two intervals $$[x_{i-1}, x_{i}]$$ and $$[x_{i}, x_{i + 1}]$$.

It remains to determine the expression of $$y''_i$$. To do this, we will calculate the first derivative and impose that it is continuous at the interface of two intervals. The first derivative is given by

$\frac{dy}{dx}=\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}-\frac{3A^2-1}{6}(x_{i+1}-x_{i})y''_i+\frac{3B^2-1}{6}(x_{i+1}-x_{i})y''_{i+1}.$

We therefore want the value of the first derivative in $$x = x_i$$ over the interval $$[x_{i-1}, x_{i}]$$ to be equal to the value of the first derivative in $$x = x_i$$ over the interval $$[x_{i}, x_{i + 1}]$$; which gives us for $$i = 2, \dots, n-1$$

$a_iy''_{i-1}+b_iy''_i+c_iy''_{i+1}=d_i,$

with

$\begin{split}\begin{array}{l} a_i = \frac{x_i-x_{i-1}}{x_{i+1}-x_{i-1}}\\ b_i = 2\\ c_i = \frac{x_{i+1}-x_{i}}{x_{i+1}-x_{i-1}}\\ d_i = \frac{6}{x_{i+1}-x_{i-1}}\left(\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}-\frac{y_{i}-y_{i-1}}{x_{i}-x_{i-1}}\right). \end{array}\end{split}$

So we have $$n-2$$ linear equations to calculate the $$n$$ unknowns $$y''_i$$ for $$i = 1, \dots, n$$. So we have to make a choice for the first and last values ​​and we will take them equal to zero. We can recognize the resolution of a system with a tridiagonal matrix. It is then easy to solve it by using the algorithm of Thomas which one recalls the principle

$\begin{split}c'_i=\left\{ \begin{array}{lr} \frac{ci}{b_i}&i=1\\ \frac{c_i}{b_i-a_ic'_{i-1}}&i=2,\dots,n. \end{array} \right.\end{split}$
$\begin{split}d'_i=\left\{ \begin{array}{lr} \frac{di}{b_i}&i=1\\ \frac{d_i-a_id'_{i-1}}{b_i-a_ic'_{i-1}}&i=2,\dots,n. \end{array} \right.\end{split}$

The solution is then obtained by the formula

$\begin{split}\begin{array}{l} y''_n = d'_n \\ y''_i = d'_i-c'_iy''_{i+1} \qquad \text{pour} \qquad i=n-1,\dots,1. \end{array}\end{split}$