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}

We cannot apply the Gaussian elimination directly to the entire system. Instead, we rewrite the system as (the first and the one before the last equation are changed):
\begin{eqnarray}
b_0 x_0 + c_0 x_1  &=& d_0 - a_0 x_n,\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}  &=& d_{n-1} - c_{n-1} x_{n},\nonumber\\
c_n x{0} + a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber
\end{eqnarray}

If we omit the last equation, the remaining system of $n$ equations can be written in the matrix form as 

\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-2} & b_{n-2} & c_{n-2}  \\
0 & 0 & 0 & \dots & 0 & a_{n-1} & b_{n-1} \end{matrix}\right)
\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-2} \\ x_{n-1} \end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-2}\\ d_{n-1}\end{matrix}\right) - \left(\begin{matrix} a_0 \\ 0 \\ \vdots \\0 \\ c_{n-1} \end{matrix}\right) x_n.
\label{eq:matrix2}\nonumber\end{equation}
Let's label the matrix and the three vectors in the last equation with $M$, $X$, $D$ and $\Delta$. The equation becomes
$$ M\, X = D - \Delta\, x_n.
$$
Due to the linearity of the problem, we may split the solution into two parts, $X^{(1)}$ and $X^{(2)}$ such that:
$$ X = X^{(1)} + X^{(2)} x_n.
$$
Each of this parts is a solution of a pure three-diagonal system of $n$ equations:
$$ M\, X^{(1)} = D,
$$

and

$$ M\, X^{(2)} =  - \Delta.
$$

These two systems we can easily solve by Gaussian elimination as it was demonstrated before. 

The next step is to compute the solution for $x_n$. For that we use the last equation of the full system that we omitted so far:

$$ c_n x_{0} + a_{n} x_{n-1} +b_{n} x_{n} = d_{n}, $$
$$ \Rightarrow c_n (X^{(1)}_{0} + X^{(2)}_{0} x_n) + a_{n} (X^{(1)}_{n-1} + X^{(2)}_{n-1} x_n) +b_{n} x_{n} = d_{n}, $$
$$ \Rightarrow x_n = \frac{d_{n} - c_n X^{(1)}_{0} - a_{n} (X^{(1)}_{n-1}}{b_n + c_n X^{(2)}_{0}+a_{n} X^{(2)}_{n-1}}. $$
The last equation gives the value of $x_n$. All other $n$ values of $x$ we find trivially by using
$$ X = X^{(1)} + X^{(2)} x_n. $$

Implementation of this algorithm in IDL is coded in threediagonal_periodic_system.pro. The function calls two times  threediagonal_system.pro to calculate $X^{(1)}$ and $X^{(2)}$, computes $x_n$ and combines them into the final solution.

At the input the code takes 4 vectors: a, b, c and d with the meaning as in this text above. Note that all 4 arrays have to be of the same size and none of the elements is allowed to be 0.

Example:

; We form a system of five equations
IDL> aa = [6., 2.,  3., 4., 1.]
IDL> bb = [3., 4., 11., 7., 2.]
IDL> cc = [1., 1.,  1., 3., 3.]

; Set the exact solution
IDL> x0 = FINDGEN(5)

; Compute the right-hand side to complete the system
IDL> dd = FLTARR(5)
IDL> dd[0] = bb[0]*x0[0] + cc[0]*x0[1] + aa[0]*x0[4]
IDL> FOR i = 1, 3 DO dd[i] = aa[i]*x0[i-1] + bb[i]*x0[i] + cc[i]*x0[i+1]
IDL> dd[4] = aa[4]*x0[3] + bb[4]*x0[4] + cc[4]*x0[0]

; Solve the system
IDL> x = threediagonal_periodic_system(aa, bb, cc, dd)

; Compare the solution of the system and the exact solution
IDL> PRINT, x
      0.00000      1.00000      2.00000      3.00000      4.00000
IDL> PRINT, x0
      0.00000      1.00000      2.00000      3.00000      4.00000