# Not-a-knot cubic spline interpolation using tridiagonal solver

I am trying to write my own cubic spline interpolator. Given the formula for the cubic spline $$S_n(x) = a_n+b_n(x-x_n)+c_n(x-x_n)^2+d_n(x-x_n)^3$$ my interpolator works perfectly for the natural boundary conditions in which $$S_0''(x_0)=0=2c_0 + 6d_0(x-x_0) \to c_0=1$$ $$S_{n-1}''(x) = 0$$

where $$n = 0, 1, 2, ... n-2, n-1.$$

However, I cannot figure out how to formulate the not-a-knot boundary conditions in such a way as to be compatible with my tridiagonal matrix. I know that the not-a-knot boundary conditions dictate that $$S_0'''(x) = S_1'''(x)$$ $$S_{n-2}'''(x) = S_{n-1}'''(x)$$

My first instinct was that I could simply set, say for the first equation, $$d_0 = d_1$$, and plug this back into $$S_n(x)$$ for $$x_0 \leq x \leq x_2$$. Of course, this did not work. After lots of searching and reworking the 1st-3rd derivatives and countless trial-and-error, I concede that I do not know what I am doing, and need a gratuitous shove in the right direction.

My current tridiagonal matrix for the natural spline looks like $$\left[ \begin{array}{cccc|c} 1 & 0 & 0 & 0 & 0 \\ h_{0} & 2(h_{0}+h_{1}) & h_{1} & 0 & 3(f[x_2,x_1]-f[x_1,x_0]) \\ 0 & h_{1} & 2(h_{1}+h_{2}) & h_{2} & 3(f[x_3,x_2]-f[x_2,x_1]) \\ 0 & 0 & 0 & 1 & 0 \\ \end{array} \right]$$

where $$h=x_{n+1}-x_n$$ and $$f[x_{n+1},x_n] =\frac{y_{n+1}-y_n}{x_{n+1}-x_n}$$ and solves for $$c_0 … c_{n-1}$$

• There are two (equivalent!) formulations of a cubic spline, where you solve for first derivatives in one, and solve for second derivatives in the other. (See e.g. this answer.) Both lead to (different!) tridiagonal systems. "Not-a-knot" just says that the first two pieces are the same cubic polynomial (and similarly for the last two pieces), so you can eliminate common terms from your first two (and last two) equations. – J. M. is a poor mathematician Mar 5 at 11:08
• When solving the tridiagonal system for not-a-knot, should the NAK system have 2 fewer equations than the natural condition? That is what seems sensible to me, as I am rationalizing it as "skipping over" the $x_1, x_{n-2}$ points. I think my confusion stems not from how the system might look, since I believe I have successfully solved how the matrix should look (and I will add that to my question at some point), but more so the implementation relative to the natural spline system. – normal chemist Mar 6 at 10:29

To borrow notation from my previous answer, let

$$S_i(x)=y_i+y_i^{\prime}\left(x-x_i\right)+c_i\left(x-x_i\right)^2+d_i\left(x-x_i\right)^3$$

where $$d_i=\dfrac{y_i^\prime+y_{i+1}^\prime-2s_i}{\left(x_{i+1}-x_i\right)^2}$$ and $$s_i=\dfrac{y_{i+1}-y_i}{x_{i+1}-x_i}$$.

Imposing "not-a-knot" conditions on the left side of the spline means that $$S_0(x)$$ and $$S_1(x)$$ are the same cubic; an equivalent condition is that the third derivatives of both pieces must be the same at $$x=x_1$$ (i.e., $$S_0^{\prime\prime\prime}(x_1)=S_1^{\prime\prime\prime}(x_1)$$). (The derivation for the right side is similar.)

Since $$S_i^{\prime\prime\prime}(x)=6d_i$$ (why?), we have the equation

$$\frac{y_0^\prime+y_{1}^\prime-2s_0}{\left(x_1-x_0\right)^2}=\frac{y_1^\prime+y_2^\prime-2s_1}{\left(x_2-x_1\right)^2}$$

Combine this with the equation (with $$h_i=x_{i+1}-x_i$$):

$$h_1 y_0^{\prime}+2(h_0+h_1)y_1^{\prime}+h_0 y_2^{\prime}=3(h_1 s_0+h_0 s_1)$$

and eliminate $$y_2^{\prime}$$; after some algebra and a little sweat, you should obtain

$$h_1 y_0^{\prime}+2(h_0+h_1)y_1^{\prime}=\frac{(3h_0+2h_1)s_0h_1}{h_0+h_1}$$

For completeness, I will also show the equivalent derivation for the form

$$S_i(x)=y_i+\beta_i\left(x-x_i\right)+\frac{y_i^{\prime\prime}}{2}\left(x-x_i\right)^2+\delta_i\left(x-x_i\right)^3$$

where $$\delta_i=\dfrac{y_{i+1}^{\prime\prime}-y_i^{\prime\prime}}{6h_i}$$. Since $$S_i^{\prime\prime\prime}(x)=6\delta_i$$ (why, again?), we have the equation

$$\frac{y_1^{\prime\prime}-y_0^{\prime\prime}}{x_1-x_0}=\frac{y_2^{\prime\prime}-y_1^{\prime\prime}}{x_2-x_1}$$

which can be combined with

$$h_0 y_0^{\prime\prime}+2(h_0+h_1)y_1^{\prime\prime}+h_1 y_2^{\prime\prime}=6(s_1-s_0)$$

to yield

$$(h_0^2-h_1^2)y_0^{\prime\prime}+(2h_0^2+3h_0 h_1+h_1^2)y_1^{\prime\prime}=6h_0(s_1-s_0)$$