# Understanding the matrix equation used to solve 2-D Poisson Equation with non-uniform grid

I am working on a module in my class that is meant to teach us about Jacobi/Gauss-Seidel methods of solving matrix equations in conjunction with Poisson's equation, and am having trouble with setting up my matrix.

$$\nabla^2 \phi = 4 \pi \rho$$

with finite difference stencil:

$$\frac{\phi_{i-1,j} - 2\phi_{i,j} + \phi_{i+1,j}}{\Delta x^2} + \frac{\phi_{i,j+1} -2\phi_{i,j} + \phi_{i,j+1}}{\Delta y^2} = 4\pi\rho_{i,j}$$

The goal is to transform this stencil into a matrix equation in two different ways. Once with Dirichlet boundary conditions and once with periodic boundary conditions.

$$\phi_{0,j} = \phi_{N,j}, \forall j\\ \phi_{i, 0} = \phi_{i,N}, \forall i\\$$ and for Direchlet on the boundaries: $$\phi_{i,j} = 1$$

I successfully programmed an algorithm using both Jacobi or Gauss Seidel method to solve a matrix, and now I am trying to understand how to put this finite difference stencil into a matrix. I decided to tackle the case of periodic boundary conditions first:

My work so far

Since I am specifically told not to assume $\Delta x = \Delta y$, I rearranged the above equation to find my matrix coefficients with the following substitutions:

$$a = 1/\Delta x^2 \\ b = 1/\Delta y^2 \\ c = a + b\\$$ gives $$\frac{1}{4\pi}[b\phi_{i,j-1} + a\phi_{i-1,j} - 2c\phi_{i,j} + a\phi_{i+1,j} + b\phi_{i,j+1} ]$$

My best attempt at arranging this equation into a matrix looks like:

$$\mathbf{A} = \begin{bmatrix} -2c & a & b & 0 & \dots & b & a \\ a & -2c & a & b & 0 & \dots & b \\ b & a & -2c & a & b & 0 & \dots \\ 0 & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & \\ b & 0 & \dots & b & a & -2c & a\\ a & b & 0 & \dots & b & a & -2c \\ \end{bmatrix} \\$$

Where the solution vector must be ordered as follows:

\begin{align} & (i_{min}, j_{min})\\ & (i_{min} + 1, j_{min}) \\ & (i_{min} + 2, j_{min}) \\ & \vdots \\ & (i_{max}, j_{min})\\ & (i_{min}, j_{min} +1) \\ & (i_{min} + 1, j_{min} + 1) \\ & \vdots \\ & (i_{max} - 1, j_{max}) \\ & (i_{max}, j_{max}) \end{align}

Note: This is my first matrix in TeX, so forgive me if it is a bit ugly.

So, the matrix equation $\frac{ab}{4\pi}\mathbf{A}\mathbf{\phi} = \mathbf{\rho}$ should yield a system of equations that makes sense to the finite difference stencil, right?

Well when I test it out, I get about half of what I'm looking for and then some stuff that doesn't make sense even when considering periodic boundary conditions. For example I will expand the matrix equation for $\rho_{2,0}$:

$$\rho_{2,0} = \frac{1}{4\pi}[b \phi_{(0,0)} + a\phi_{(1,0)} -2c\phi_{(2,0)} + a\phi_{(3,0)} + b\phi_{(4,0)}]$$

but this should, by the finite difference stencil be:

$$\rho_{2,0} = \frac{1}{4\pi}[b \phi_{(2,-1)} + a\phi_{(1,0)} -2c\phi_{(2,0)} + a\phi_{(3,0)} + b\phi_{(2,1)}]$$

So it looks almost right, but no matter how I rearrange my matrix some terms are wrong. I have tried multiple arrangements, but this one gets the most correct terms when I start expanding.

Where am I going wrong?

Some remarks that are to large for a comment, but not really an answer.

An important information is missing, but I assume that you are working on the unit square $Ω=I_x×I_y=[0,1]^2$. And $I_x$ is divided in equidistant subintervals $$0=t^x_0≤t^x_1≤…≤t^x_{m^x+1}=1,$$ with $t^x_{i+1}-t^x_{i}=Δx$, and $I_y$ is divided in the same way.

(I personally like the notation with $m$ representing the inner points and $0$ and $m+1$ are boundary points, but use whatever you like best.)

So your grid looks something like this:

(0,1)                (1,1)
15--16--17--18--19
|   |   |   |   |
10--11--12--13--14
|   |   |   |   |
5---6---7---8---9
|   |   |   |   |
0---1---2---3---4
(0,0)                (1,0)


The numbering you introduced is often called "lexicographical numbering" in the literature. Here in this example $m_x = 3$, $m_y=2$. The matrix size is $n=(m_x+2)(m_y+2)$.

Extracting $\frac{ab}{\cdot}$ from the equation does not give you anything, except a lot of confusion. You wrote this: $$\frac{ab}{4\pi}[a^{-1}\phi_{i,j-1} + b^{-1}\phi_{i-1,j} - 2\phi_{i,j}(a^{-1} + b^{-1}) + b^{-1}\phi_{i+1,j} + a^{-1}\phi_{i,j+1}]$$ which is nothing else than $$\frac{1}{4\pi}[b\phi_{i,j-1} + a\phi_{i-1,j} - 2\phi_{i,j}(b + a) + a\phi_{i+1,j} + b\phi_{i,j+1}]$$

which looks more like the original stencil, as $b$ is simply the distance between two points in $y$-direction. And that is indicated by a change in the $j$-index.

I don't want to introduce a new notation here, but if I were you I would have used $Δ\tilde{x}=\frac{1}{Δx^2}$ to save space (if necessary).

About the stencil: The stencil represents the matrix row of a node in the centre of the grid. For nodes at the boundary you have the boundary conditions (see below.)
E.g. for row $k:=7$, representing point 7 you have

(0,1)                (1,1)
15--16--17--18--19
|   |   |   |   |
10--11--12--13--14
|   |   ∥   |   |
5---6===7===8---9
|   |   ∥   |   |
0---1---2---3---4
(0,0)                (1,0)


which results in: $$\frac{1}{4\pi}[b\phi_{2} + a\phi_{6} - 2\phi_{7}(b + a) + a\phi_{8} + b\phi_{12}]$$

This will give you the matrix row:

$$\begin{matrix} … & 0 &b & 0&…&0& a & 2(a+b)&a & 0 & … & 0 &b & 0&… \\ & &2 & &&& 6 & 7&8 & & & &12 & & \\ & &k-(m_x+2) & &&& k-1 & k&k+1 & & & &k+(m_x+2) & & \end{matrix}$$

First row: Value
Second row: column-index for the specific row 7
Third row: column-index for an arbitrary row $k$

Zero-Boundary-Condition
As I stated above the matrix-size is $n=(m_x+2)(m_y+2)$, but for zero boundary conditions you can reduce the system to $\tilde{n}=m_xm_y$.
E.g. the equation for point 2 reads: $$φ_{3,0}=0$$

Thus this translates to a row: $$\left[\begin{matrix} 0& …& 0 & 1 & 0 & … &0\end{matrix}\right]=0$$

This row does not add any value to the matrix system and the row can be erased from the system. And as the zero boundary condition also introduces a $0$ in the row for node 7 above: $$\left[\begin{matrix} … & 0 &\color{red}{0} & 0&…&0& a & 2(a+b)&a & 0 & … & 0 &b & 0&… \end{matrix}\right]$$ You can also erase column 7 from the system.

I hope that this is can be understood.

Periodic Boundary Condition

You did not state how exactly your boundary condition looks like, but this could be something like: $$φ(x,0) = φ(x,1)$$ which can be formed to $$φ(x,0) - φ(x,1) = 0,$$ which again can be represented in the equation system.

• I will make some edits based on your help, and will see where that leads me. I have been banging my head against a wall on this one for a while. I expect that this is a nudge in the right direction (which is all the answer I was looking for anyway). – rocksNwaves Apr 17 '18 at 11:24
• Since a "b" coefficient comes with a shift in the $y$ direction, it doesn't seem intuitive to me that I should move a distance of $k\pm (m_x +2)$ along a row to place that coefficient. On writing it out, I see that for interior points it works out great. My intuition is that the boundary conditions take care of the rest (because if not, then what else?) but I am afraid that my understanding of the math behind that intuition is very lacking. – rocksNwaves Apr 17 '18 at 13:59
• I get it! You take the grid and "unwind it into a string of values, according to the order of the solution vector, then simply count! – rocksNwaves Apr 17 '18 at 15:54
• Yes, you basically answered your question already. In this answer to a different question, you can also find an example for the matrix, but on an equidistant grid. As you said, you need to describe the coordinates $(i,j)$ for each node in the row-index $k$ representing each row. The mapping would be simply $k= i*(m_x+2)+j$. – P. Siehr Apr 17 '18 at 16:07