# How to solve a problem with structure similar to a finite difference discretization of the 2D Poisson equation, but with non-symetric coefficients?

Recently, I've been asking about methods to solve a finite difference discretization of the 2D Poisson equation (see here and here) of the form:

$$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$$

The left hand side of this system is often referred to as the 5-point discrete Laplacian. Suppose that left hand side of the linear system is not quite the traditional five-point discrete Laplacian, but instead has the form:

$$\sum_{k,l} a^{ij}_{kl} U_{kl} = a^{ij}_{i-1,j}U_{i-1,j} + a^{ij}_{i+1,j}U_{i+1,j} -a^{ij}_{i,j}U_{i,j} + a^{ij}_{i,j-1}U_{i,j-1} + a^{ij}_{i,j+1}U_{i,j+1} = f_{i,j}$$

with non-negative off-diagonal ($a^{ij}_{k,l} \ge 0$ for all off-diagonal entries $(i,j) \ne (k,l)$) and zero row sums, $$\sum_{kl} a^{ij}_{kl} = 0$$ for all $i,j$. From these conditions, the matrix $A = (a^{ij}_{kl})$ is the negative of an $M$-matrix.

The problem is sparse and diagonally dominant (I would characterize it as "weakly diagonally dominant" because the diagonal term is exactly equal to the sum of the off diagonal elements of the same row). Since it exhibits a similar structure to the five-point Laplacian (except perhaps the possibility of being symmetric?), I figured that there may be some similarities between methods to solve this system of equations and those of the five-point discrete Laplacian. Specifically, I was thinking that one of the following may be true:

1. We can compact the sparsity of the matrix into a thinly banded system (since there are only five variables per equation) and use a direct solver.
2. We can use some sort of Fast Fourier Transform to solve it. (I think this one is highly unlikely though.)
3. We can use some sort of multigrid approach to solve the system iteratively.

In the end, the systems will be very large and a parallel implementation will be necessary. Any feedback on fast methods to approach this problem in parallel would be greatly appreciated.

• @Paul Your stencil notation does need correcting. In the general case, you need four indices, in which case the matrix equation reads $\sum_{kl} a_{ij}^{kl} U_{kl} = f_{ij}$ for each grid point $(i,j)$. Then you can ask for zero row sums $\sum_{kl} a_{ij}^{kl} = 0$ and non-negative off-diagonal entries $a_{ij}^{kl} \ge 0$ for $(k,l) \ne (i,j)$. – Jed Brown Jan 19 '12 at 21:32

## 1 Answer

1. Please take the advice given by the answers to your several previous questions and stop suggesting that banded direct solvers would be appropriate or this problem. They are not, and sparse direct solvers do not attempt to find low-bandwidth orderings, they choose orderings that reduce fill (like nested dissection).

2. FFT is not appropriate for the case of general coefficients.

3. Geometric or (depending on problem structure) algebraic multigrid should work well for this problem, provided you handle a few possibilities correctly. Note that if the coefficients represent extreme anisotropy, then the continuum equations become similar to a bunch of decoupled 1D problems, so isotropic coarsening cannot be expected to work (at least without line smoothers). If the equations represent transport, then you may have to use upwinded interpolants and/or relaxation. Note that second order linear discretizations of transport produce oscillations, so you would likely be using the upwinded form already. If the coefficients are highly variable, you may need to construct low-energy interpolants.

In general, it is much better to say which equations you are solving in which regime with which discretization than it is to make superficial observations about the matrix structure.

• Your answers have been extremely helpful. Do you know of any good paper/book for the FFT in the case of general coefficients, in the context of this problem? – Paul Jan 19 '12 at 19:49
• The equations that I'm solving are derived from solving for the potentials at the nodes of an electric resistor network on a rectangular lattice. I'm using kirchoff's laws, and imposing a constant electric potential at both ends of the network. – Paul Jan 19 '12 at 20:10
• Sorry, that was a typo. FFT is not appropriate for general coefficients. It works for arbitrary forcing, and there are fast direct methods that work for simple coefficient structure, but they would be too inefficient in the general case. – Jed Brown Jan 30 '12 at 13:56
• That's exactly what I suspected. – Paul Jan 30 '12 at 13:58