Cubic spline interpolation: Periodic splines

If the unknown function is periodic, then the periodicity may be used to specify the boundary conditions. Before (in the case of the natural splines) we had $n+1$ value of x and y and as a result we obtained values of the cubic polynomial at $n$ intervals.

Now our function is periodic with period $p$, thus $S(x) = S(x+p)$. In practice, we have an additional data point $x_{n+1}$, such that
$$x_{n+1} = x_0 + p,$$
and
$$y_{n+1} = y_0.$$

The total number of data points is $n+2$ and the total number of segments (and cubic splines) between them is $n+1$.


The smoothness requirement at the nodes is exactly the same as before, so we can rewrite the equations (21):
\begin{equation}
q_i = c_{i-1}h_{i-1}+c_{i}w_{i} + c_{i+1}h_{i},
\label{eq:qwsystem}
\end{equation}
where
\begin{equation}
w_i = 2(h_i + h_{i-1}) = 2(x_{i+1} - x_i + x_i - x_{i-1}) = 2(x_{i+1}-x_{i-1}),
\label{eq:w}
\end{equation}
and
\begin{equation}
q_i = y_{i+1}\frac{3}{h_i} + y_i \left(-\frac{3}{h_i}-\frac{3}{h_{i-1}} \right) +y_{i-1}\frac{3}{h_{i-1}},
\label{eq:q}
\end{equation}

where $i = 0, \dots, n$.

The equation for $i = 0$ is

\begin{equation}
q_0 = c_{-1}h_{-1} + c_{0}w_{0} + c_{1}h_{0}.
\end{equation}

Due to the periodicity, the -1-st interval is exactly the same as n-th interval, so we may write

\begin{equation} q_0 = c_{n}h_{n} + c_{0}w_{0} + c_{1}h_{0}. \end{equation}

The equation for $i = n$ reads:

\begin{equation}
q_n = c_{n-1}h_{n-1}+c_{n}w_{n} + c_{n+1}h_{n}.
\end{equation}

Again, the periodicity wraps up the system, n+1-st segment is equal to the 0-th segment, so:

\begin{equation} q_n = c_{n-1}h_{n-1}+c_{n}w_{n} + c_{0}h_{n}. \label{eq:qn} \end{equation}

Here we implicitly introduced the boundary condition for the periodic splines that requires that the interpolation function remains smooth at the end-points. The condition can be formulated as:

Periodic splines:
\begin{equation} S''_0(x_n+1) = S_{n}''(x_n), \;\;\;\;\;\;\;S_{n}''(x_n) = 2c_n = 0. \label{eq:natural} \end{equation}

In the matrix form, the system of equation can now be written as
\begin{equation}\left(\begin{matrix} q_0 \\ q_1 \\ \vdots \\ q_{n-1}  \\q_{n}\end{matrix}\right)=
\left(\begin{matrix} w_0 & h_0 & 0   & \dots & 0   & 0  & h_n \\
h_0 & w_1 & h_1 & \dots & 0   & 0  & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots   & \vdots  & \vdots  \\
0 & 0 & 0 & \dots & h_{n-2} & w_{n-1} & h_{n-1}  \\
h_n & 0 & 0 & \dots & 0 & h_{n-1} & w_n \end{matrix}\right)
\left(\begin{matrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1}  \\c_{n}\end{matrix}\right),
\label{eq:matrix1}\end{equation}.

The matrix is strictly diagonally dominant, so the system has a unique solution and in can be solved by Gaussian elimination. There is an earlier post describing how to solve this system in an elegant way (including and IDL implementation of the algorithm).

My IDL implementation of the periodic cubic splines for a given set of data points (x, y) is coded in the function csp_interpolation.pro. To evaluate the interpolation function (once the coefficients are known), one may use csp_evaluate.pro.

Example:

; The exact function for comparison.
nxa = 241
period = 2.
xa = DINDGEN(nxa)/100-1.2
ya = SIN(2.*!pi*xa/period)  
     
; The 'measurement' samples the function at 9 points (the array must be
; monotonically increasing).
x = DOUBLE([-1., -0.8, -0.6, -0.45,  0.0, 0.3, 0.5, 0.6, 1.])
nx = N_ELEMENTS(x)
      
; Force the exact periodicity, y[0] = y[nx-1]
y = [SIN(2*!pi*x[0:nx-2]/period), SIN(2*!pi*x[0]/period)]
        
; Compute the coefficients of the approximative function S and evaluate S 
; at the same grid used to compute the exact function.
nx0 = nxa
cp = CSP_INTERPOLATION(x, y)
x0 = FINDGEN(120)/4.-3
y0 = CSP_EVALUATE(x0, x, cp) 

Fig. Pluses show the nine data points in which the sine function is sampled. All the samples are from the [-1, 1] interval. The black solid line shows the exact function. The periodic cubic splines are computed using the data points and the interpolation function is evaluated at the entire interval. The periodicity of the exact function is of course fully replicated.