# TDMA with 3rd order upwind scheme

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?

• what paper are you referring to, exactly? – GoHokies Sep 5 '18 at 11:12
• Can you write what the structure of the matrix looks like? And the equations that you’re solving? – Charles Sep 5 '18 at 14:28

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:

1. Explicitly derive the expressions for this four-diagonal case, similar to the derivation of the Thomas algorithm.
2. The derivation in the previous "approach" can go along Crout's method. For example, see Numerical Recepies (Sections 2.3, 2.4).
3. Use block-tridiagonal algorithm (Thomas algorithm for block matrices) with block size of $$2\times 2$$.

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.