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: enter image description here

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. enter image description here

  • $\begingroup$ Is $C$ scalar or vector? $\endgroup$
    – Kirill
    Jul 10, 2015 at 11:19
  • $\begingroup$ Can you invert $C^{n+1}$? If yes, the isn't the answer $A^{n+1}=C^{n}.(C^{n+1})^{-1}$ $\endgroup$
    – dpmcmlxxvi
    Jul 10, 2015 at 16:28
  • $\begingroup$ $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. $\endgroup$
    – user5510
    Jul 10, 2015 at 18:21
  • $\begingroup$ @Pupil Sorry I meant to write can you invert $A^{n+1}$. $\endgroup$
    – dpmcmlxxvi
    Jul 10, 2015 at 22:16
  • 1
    $\begingroup$ 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). $\endgroup$ Jul 11, 2015 at 8:58

1 Answer 1


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.

  • $\begingroup$ Where are you suggesting to convert those indices? In the final $A.C^{n+1}=C^{n}$ equation they're already in that form. $\endgroup$
    – user5510
    Jul 10, 2015 at 18:24
  • $\begingroup$ @Pupil, no, that's your time index. You need to linearize over space, too. $\endgroup$
    – Bill Barth
    Jul 10, 2015 at 18:25
  • $\begingroup$ @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$? $\endgroup$
    – Bill Barth
    Jul 10, 2015 at 18:31
  • 1
    $\begingroup$ @Pupil, everything you wrote for the $p$s appears to be a scalar. What's your problem, again? $\endgroup$
    – Bill Barth
    Jul 12, 2015 at 0:30
  • 2
    $\begingroup$ @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. $\endgroup$
    – Bill Barth
    Jul 13, 2015 at 20:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.