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
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
. This norm (called the grid norm) is the 2-norm multiplied by
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.
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
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.
for use with SOR is calculated by finding the smaller root of the quadratic equation in
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
. Substituting this in the expression for
. 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
; 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
is the negative of the strictly upper triangle of the
is the diagonal of
The following table gives the expression
for each method, based on the above values for
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.
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
 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.
 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.
 S. J. Farlow, Partial Differential Equations for Scientists and Engineers,
New York: Dover, 1993.