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$

1 Answer 1

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$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .