Perturbed cases

We consider the following system:

\[\begin{split}\begin{bmatrix} c_0 & u_0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & l_0 \\ l_1 & c_1 & u_1 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & l_2 & c_2 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & c_{i-1} & u_{i-1} & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & \cdots & l_{i } & c_{i } & u_{i } & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & l_{i+1} & c_{i+1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & c_{n-3} & u_{n-3} & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & l_{n-2} & c_{n-2} & u_{n-2} \\ u_{n-1} & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & l_{n-1} & c_{n-1} \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{i-1} \\ x_{i } \\ x_{i+1} \\ \vdots \\ x_{n-3} \\ x_{n-2} \\ x_{n-1} \end{bmatrix} = \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ \vdots \\ q_{i-1} \\ q_{i } \\ q_{i+1} \\ \vdots \\ q_{n-3} \\ q_{n-2} \\ q_{n-1} \end{bmatrix},\end{split}\]

where, compare with the original tri-diagonal system, the top-right and the bottom-left corners have non-zero values (\(l_0\) and \(u_{n - 1}\), respectively), which appears in particular when periodic boundary conditions are imposed.

Since we have

\[ \begin{align}\begin{aligned}c_0 x_0 + u_0 x_1 & = q_0 - l_0 x_{n - 1},\\l_{n - 2} x_{n - 3} + c_{n - 2} x_{n - 2} & = q_{n - 2} - u_{n - 2} x_{n - 1},\end{aligned}\end{align} \]

at the \(0\)-th and the \(n - 2\)-th rows, we reformulate the system as

\[\begin{split}\begin{bmatrix} c_0 & u_0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ l_1 & c_1 & u_1 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & l_2 & c_2 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & c_{i-1} & u_{i-1} & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & \cdots & l_{i } & c_{i } & u_{i } & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & l_{i+1} & c_{i+1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & c_{n-4} & u_{n-4} & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & l_{n-3} & c_{n-3} & u_{n-3} \\ 0 & 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & l_{n-2} & c_{n-2} \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{i-1} \\ x_{i } \\ x_{i+1} \\ \vdots \\ x_{n-4} \\ x_{n-3} \\ x_{n-2} \end{bmatrix} = \begin{bmatrix} q_0 - l_0 x_{n-1} \\ q_1 \\ q_2 \\ \vdots \\ q_{i-1} \\ q_{i } \\ q_{i+1} \\ \vdots \\ q_{n-4} \\ q_{n-3} \\ q_{n-2} - u_{n-2} x_{n-1} \end{bmatrix}.\end{split}\]

Note that we kick out \(x_{n-1}\) from the left-hand side, and the resulting system is now tri-diagonal with the size of \(n - 1\). For notational simplicity, hereafter we write this system as

\[A_{ij} x_j = b_i.\]

Although this is a purely tri-diagonal system (i.e., without perturbed elements), we cannot perform forward sweep as an unknown \(x_{n - 1}\) exists in the right-hand-side vector. To circumvent this issue, we split the system into two sub problems:

\[ \begin{align}\begin{aligned}A_{ij} x_j^{\prime} & = b_i^{\prime},\\A_{ij} x_j^{\prime\prime} & = b_i^{\prime\prime},\end{aligned}\end{align} \]

where

\[\begin{split}b_i^{\prime} = \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ \vdots \\ q_{i-1} \\ q_{i } \\ q_{i+1} \\ \vdots \\ q_{n-4} \\ q_{n-3} \\ q_{n-2} \end{bmatrix}, b_i^{\prime\prime} = \begin{bmatrix} - l_0 \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ - u_{n-2} \end{bmatrix},\end{split}\]

which satisfy

\[ \begin{align}\begin{aligned}& x_j = x_j^{\prime} + x_{n-1} x_j^{\prime\prime},\\& b_i = b_i^{\prime} + x_{n-1} b_i^{\prime\prime},\end{aligned}\end{align} \]

and thus the solution \(x_j\) is found as the superposition of the solution of these two sub systems.

To find \(x_{n - 1}\), we look at the \(n - 1\)-th row of the original system:

\[u_{n - 1} x_0 + l_{n - 1} x_{n - 2} + c_{n - 1} x_{n - 1} = q_{n - 1}.\]

By assigning

\[ \begin{align}\begin{aligned}x_0 & = x_0^{\prime} + x_{n - 1} x_0^{\prime\prime},\\x_{n - 2} & = x_{n - 2}^{\prime} + x_{n - 1} x_{n - 2}^{\prime\prime},\end{aligned}\end{align} \]

we notice

\[\left( u_{n - 1} x_0^{\prime} + l_{n - 1} x_{n - 2}^{\prime} \right) + \left( u_{n - 1} x_0^{\prime\prime} + l_{n - 1} x_{n - 2}^{\prime\prime} + c_{n - 1} \right) x_{n - 1} = q_{n - 1},\]

and thus

\[x_{n - 1} = \frac{ q_{n - 1} - u_{n - 1} x_0^{\prime} - l_{n - 1} x_{n - 2}^{\prime} }{ c_{n - 1} + u_{n - 1} x_0^{\prime\prime} + l_{n - 1} x_{n - 2}^{\prime\prime} }.\]

To summarize, the perturbed system is solved as follows.

1. Prepare the right-hand-side vector of a sub system \(b_i^{\prime\prime}\):

   src/tri_diagonal_solver.c

   for (size_t i = 0; i < size - 1; i++) {
     w[i] =
         i ==        0 ? - l[       0]
       : i == size - 2 ? - u[size - 2]
       : 0.;
   }
   

2. Solve two sub systems:

   Forward sweep:

   src/tri_diagonal_solver.c

   v[0] = u[0] / c[0];
   q[0] = q[0] / c[0];
   w[0] = w[0] / c[0];
   

   src/tri_diagonal_solver.c

   for (size_t i = 1; i < size - 2; i++) {
     // assume positive-definite system
     //   to skip zero-division checks
     const double val = 1. / (c[i] - l[i] * v[i - 1]);
     v[i] = val *  u[i];
     q[i] = val * (q[i] - l[i] * q[i - 1]);
     w[i] = val * (w[i] - l[i] * w[i - 1]);
   }
   

   src/tri_diagonal_solver.c

   // NOTE: do the same thing but consider singularity (degeneracy)
   const double val = c[size - 2] - l[size - 2] * v[size - 3];
   if (DBL_EPSILON < fabs(val)) {
     q[size - 2] = 1. / val * (q[size - 2] - l[size - 2] * q[size - 3]);
     w[size - 2] = 1. / val * (w[size - 2] - l[size - 2] * w[size - 3]);
   } else {
     // singular
     q[size - 2] = 0.;
     w[size - 2] = 0.;
   }
   

   Backward substitution:

   src/tri_diagonal_solver.c

   for (size_t i = size - 3; ; i--) {
     q[i] -= v[i] * q[i + 1];
     w[i] -= v[i] * w[i + 1];
     if (0 == i) {
       break;
     }
   }
   

3. Couple the result of two sub systems:

   src/tri_diagonal_solver.c

   const double num = q[size - 1] - u[size - 1] * q[0] - l[size - 1] * q[size - 2];
   const double den = c[size - 1] + u[size - 1] * w[0] + l[size - 1] * w[size - 2];
   q[size - 1] = fabs(den) < DBL_EPSILON ? 0. : num / den;
   for (size_t i = 0; i < size - 1; i++) {
     q[i] = q[i] + q[size - 1] * w[i];
   }