I am writing a code snippet in Python to do an interpolation using cubic splines. I have first done the math, and then attempted to implement the pseudo code in Python. However, I think i might have messed up with the running index or a coefficient.
Would someone please be kind enough to check my math? The resulting curve is not smooth, does not fit the interior points and is all over the place.
Let $(x_{0},y_{0}),(x_{1},y_{1}),\ldots,(x_{n},y_{n})$ be $n+1$ points in $\mathbb{R}^2$. A spline is a piece-wise polynomial function of the form –
$$S(x)=\begin{cases}
S_{0}(x),& \text{if } x_{0}\le{x}\lt x_{1}\\
\vdots\\
S_{i}(x),& \text{if } x_{i}\le{x}\lt x_{i+1}\\
\vdots\\
S_{n-1}(x),& \text{if } x_{n-1}\le{x}\lt x_{n}\\
\end{cases}$$
$S_{i}(x)$ is a cubic polynomial with $4$ four coefficients, $\forall{i}$. There are $n$ intervals and a total of $4n$ unknowns. So, we need $4n$ conditions.
Suppose $S_{i}(x)$ has the form $S_{i}(x)=A_{i}(x-x_{i})^3+B_{i}(x-x_{i})^2+C_{i}(x-x_{i})+D_{i}$. The first and second derivatives of the cubic polynomials are:
$\begin{aligned}
S_{i}(x)&=A_{i}(x-x_{i})^3+B_{i}(x-x_{i})^2+C_{i}(x-x_{i})+D_{i}\\
S_{i}'(x)&=3A_{i}(x-x_{i})^2+2B_{i}(x-x_{i})+C_{i}\\
S_{i}''(x)&=6A_{i}(x-x_{i})+2B_{i}\\
\end{aligned}$
Also,
$\begin{aligned}
S_{i}(x_{i})&=D_{i}\\
S_{i}'(x_{i})&=C_{i}\\
S_{i}''(x_{i})&=2B_{i}\\
\end{aligned}$
We define $h_{i}=x_{i}-x_{i-1}$. We have:
$\begin{aligned}
S_{i-1}(x_{i})&=A_{i-1}h_{i}^3+B_{i-1}h_{i}^2+C_{i-1}h_{i}+D_{i-1}\\
S_{i-1}'(x_{i})&=3A_{i-1}h_{i}^2+2B_{i-1}h_{i}+C_{i-1}\\
S_{i-1}''(x_{i})&=6A_{i-1}h_{i}+2B_{i-1}\\
\end{aligned}$
Four properties of cubic splines
The spline should satisfy meet the below criteria –
-
The function $S(x)$ will interpolate all data points. $S(x)$ must be continuous. And so in each interval, $S_{i}(x_{i})=y_{i}$ and $S_{i-1}(x_{i})=y_{i}$.
-
The curve $S(x)$ should be smooth without jumps. $S'(x)$ must be continuous on the interval $[x_{i},x_{i+1}]$. Therefore, the slopes at each interior points must match. $S_{i}'(x_{i})=S_{i-1}'(x_{i})$.
-
The curve $S(x)$ should be not have any abrupt changes in its bentness or convexity. $S''(x)$ will be continuous on the interval $[x_{i},x_{i+1}]$. $S_{i}''(x_{i})=S_{i-1}''(x_{i})$.
-
A choice of one of the following two conditions at the end points $(x_{0},y_{0})$ and $(x_{n},y_{n})$
(a) The natural spline: $S'_{0}(x_{0})=0=S'_{n-1}(x_{n})$
(b) The clamped cubic spline : $S'_{0}(x_{0})=f'(x_{0})$ and $S'_{n-1}(x_{n})=f'(x_{n})$ where $f$ is presumably the function, we are trying to approximate.
Let's determine the $4n$ conditions.
Evaluation of the coefficients $A_i,B_{i},C_{i},D_{i}$
1) The first condition yields
$A_{i-1}h_{i}^3+B_{i-1}h_{i}^2+C_{i-1}h_{i}+D_{i-1}=D_{i}$
2) The second condition yields
$3A_{i-1}h_{i}^2+2B_{i-1}h_{i}+C_{i-1}=C_{i}$
3) The third condition yields
$6A_{i-1}h_{i}+2B_{i-1}=2B_{i}$
The above equations can be somwhat simplified, if we substitute $S_{i}''(x_{i})=2B_{i}=z_{i}$. Thus, we have $B_{i}=z_{i}/2$.
1) The last equation becomes :
$6A_{i-1}h_{i}=2(z_{i}/2)-2(z_{i-1}/2)=z_{i}-z_{i-1}$.
$A_{i-1}=\frac{z_{i}-z_{i-1}}{6h_{i}}$
2) The first equation becomes :
$\begin{aligned}
C_{i-1}h_{i}&=y_{i}-y_{i-1}-A_{i-1}h_{i}^3-B_{i-1}h_{i}^2\\
C_{i-1}&=\frac{y_{i}-y_{i-1}}{h_{i}}-(A_{i-1}h_{i}^2+B_{i-1}h_{i})\\
C_{i-1}&=\frac{y_{i}-y_{i-1}}{h_{i}}-\left(\frac{z_{i}-z_{i-1}}{6h_{i}}h_{i}^2+\frac{z_{i-1}}{2}h_{i}\right)\\
C_{i-1}&=\frac{y_{i}-y_{i-1}}{h_{i}}-h_{i}\left(\frac{z_{i}+2z_{i-1}}{6}\right)\\
\end{aligned}$
We define $b_{i}=\frac{y_{i}-y_{i-1}}{h_{i}}$.
In all of the above the equations, the running index $i$ goes from $1$ to $n$. Thus, we now have our equations for determining the coefficients.
$A_{i-1}=\frac{z_{i}-z_{i-1}}{6h_{i}}$
$B_{i-1}=\frac{z_{i-1}}{2}$
$C_{i-1}=b_{i}-h_{i}\left(\frac{z_{i}+2z_{i-1}}{6}\right)$
$D_{i}=y_{i}$
The system of equations in $z_{0},z_{1},z_{2},\ldots,z_{n-1}$
If we substitute these values in the second equation $3A_{i-1}h_{i}^2+2B_{i-1}h_{i}+C_{i-1}=C_{i}$, we should get a recurrence relation between $z_{i}$ –
$\begin{aligned}
3\frac{z_{i}-z_{i-1}}{6h_{i}}h_{i}^{2}+2\frac{z_{i-1}}{2}h_{i}+b_{i}-h_{i}\left(\frac{z_{i}+2z_{i-1}}{6}\right)&=b_{i+1}-h_{i+1}\left(\frac{z_{i+1}+2z_{i}}{6}\right)\\
\left(\frac{2z_{i}+z_{i-1}}{6}\right)h_{i}+\left(\frac{z_{i+1}+2z_{i}}{6}\right)h_{i+1}&=b_{i+1}-b_{i}\\
h_{i+1}z_{i+1}+2z_{i}(h_{i}+h_{i+1})+z_{i-1}h_{i}&=6(b_{i+1}-b_{i})
\end{aligned}$
for $i=1,2,3,\ldots,n-1$
This system of linear equations in $z_{0},z_{1},z_{2},\ldots,z_{n-1}$ can be represented in the matrix form as –
$\begin{bmatrix}
h_{1} & 2(h_{1}+h_{2}) & h_{2} & 0 & \ldots & 0 & 0 \\
0 & h_{2} & 2(h_{2}+h_{3}) & h_{3} & \ldots & 0 & 0 \\
0 & 0 & h_{3} & 2(h_{3}+h_{4}) & \ldots & \\
\vdots & & & & \ddots & \\
0 & 0 & 0 &\ldots & h_{n-1} & 2(h_{n-1}+h_{n}) & h_{n}
\end{bmatrix}
\begin{bmatrix}
z_{0}\\
z_{1}\\
z_{2}\\
\vdots\\
z_{n-1}
\end{bmatrix}
=\begin{bmatrix}
6(b_{2}-b_{1})\\
6(b_{3}-b_{2})\\
6(b_{4}-b_{3})\\
\vdots\\
6(b_{n}-b_{n-1})
\end{bmatrix}$
This matrix has $n-1$ rows and $n+1$ columns. So, we need two additional conditions.
For natural splines, $z_{0}=0=z_{n}$. The first column and the last column in the above system of linear equations can be eliminated, resulting in,
$\begin{bmatrix}
2(h_{1}+h_{2}) & h_{2} & 0 & \ldots & 0 & 0 \\
h_{2} & 2(h_{2}+h_{3}) & h_{3} & \ldots & 0 & 0 \\
0 & h_{3} & 2(h_{3}+h_{4}) & \ldots & \\
\vdots & & & \ddots & \\
0 & 0 & 0 &\ldots & h_{n-1} & 2(h_{n-1}+h_{n})
\end{bmatrix}
\begin{bmatrix}
z_{1}\\
z_{2}\\
\vdots\\
z_{n-1}
\end{bmatrix}
=\begin{bmatrix}
6(b_{2}-b_{1})\\
6(b_{3}-b_{2})\\
6(b_{4}-b_{3})\\
\vdots\\
6(b_{n}-b_{n-1})
\end{bmatrix}$
Best Answer
I have checked my math and python code with regular and boundary cases. After, a few corrections, the math and the code are working like a charm! The intuition, the math behind cubic splines and the python code snippet can be found in this Jupyter notebook.