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},\]


\[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


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\)



\[\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}\]