I am using the Finite Difference Method to solve Poisson's equation

$$\frac{\partial \phi}{\partial z^2} = \frac{\rho}{\epsilon}$$

To do it is discretized according to the Finite Difference Approximation of the second order derivative yielding the following set of equations for each grid point

$$\frac{1}{2\Delta z^2}(\phi_{i+1}+\phi_{i-1}-2\phi_{i}) = \frac{\rho_{i}}{\epsilon} $$

My question is the following: Is there any simple way to see how many points and thus what spacing $\Delta z$ one should choose to produce reliable results? Obviously you should choose more than 5 but should I choose 500, 5000 or 50000? Asked another way: Over what length scale will the electrostatic potential vary significantly? I suppose if this is known one should make the spacing smaller than this length scale.

I once did a similar thing for the Schrödinger Equation. In this case it was straightforward to see what the characteristic wavelength of the solutions were and then to choose the spacing much smaller than this.

  • $\begingroup$ Besides convergence you'll need stability too, in many cases you'll end up with an integral with an exponential of a matrix so I think you want all the eigenvalues negative for that to not explode. Not 100% sure if it was the eigenvalues or determinant or something with the SVD but it might be good to remember if refining the mesh isn't enough and you wonder what could be the source of failure. $\endgroup$
    – Emil
    Feb 13, 2017 at 23:05

1 Answer 1


You're looking for the convergence of the simulation. For the simplest FD scheme, you have $$ u(x)\approx\frac{f(x+h) - f(x-h)}{2h} $$ which is $\mathcal{O}(h^2)$, as can be easily verifies with a Taylor expansion. Here, $h$ is going to be inversely proportional to the number of grid points, $N$: $$ h\propto\frac{1}{N} $$ where the equality holds if the width of the domain is 1 unit. So if you increase the number of grid points, you'll decrease $h$ and thus the error $h^2$ decreases even more.

To verify the convergence of your simulation, you should, ideally, compare your simulation to a known solution. If that does not work, then you can compare your solution at three resolutions (typically $N/2$, $N$ and $2N$). With these solutions, you would look at the L-2 norm: $$ E(N)=\frac{\vert u_{N} - u_{N/2}\vert}{\vert u_{2N} - u_{N}\vert} $$ where $$ \vert y\vert = \sqrt{ \frac{1}{N}\sum_{j=1}^ny\left(jh\right)^2} $$ is the $L_2$ norm. The value $E(N)/E(N/2)$ should be approximately 2. You'll probably want to pick $N$ to be 100, 200, 400, 800 and 1600 to start with and work from there (you may need to go higher).

As you do this, you'll note that $E(N)$ will reach some minimum and then start increasing due to rounding errors. It's going to be before this increase that you'll want to stop increasing $N$ and use a value in that range. You may also find that $E(N)$ is sufficiently small for your purposes and it's fine to stop there.

It is also possible, though not likely in 1D, that the number of grid points becomes too large to handle on your computer efficiently and you are not at an error level you are comfortable with. In this case, you should try a higher order scheme (Wikipedia has a nice table of coefficients for various finite difference schemes) and repeat the exercise.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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