I'm trying to implement a model I found in a paper, but there is something I do not understand. The authors say they use TDMA to solve their equations; however, they use a 3rd order upwind biased approximation to discretize the convective terms. This leads to an almost tridiagonal matrix (with coeff. for the terms $i-1$, $i$, and $i+1$), but with an extra diagonal line for the $i-2$ terms.
Is it possible to use a TDMA solver with a matrix like this? Is there some procedure to tridiagolise?
I would assume, that your $N\times N$ system has the following form:
$$ \underbrace{\begin{pmatrix} b_1 & c_1 & & & & 0 \\ a_1 & b_2 & c_2 & & & \\ s_1 & a_2 & b_3 & c_3 & & \\ & s_2 & a_3 & b_4 & \ddots & \\ & & \ddots & \ddots & \ddots & c_{N-1}\\ 0 & & & s_{N-2} & a_{N-1} & b_N \end{pmatrix}}_{M} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \vdots\\ x_{N} \end{pmatrix}= \begin{pmatrix} d_1\\ d_2\\ d_3\\ d_4\\ \vdots\\ d_{N} \end{pmatrix} \tag{1} \label{fourDiagSystem} $$ So, our matrix $M$ is a four-diagonal matrix (lower bandwidth = 2, upper bandwidth = 1), where $\{s_1,\ldots,s_{N-2} \}$ is the "additional" (to the classic triadiagonal case) lower diagonal.
So, if all $s_i$ are zeroes, this matrix can be solved using the classical tridiagonal algorithm (Thomas). There are several ways to look at this problem:
Since the question asks directly about using Thomas algorithm, I would go only into the third point.
Assume (without loss of generality) $N$ being even. Then, $M$ from $\eqref{fourDiagSystem}$ can be rewritten, as follows: $$ M=\begin{pmatrix} \underbrace{\begin{pmatrix} b_1 & c_1 \\ a_1 & b_2 \end{pmatrix}}_{B_{1}} & \underbrace{\begin{pmatrix} 0 & 0 \\ c_2 & 0\end{pmatrix}}_{C_{1}} & & \\ \underbrace{\begin{pmatrix} s_1 & a_2 \\ 0 & s_2 \end{pmatrix}}_{A_{1}} & \underbrace{\begin{pmatrix} b_3 & c_3 \\ a_3 & b_4 \end{pmatrix}}_{B_{2}} & \underbrace{\begin{pmatrix} 0 & 0 \\ c_4 & 0\end{pmatrix}}_{C_{2}} & \\ & \underbrace{\begin{pmatrix} s_3 & a_4 \\ 0 & s_4 \end{pmatrix}}_{A_{2}} & \underbrace{\begin{pmatrix} b_5 & c_5 \\ a_5 & b_6 \end{pmatrix}}_{B_{2}} & \ddots \\ & \ddots & \ddots & \underbrace{\begin{pmatrix} 0 & 0\\ c_{N-2} & 0 \end{pmatrix}}_{C_{N/2-1}} \\ & & \underbrace{\begin{pmatrix} s_{N-3} & a_{N-2}\\ 0 & s_{N-2} \end{pmatrix}}_{A_{N/2-1}} & \underbrace{\begin{pmatrix} b_{N-1} & c_{N-1}\\ a_{N-1} & b_{N} \end{pmatrix}}_{B_{N/2}} \end{pmatrix} $$
or, removing the clarifying annotations, $\eqref{fourDiagSystem}$ becomes:
$$ \underbrace{\begin{pmatrix} B_1 & C_1 & & &0\\ A_1 & B_2 & C_2 &\\ & A_2 & B_3 & \ddots &\\ & & \ddots & \ddots & C_{N/2-1}\\ 0 & & & A_{N/2-1} & B_{N/2} \end{pmatrix}}_{M_\text{block}} \begin{pmatrix} Y_1\\ Y_2\\ Y_3\\ \vdots\\ Y_{N/2} \end{pmatrix}= \begin{pmatrix} D_1\\ D_2\\ D_3\\ \vdots\\ D_{N/2} \end{pmatrix} \tag{2} \label{triBlockDiagSystem} $$
Now, in $\eqref{triBlockDiagSystem}$, $M_\text{block}$ is simply a $2\times 2$ block representation of $M$ in $\eqref{fourDiagSystem}$. $$ B_i=\begin{pmatrix} b_{2i-1} & c_{2i-1}\\ a_{2i-1} & b_{2i} \end{pmatrix},\quad A_j=\begin{pmatrix} s_{2j-1} & a_{2j}\\ 0 & s_{2j} \end{pmatrix},\quad C_j=\begin{pmatrix} 0 & 0\\ c_{2j} & 0 \end{pmatrix},\\ i=1,\ldots,\frac{N}{2},\quad j=1,\ldots,\frac{N}{2}-1 $$
Now, you can apply the block-Thomas algorithm to the $\eqref{triBlockDiagSystem}$, in which you can also take advantage of the special structure of your $C_j$ and $A_j$ blocks.
For the details on block-Thomas, you can refer to Chapter 2.5 of MA/CSC 580: Numerical Analysis I course by R. E. White from North Carolina, but they are effictively standard non-blocked Thomas where operations with scalars are changed to their counterparts with matrices (scalar-scalar product and scalar division become matrix-matrix product and inverse, or, more stable small, $2\times 2$ in your case, system solve).
Also, see the discussion in this related CompSci question on block-tridiagonal algorithm.
