Discrete Poisson Equation with Pure Neumann Boundary Conditions

Question

I'm trying to implement the Helmholtz-Hodge Decomposition in 2D, which states that a vector field is composed by a rotational free component, a divergence free component and a harmonic component.

This leads me to a Poisson Equation:

$\begin{equation} \left\{ \begin{array}{lcl} \nabla \cdot \mathbf{F} & = & \Delta \varphi \\ (\nabla \cdot J) \mathbf{F} & = & - \Delta \psi \end{array} \right. \end{equation}$

I'm trying to solve it using finite-differences.

To comput gradients, I'm using:

$\left( \begin{array}{ccc} -1 & 0 & 1 \end{array} \right)$

For dx.

$\left( \begin{array}{c} -1 \\ 0 \\ 1 \end{array} \right)$

For dy.

For the Laplacian, I'm using:

$\left( \begin{array}{ccc} 0 & -1 & 0 \\ -1 & +4 & -1 \\ 0 & -1 & 0 \end{array} \right)$

Then, for a 2 x 2 input, for example, I will have an equation of the form:

$\left( \begin{array}{c} \nabla \cdot \mathbf{F}_{11}\\ \nabla \cdot \mathbf{F}_{12}\\ \nabla \cdot \mathbf{F}_{21}\\ \nabla \cdot \mathbf{F}_{22}\\ \end{array} \right) $= $\left( \begin{array}{cccc} 4. & -1. & 0. & -1.\\ -1. & 4. & -1. & 0.\\ 0. & -1. & 4. & 0.\\ -1. & 0. & 0. & 4 \end{array} \right) $ $\left( \begin{array}{c} \varphi_{11}\\ \varphi_{12}\\ \varphi_{21}\\ \varphi_{22} \end{array} \right) $

To solve $\nabla \cdot \mathbf{F} = \Delta \varphi$.

Using this approach, I could solve the system and obtain the vector field components. However, I'm having problem at the boundaries.

Researching, I read about the Pure Neumann Boundary Conditions. However, I don't know how to include them in the system of equations. I'm wondering how to do that.

Thank you.

Answer 1

Previous comments gave you good suggestions, I try to add some more.

Firstly, your example 2x2 really does not correspond to zero Neumann boundary conditions. In fact, one can show that your choice leads to zero Dirichlet boundary conditions.

If you domain is a rectangle, then implementing the zero Neumann condition can be done easily. You have to include the values at boundaries to unknowns (you might done that already). For four nodes in the corners of your rectangle, you have to change in the corresponding row of matrix two -1 to -2, for all other rows corresponding to boundary nodes change only one -1 to -2, it must be exactly the single entry that is inside of rectangle so not on boundary. As noted by someone else previously, the sum off all entries for each row in the matrix must be equal to zero, that is another check for you if you are doing it correctly.

The idea behind finite difference discretization for zero Neumann boundary conditions is that you imagine that you have one row or one column of nodes next to each side of your rectangle with additional (artificial) unknowns. It means that you may use the standard stencil with 4 on the diagonal and four time -1 offdiagoals. Additional unknowns can be eliminated by adding discrete equations that correspond to zero Neumann boundary condition, e.g. $\partial_x \varphi = 0$, that takes in the discrete form (hopefully the notation is clear) e.g. $(\varphi_{1 j} - \varphi_{-1 j})/dx = 0$, so the artificial unknowns are equal $\varphi_{-1 j} = \varphi_{1 j}$. Plug this into the discrete equation for unknowns $\varphi_{0 j}$ and you get what I did suggest. By the way this corresponds to second order accurate approximation that you should prefer.

Another point is that applying zero Neumann boundary condition for Poisson equations gives you underdetermined system where there are formally infinitely many solutions, but they differ only by constant (it means if you subtract two solution functions you get constant function). Once you fix one value e.g. to 0 (as proposed by someone else), you choose one of those infinitely many solution, so your problem has then unique solution. It is like you define Dirichlet boundary condition (e.g. $\varphi =0$) in one single point at boundary. Of course, then such node must be excluded from your unknowns, so you do not have a row in your matrix for this node.

Sometimes your linear solver can take care for this, but this is another story, the suggestion above may be much simpler.

Answer 2

Pure Neumann boundary conditions give a non-unique solution. You will need to set some other condition to make it unique. The simplest would be to choose some solution value somewhere and set it to any constant, I would choose a corner or the middle and set it to 0. Another common approach is to apply a zero mean condition that the integral of the solution of the entire solution is zero.

Answer 3

Your difference stencil is wrong for nodes at the boundary (which in your 2x2 case happens to be every row). A simple way to see this is to realize that if you applied the Laplace operator to a constant function, the result is zero. Similarly, if you applied the Laplace matrix to a vector where all components have the same value (the discretized version of a constant function) then you need to get the zero vector. This implies that the sum of elements in each row must be zero -- which it is for the 5-point stencil you are using at interior nodes, but it isn't the way you wrote it down. You just need to do something different in those rows that correspond to boundary nodes.