0
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – rocksNwaves Apr 17 '18 at 11:24
  • $\begingroup$ 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. $\endgroup$ – rocksNwaves Apr 17 '18 at 13:59
  • $\begingroup$ 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! $\endgroup$ – rocksNwaves Apr 17 '18 at 15:54
  • $\begingroup$ 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$. $\endgroup$ – P. Siehr Apr 17 '18 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.