7
$\begingroup$

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 -

  1. 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}$.

  2. 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})$.

  3. 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})$.

  4. 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}$


$\endgroup$
6
$\begingroup$

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.

Interpolation using cubic splines

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.