Solving the 2D Poisson PDE by Eight Different Methods

This Demonstration considers solutions of the Poisson elliptic partial differential equation (PDE) on a rectangular grid. Eight numerical methods are based on either Neumann or Dirichlet boundary conditions and nonuniform grid spacing in the and directions. Different source functions are considered.


  • [Snapshot]
  • [Snapshot]
  • [Snapshot]


The 2D PDE solved is . You can specify either Neumann or Dirichlet boundary conditions. At least one of the edges must have a Dirichlet boundary condition. You can choose the initial guess of the solution to be zero or a random number. Iteration toward a solution stops when , where is the tolerance you specify using the slider located at the bottom of the solver panel and is the source vector. is the 2-norm of the residual grid at the end of step , where . This norm (called the grid norm) is the 2-norm multiplied by , where is the grid spacing (for the case of uniform grid spacing). You can see a plot of as a function of the step number by selecting the convergence option from the plot pull-down menu in the top row.
Using superscript as the step number, the following table summarizes the finite difference scheme for Jacobi, Gauss–Seidel, and successive overrelaxation (SOR) for the case of uniform grid spacing.
The following table summarizes the finite difference scheme for nonuniform grid spacing.
The SOR method is the same as Gauss–Seidel except for the use of the parameter to accelerate the convergence. The SOR with Chebyshev acceleration is the same as Gauss–Seidel red/black except that the parameter is updated at each step according to the following scheme: initially it is set to 1, then after the first sweep over the red squares , and from then on , where is the Jacobi spectral radius.
The Gauss–Seidel red/black is a variation of the Gauss–Seidel, where the Gauss-Seidel scheme is first applied to the even grid points and then applied to the odd points. A grid point is considered even when is even and odd otherwise.
In the above table, is the SOR parameter. You can choose to enter a numerical value for this parameter using the input labeled "user" or select "optimal" so that the program calculates the optimal value. For SOR with Chebyshev acceleration, the optimal from the UI is not used, but is calculated at each step as described above.
The optimal for use with SOR is calculated by finding the smaller root of the quadratic equation in given by , where and is the grid spacing. This gives . This can also be written as . This formula is known to be optimal for uniform grid spacing and equal sides and with homogeneous Dirichlet or Neumann boundary conditions.
Performance of the Numerical Methods
The following table summarizes the performance of the numerical schemes as a function of the number of grid points on a small test problem. The force function used to generate this table is with convergence tolerance set to . The initial guess is zero with zero Dirichlet boundary conditions on the four edges of the square with centered grid and a square of side length one. The numbers inside the parentheses are the number of steps to converge and the CPU time used by the solver. CPU time shown does not include the time used for plotting.
The above table shows that among the stationary methods, SOR with Chebyshev acceleration used the smallest number of steps to converge.
Neumann boundary conditions: For an edge with Neumann boundary conditions, the method of ghost points was used to derive the finite difference scheme. To illustrate this method and assuming the west edge has Neumann boundary condition defined by , a ghost point is added to the left of the edge at location , and the first derivative at the west edge point is approximated using a centered difference formula, giving . In this formula, the first column in the grid has an index of . The five-point Laplacian is applied for (a point on the row at the west edge), giving . The above equation for gives . Substituting this in the expression for gives . This formula is used to evaluate on the west edge. Corner grid points are approximated by the average of the adjacent grid points on the two edges meeting at the corner. Similar expressions are used to update the other three edges for Neumann boundary conditions. The table below summarizes the expressions for the finite difference scheme at each edge for the case of Neumann boundary conditions with uniform grid spacing.
Matrix splitting method formulas: Only Dirichlet boundary conditions are supported for the matrix splitting method for Jacobi, Gauss–Seidel, and SOR. Let and ; then the iterative solution becomes . The following table gives the formulas for the splitting method for Jacobi, Gauss–Seidel, and SOR. In this table, is the negative of the strictly lower triangle of the matrix, is the negative of the strictly upper triangle of the matrix, and is the diagonal of .
The following table gives the expression for each method, based on the above values for and .
Some Observations and Possible Issues
Error reduction can be seen to be most rapid in the first few steps of iterations. This makes some of the methods suitable to be used as smothers for the multigrid method.
The formula and the Chebyshev acceleration scheme used are optimal only for homogeneous Neumann and homogeneous Dirichlet boundary conditions. For non-homogeneous boundary conditions, the formula used might not produce the fastest convergence. You can choose to enter a different value.
[1] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, Philadelphia, PA: Society for Industrial and Applied Mathematics, 2007.
[2] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in Fortran 77: The Art of Scientific Computing, 2nd ed., New York: Cambridge University Press, 1992.
[3] S. J. Farlow, Partial Differential Equations for Scientists and Engineers, New York: Dover, 1993.
    • Share:

Embed Interactive Demonstration New!

Just copy and paste this snippet of JavaScript code into your website or blog to put the live Demonstration on your site. More details »

Files require Wolfram CDF Player or Mathematica.