Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am interested in solving the Poisson equation using the finite-difference approach. I would like to better understand how to write the matrix equation with Neumann boundary conditions. Would someone review the following, is it correct?

The finite-difference matrix

The Poisson equation,

$$ \frac{\partial^2u(x)}{\partial x^2} = d(x) $$

can be approximated by a finite-difference matrix equation,

$$ \frac{1}{(\Delta x)^2} \textbf{M}\bullet \hat u = \hat d $$

where $\textbf{M}$ is an $n \times n$ matrix and $\hat u$ and $\hat d$ are $1 \times n$ (column) vectors,

Finite-difference matrix of the Poisson equation

Adding a Neumann boundary condition

A Neumann boundary condition enforces a know flux at the boundary (here we apply it at the left-hand side where the boundary is at $x=0$),

$$ \frac{\partial u(x=0)}{\partial x} = \sigma $$ writing this boundary condition as a centred finite-difference,

Error in equation. NB. I originally made an error here, sign error and didn't divide by 2. The following has been corrected. $$ \frac{u_2 - u_0}{2\Delta x} = \sigma $$

Note the introduction of a mesh point outside the original domain ($u_0$). This term can be eliminated by introducing the second equation, $$ \frac{u_0 - 2u_1 + u_2}{(\Delta x)^2} = d_1 $$

The equation arrises from having more information because of the introduction of the new mesh point. It allows us to write the double derivative of the $u_1$ as the boundary in terms of $u_0$ using a centred finite-difference.

The part I'm not sure about

Combining these two equations $u_0$ can be eliminated. To show the working, let's first re-arrange for the unknown,

$$ u_0 = -2\sigma\Delta x + u_2 \\ u_0 = (\Delta x)^2 d_1 + 2 u_1 - u_2 $$

Next they are set equal and rearranged into the form,

$$ \frac{u_2 - u_1}{(\Delta x)^2} = \frac{d_1}{2} + \frac{\sigma}{\Delta x} $$

I chose this form because it is the same form as the matrix equation above. Notice that the $u$ terms are divide by $(\Delta x)^2$ both here and in the original equation. Is this the correct approach?

Finally, using this equation as the first row of the matrix,

Poisson equation with a Neumann boundary condition on the left-hand side (corrected)

Some final thoughts,

  1. Is this final matrix correct?
  2. Could I have used a better approach?
  3. Is there a standard way of writing this matrix?
share|improve this question
There are two errors in your calculation: a centered finite difference has to be divided by $2 \Delta x$. Second, $u_0 = - \sigma \Delta x + u_2$ is also wrong, because the minus must be a plus. – vanCompute Feb 22 '13 at 10:02
This is worked out quite nicely in LeVeque's finite difference text, chapter 2. – David Ketcheson Feb 23 '13 at 11:42
Thank you for your helpful comments. – boyfarrell Feb 25 '13 at 2:28
These issues are also well explained in… – Evgeni Sergeev May 25 '13 at 4:40
could you please see this… – usumdelphini Aug 4 '14 at 10:50
up vote 5 down vote accepted

I think you are on the right way. If you correct your errors, it will look very similar to

share|improve this answer
Thank you for the comments and corrections. I have updated my question with corrections. I believe it is correct now. – boyfarrell Feb 25 '13 at 2:25

Great observation to see that $u_0$ can be eliminated.

Step back and think about the problem for a second. Specifying a Laplace equation fundamentally states that each point is the average of its neighbors. This is commonly visualized as a rubber sheet, and helps me to think about these things. (Poisson is similar w/ more or less stretchy points)

When you specify the value of the solution surface at the outermost edges you are "pinning" the sheet down in space at those points. When you specify the sheet by its derivative at the edges, there are any number of solutions that fulfill the equation that are translations the sheet in space whilst maintaining the same actual shape and thus derivatives.

In a practical sense however this can be troublesome. Matrices are ill-conditioned, and solvers act unpredictably. The most common thing I have seen done is simply "pinning" the solution to a fixed point by specifying $u_0 = 0$ or to a constant relevant to the solution space.

share|improve this answer
So generally the Poisson equation is solved with at least one Dirichlet boundary condition, so that a unique solution can be found? I guess it makes sense that the Neumann boundary conditions only make sense when source and sinks are included, otherwise there are an infinite number of solutions. However, if I take the diffusion equation instead, sometime Neumann boundary conditions are required for the correct physics (e.g. no flux of the quantity through a boundary when du/dx=0). This is what I'm really interested in. Is the above method the correct approach for apply Neumann BCs? – boyfarrell Feb 22 '13 at 6:32
You cannot apply Neumann BCs at all sides of you paper. If you do, you will not have a unique solution. It must be pinned at least at one side. – vanCompute Feb 22 '13 at 10:06
@meawoppl: How does one specify the fixed point while also doing a direct matrix solve? – jvriesem Jul 30 '15 at 22:58
Typically just assigning a point to a constant, by means of setting just one term in a row to 1, the rest zero, and a value on the RHS that corresponds to the solution plane you want to see. – meawoppl Jul 31 '15 at 21:59

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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