Cubic spline interpolation: Natural splines

Last year I published a post on the data interpolation and smoothing using the cubic splines using IDL. (Soon I'll do an update of that post as well.) Here I go a step back and show how the classical spline interpolation works. As before, the derivation follows the book of D.S.G.Pollock (highly recommended to any student interested in the analysis of time-series). I'll first go through some theory general for all the cubic splines, than through the derivation of the natural cubic splines and at the end I'll show some examples. In the next post I'll describe the solution of the clamped cubic splines and do the same examples using that variation of the method.

Three-diagonal periodic system

Three-diagonal periodic system of linear equations can be solved in a similar way as the pure three-diagonal system. The difference is that now we have additional two coefficients in the matrix of the system ($a_0\ne0$, $c_n \ne0$). This form of a system is common when the boundary conditions are periodic (e.g. periodic cubic splines). The solution is again based on the Gaussian elimination, but because of the additional matrix elements, it requires some more steps. The system is strictly diagonal dominant. It can be proved that it has a unique solution.

There is n+1 equation. Our goal is to solve it for n+1 value of $x_i$. The matrix of the system (size $n+1\times n+1$) is of the following form (additional coefficients are shown in bold font):

\begin{equation}\left(\begin{matrix} b_0 & c_0 & 0   & \dots & 0   & 0  & \mathbf{a_0} \\
a_1 & b_1 & c_1 & \dots & 0   & 0  & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots   & \vdots  & \vdots  \\
0 & 0 & 0 & \dots & a_{n-1} & b_{n-1} & c_{n-1}  \\
\mathbf{c_n} & 0 & 0 & \dots & 0 & a_{n} & b_n \end{matrix}\right)
\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1}  \\x_{n}\end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1}  \\d_{n}\end{matrix}\right)
\label{eq:matrix}\nonumber\end{equation}

or in the form of the individual equations ($n+1$ equation with $n+1$ unknown $x_i$):

\begin{eqnarray}
b_0 x_0 + c_0 x_1 + a_0 x_n &=& d_0,\nonumber\\
a_1 x_0 + b_1 x_1 + c_1 x_2 &=& d_1,\nonumber\\
\dots &=& \dots, \nonumber\\
a_i x_{i-1}+ b_i x_i + c_i x_{i+1} &=& d_i,\nonumber\\
\dots &=& \dots ,\nonumber\\
a_{n-1} x_{n-2} +b_{n-1} x_{n-1} + c_{n-1} x_{n} &=& d_{n-1},\nonumber\\
c_n x_{0} + a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber
\end{eqnarray}

Three-diagonal system of linear equations

In this post I'll show a fast and accurate solution for the three-diagonal system of linear equations based on Gaussian elimination. In the matrix notation the system has the following form:

\begin{equation}\left(\begin{matrix} b_0 & c_0 & 0   & \dots & 0   & 0  & 0 \\
a_1 & b_1 & c_1 & \dots & 0   & 0  & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots   & \vdots  & \vdots  \\
0 & 0 & 0 & \dots & a_{n-1} & b_{n-1} & c_{n-1}  \\
0 & 0 & 0 & \dots & 0 & a_{n} & b_n \end{matrix}\right)
\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1}  \\x_{n}\end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1}  \\d_{n}\end{matrix}\right)
\label{eq:matrix00}\nonumber\end{equation}

or in the form of the individual equations ($n+1$ equation with $n+1$ unknown $x_i$):

\begin{eqnarray}
b_0 x_0 + c_0 x_1 &=& d_0,\nonumber\\
a_1 x_0 + b_1 x_1 + c_1 x_2 &=& d_1,\nonumber\\
\dots &=& \dots, \nonumber\\
a_i x_{i-1}+ b_i x_i + c_i x_{i+1} &=& d_i,\nonumber\\
\dots &=& \dots ,\nonumber\\
a_{n-1} x_{n-2} +b_{n-1} x_{n-1} + c_{n-1} x_{n} &=& d_{n-1},\nonumber\\
a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber
\end{eqnarray}

"The Sun is miasma of incandescent plasma"

"Why Does the Sun Shine?" is a funny song on the solar physics originally written by Singer and Zaret. Here it is covered by They Might Be Giants, an alternative band from the 80's, best know for their super-mega-hit "Istanbul (Not Constantinople)". On the same album of children's music TMBG published another song - "Why Does the Sun Really Shine?", an errata to the original one ("the Sun is miasma of incandescent plasma" :) ). Both songs can be heard on YouTube:

Kourganoff graphs

Scharzschild equation provides a tool to easily find and analyze the solution of the radiative transfer cases for some simple source functions. First of all, the lambda operator is clearly linear, so\begin{equation}
\Lambda[a\,f(t) + b\,g(t)] = a\,\Lambda[f(t)] + b\,\Lambda[g(t)].
\end{equation}

Formal solution and Schwarzschild equation

The radiation transport equation (RTE) for time-independent plane-parallel atmosphere has well known form (e.g. Rutten, 2003):
\begin{equation}
\mu\,\frac{\mathrm{d}I_{\nu}}{\mathrm{d} \tau_{\nu}} = I_{\nu} - S_{\nu},
\end{equation}
 where $I_{\nu}$ is specific intensity, $S_{\nu}$ source function, $\tau_{\nu}$ is the radial optical depth and $\mu$ is the cosine of the angle between the line of sight and vertical axis. Mathematically, this is a linear first order differential equation with constant coefficients. Physically, it describes propagation of radiation through an interacting medium. Properties of the medium are hidden in the source function and the optical depth. However, to make the first step toward the solution of RTE, we will assume that the optical depth is independent variable and that the source function is a know function of it.

Gauss quadrature and Legendre polynomials

Numerical integration is among the most common tasks in astrophysics. Simple formulae, like trapezoidal, Newton-Cotes or Simpson, are often not enough accurate. Gauss quadrature provides a more accurate solution, but its implementation is a bit more difficult. Here I will demonstrate the implementation of that algorithm in IDL. The algorithm includes finding the coefficients of Legendre polynomials and their zeros. Both tasks are explicitly solved in the following as well.

But let's go step by step. Our task is to solve numerically the integral:
$$I = \int_a^b f(x)\, \mathrm{d} x.$$