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:
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
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.