# $AX=B$: How to solve for $X$ if elements of matrix A are matrices

Objective: I am trying to solve for $C$ in 2D space (x,y) and time from following PDE. $$\text{PDE: }\frac{\partial C}{\partial t} + \nabla\left(v.C - D\nabla{C} \right)= \alpha.C$$ Method: I discretized the above equation in 2D using finite-difference implicit scheme and as a result I get the following discretized equation.

$$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}$$ where, $p_1, p_2, p_3, p_4, p_5$ are constant matrices in time. They are matrices instead of scalars because I am assuming $v$ and $D$ vary in space $(x,y)$ but constant in time i.e. $v(x,y)$ and $D(x,y)$. The variables $p_1, p_2, p_3, p_4, p_5$ are shown in the snapshot below: I can solve for $C$ at each time by solving the system of equations as $A^{n+1}.C^{n+1}=C^{n}$.

Issue: Since $p_1, p_2, p_3, p_4, p_5$ are matrices, I am not sure how to develop the matrix $A$. As an example for a system of $3\times3$ grid size, matrix $A$ will have the form as shown below. In that case how can I form matrix $A$ from $p_1, p_2, p_3, p_4, p_5$, each of which is a matrix of size $3\times3$? I would really appreciate if someone can guide how to solve this. Thanks. • Is $C$ scalar or vector? – Kirill Jul 10 '15 at 11:19
• Can you invert $C^{n+1}$? If yes, the isn't the answer $A^{n+1}=C^{n}.(C^{n+1})^{-1}$ – dpmcmlxxvi Jul 10 '15 at 16:28
• $C^{n+1}$ is an unknown vector that I'm trying to compute. It contains value of each index at time $n+1$. $C^{n}$ at right hand side is a known vector. – user5510 Jul 10 '15 at 18:21
• @Pupil Sorry I meant to write can you invert $A^{n+1}$. – dpmcmlxxvi Jul 10 '15 at 22:16
• As Bill Barth pointed out in his answer, the fundamental problem is that your discretization is not correct (or at least, not standard). For every $n$, you should treat $C^n$ as a vector of unknowns, so you end up with an equation of the form $AC^{n+1} = B C^{n}$, where $A,B$ are matrices (of finite differences, involving the coefficients $\nu$, $D$, $\alpha$). Deriving the right entries and getting them in the right place is a bit tedious but standard, and I would recommend looking at a standard textbook on finite difference methods (such as Morton/Myers or one of Uri Ascher's books). – Christian Clason Jul 11 '15 at 8:58

What is usually done when you have a matrix-like variable due to the dimensions of your problem and the indexing is to linearize the index. I.e., if you have $C_{i,j}$ you would replace that with $\hat{C}_k=C_{i,j}$ where $k=i\times N+j$ and $N$ is the number of points on one side of your $N\times N$ grid of points that you are doing finite differences on. This makes it easy to convert your kind of strange linear system into the standard form.
This means that you might end up representing $C$ as a single-indexed array in your program instead of treating it with double-indexing. There are other ways to accomplish this once you've written down the right linear system, but this method is compatible with most linear solver libraries which expect vectors of unknowns and right-hand sides.
• Where are you suggesting to convert those indices? In the final $A.C^{n+1}=C^{n}$ equation they're already in that form. – user5510 Jul 10 '15 at 18:24
• @Pupil, or you need to express the forms of the $p_{*}$ matrices so that we can understand your problem better. Do they depend on $i$ and $j$? – Bill Barth Jul 10 '15 at 18:31
• @Pupil, everything you wrote for the $p$s appears to be a scalar. What's your problem, again? – Bill Barth Jul 12 '15 at 0:30
• @Pupil, no, no, no, at each point $i,j$, they are scalars just like $C$. It's only a "matrix" when you let them run over the indices. But the indices coincide. Write it all down with a linearized index like I said, and you'll see it. You'll have an equation for $C_k$ with some $k+1$ and $k-1$ terms. That gives you the general form of a matrix row. – Bill Barth Jul 13 '15 at 20:51