# ADR equation implicit solution: Penta-diagonal matrix for a 2D $N\times N$ system

Objective: I am trying to simulate the following advection-diffusion-reaction equation in 2D space (x,y) and time. \begin{align} \text{ADR Equation: }\frac{\partial C}{\partial t} + \nabla\left(v.C - D\nabla{C} \right)= \alpha.C \end{align} I discretized the above ADR equation in 2D using finite-difference implicit scheme and as a result I get the following discretized equation.

\begin{align} p_1C^{n+1}_{i,j-1}+p_2C^{n+1}_{i-1,j}+p_3C^{n+1}_{i,j}+p_4C^{n+1}_{i+1,j}+p_5C^{n+1}_{i,j+1} = C^{n}_{i,j} \end{align} where, $p_1, p_2, p_3, p_4, p_5$ are constants in time.

I want to solve this as a system of equations using $A^{n+1}.C^{n+1}=C^{n}$, with no-flow i.e. $C=0$ outside the boundary domain. Here, $A^{n+1}$ would be a penta-diagonal, symmetric (not sure about this) and a diagonally dominant matrix. I have derived matrix $A$ for $2\times2$, $3\times3$ and $4\times4$ systems. For example, below you can see matrix $A$ for $3\times3$ and $4\times4$ systems, respectively.

Issue:

1. I am not sure if the form of matrix $A$ I have derived is correct because as per my understanding it should be symmetric, however, it's not as per my derivations.
2. Owing to an unsymmetric form of the matrix $A$ I need help to efficiently form $A$ for $N\times N$ system.

Would appreciate if someone could use their awesome numerical skills to answer these issues.

• Doesn't pentagonal usually mean that the nonzero elements are in the shape of a pentagon inside the matrix? It's not common, so I'm not sure what you mean. The matrix would normally have nonzero diagonals $\{-n,-1,0,1,n\}$, like $T\otimes I+I\otimes T$ in terms of Kronecker product, but that's not pentagonal. – Kirill Jul 8 '15 at 5:51
• @Kirill: Please ignore the term pentagonal. But, matrix A that I am talking about looks like the one shown in eqn. 17 here geodynamics.usc.edu/~becker/teaching/557/problem_sets/… – user5510 Jul 8 '15 at 6:13
• Or the A matrix for eqn. 6.4.45 in this link: www3.msiu.ru/~belova/compmod/Crank_nic.pdf – user5510 Jul 8 '15 at 6:24
• Or maybe the correct word is 'penta-diagonal' – user5510 Jul 8 '15 at 6:29
• @Pupil, You will not get symmetric matrix unless you use symmetric discretization schemes for your convection term (eg Central difference, compact scheme etc.). In most of the practice we rarely use central schemes for convection because of its numerical dispersion nature . – AGN Jul 18 '15 at 0:50

This is mostly just a matter of assigning an index to each vertex. If you have $n\times n$ unknowns $C_{i,j}$, $0\leq i,j\leq n-1$, then you would typically map the variable corresponding to the vertex $(i,j)$ to something like an index $I(i,j) = i + n j$. So each variable $C_{i,j}$ that is indexed by a pair of numbers $(i,j)$ is equivalent to the variable $C_{I(i,j)}$ indexed by a single number.

With this mapping, you rewrite the equations as $$p_1 C_{I(i,j)} + p_{2:3} C_{I(i\pm1,j)} + p_{4:5} C_{I(i,j\pm1)} = \cdots,$$

This means that the $I(i,j)$-th row (the $(i+jn)$-th row) of the matrix $A^{n+1}$ has the entries $p_{1,2,3,4,5}$ in columns $I(i,j)$, $I(i\pm1,j)$, $I(i,j\pm1)$.

With this indexing method the nonzero entries will all fall into the diagonals $\pm n$ (coefficients $p_4$, $p_5$), $\pm1$ ($p_2$, $p_3$), and $0$ ($p_1$). The main thing to bear in mind, when constructing such a matrix knowing its diagonal elements, is that the diagonal elements for non-existent boundary variables such as $I(i,n)$, $I(i,-1)$, $I(-1,j)$, $I(n,j)$ should be zeroed out.

So for a $3\times 3$ grid, you would get the matrix (I dropped the minus signs here) $$\begin{pmatrix} p_1 & p_2 & 0 & p_4 & 0 & 0 & 0 & 0 & 0 \\ p_3 & p_1 & p_2 & 0 & p_4 & 0 & 0 & 0 & 0 \\ 0 & p_3 & p_1 & 0 & 0 & p_4 & 0 & 0 & 0 \\ p_5 & 0 & 0 & p_1 & p_2 & 0 & p_4 & 0 & 0 \\ 0 & p_5 & 0 & p_3 & p_1 & p_2 & 0 & p_4 & 0 \\ 0 & 0 & p_5 & 0 & p_3 & p_1 & 0 & 0 & p_4 \\ 0 & 0 & 0 & p_5 & 0 & 0 & p_1 & p_2 & 0 \\ 0 & 0 & 0 & 0 & p_5 & 0 & p_3 & p_1 & p_2 \\ 0 & 0 & 0 & 0 & 0 & p_5 & 0 & p_3 & p_1 \end{pmatrix}$$

• Thanks for your answer Kirill. One thing that concerns me is that this matrix doesn't looks symmetric. In fact, that was one reason I asked this question here because when I personally tried I was getting exactly this matrix, though in N X N domain. – user5510 Jul 8 '15 at 8:18
• @Pupil It has to be $n^2\times n^2$ because that's how many variables there are. It's not really going to be symmetric because the operator $\nabla\cdot (v C-D\nabla C)$ you're discretizing is itself not symmetric. I don't see how you can expect this matrix to be symmetric. – Kirill Jul 8 '15 at 8:52
• Yeah, the matrix will be N^2 X N^2, however, by domain I meant the extent of the studied area in X and Y domain, which would be N grids in X direction and N grids in Y direction. – user5510 Jul 8 '15 at 9:03
• @Pupil The general form of the matrix $A$ is straightforward, and is already explained in my answer. For constructing the matrix, one usually uses something like spdiags, or diags for specifying diagonal elements only. You can also just construct it by element-by-element. – Kirill Jul 9 '15 at 4:39