We are interested in solving the Poisson equation, which can be expressed
in a domain \(\Omega\subset \mathbb{R}^{n}\) with a boundary \(\partial\Omega\) as
where \(a\) is a known constant, \(\nabla^{2}=\nabla\cdot\nabla(\cdot)\) is the Laplacian operator, \(\phi=\phi(\vec{x})\)
is our scalar variable that is a function of space, and \(f=f(\vec{x})\) is a forcing function. The second line represents a Dirichlet
condition on the boundary \(\partial\Omega\), where \(g=g(\vec{x})\) is a prescribed function along the boundary.
Let \(\Omega\) represent a 2D square domain given by \([0, 1] \times [0, 1]\). We choose to discretize the above equation using a 2nd-order
finite difference approximation on a cartesian mesh composed of a number \(N\) nodes in the x- and y-directions
with uniform spacing \(h\). The value of \(\psi\) on the 2D cartesian mesh can then be approximated
for each node \(\{i, j\}\) in the interior of the computational domain as
where the subscripts \(i\) and \(j\) represent the indices of the current node in the computational domain for the
\(x\)- and \(y\) -directions, respectively. Assuming a Dirichlet boundary condition specifies the values for \(\phi\) on the
boundary of the domain, we can consider the values for \(\phi_{1,j}\), \(\phi_{N,j}\) , \(\phi_{i,1}\), and \(\phi_{i,N}\)
known. The result is a system of \((N − 2) \times (N − 2)\) linear equations for the unknown
values of \(\phi_{i,j}\) in the interior of the domain.
2D uniform mesh featuring a discretization with a 5-point stencil.
We have freedom in deciding how to organize and ultimately solve the system of linear equations. Given a natural ordering of the unknowns and the 5-point stencil, we can express the full
linear system:
\[A\Phi=\mathbf{f}\]
where the matrix \(A\) is banded as a result of the structured mesh discretization and is given by
At this point, a number of methods can be employed for solving the linear system in Eqn. 3, but we will
focus on the classical iterative methods in this project. The Gauss-Seidel will serve as an example in the
discussion below. Starting with a linear system
Geometric Multigrid
For typical iterative numerical solution methods, high-frequency (local) errors in the solution are well-damped, while lower frequency (global) errors are poorly damped. Therefore, the low-frequency errors are
difficult to eliminate, which leads to slower solver convergence, especially on fine meshes. The key idea behind
multigrid is that effective rates of convergence at all scales can be maintained in a solver by leveraging a
sequence of grids at various resolutions. With geometric multigrid, multiple levels of physical grids with
varying resolution are used to provide better approximations of the solution with each step of an iterative
solution method (i.e., a multigrid cycle).
To illustrate the basic components of linear multigrid for elliptic problems, define the error in the
solution to be the difference between the solution \(\Phi\) and the approximation to the solution \(\tilde{\Phi}\), or
\[\mathbf{e}={\Phi} -\tilde{\Phi}\]
where \(e\) is the error vector (one value per node in the computational mesh). We can also define a residual
vector \(r\), which is a measure of how well the discretized governing equations are being satisfied by our
numerical solution procedure, as
which relates the error in the solution to the residual. Eqn. \(A\mathbf{e}=\mathbf{r}\) allows us to compute a measure of the error
on coarser mesh levels after transferring the values of the residual from the fine mesh level onto the coarse
level (restriction). After calculating e on a coarse level, we can form a correction to the solution on the fine
mesh as
\[{\Phi}=\tilde{\Phi}+\mathbf{e}\]
upon transferring the error up to a fine mesh from the coarse mesh level below (prolongation). Furthermore,
we can apply these ideas recursively over an entire set of grids of various resolutions to complete a full
multigrid cycle, such as the V-cycle detailed in Alg. 1.
During a multigrid V-cycle, the solution is first approximated using several smoothing iterations with a
method like Gauss-Seidel on the finest mesh (pre-smoothing), and then the residual is transferred to the first
coarse level, where additional smoothing iterations occur. This restriction followed by smoothing continues
recursively until the coarsest mesh level is reached (the downstroke of the cycle). After performing some
smoothing iterations on the coarsest level, a correction for the solution values is transferred to the finer mesh
level above. This upward stroke of the cycle with prolongation and smoothing continues until a correction
is finally applied to the solution on the finest mesh. Typically, several final smoothing iterations (post-smoothing) are performed on the finest mesh before moving on to the next multigrid cycle. The downstroke
and upstroke of the cycle form a V-shape when viewed graphically, as in Fig. 3. Other cycles are possible,
and W-cycles are common, for instance.
For our discretized Poisson problem, we can express the value of the residual \(\mathbf{r}\) from Eqn. \(\mathbf{r}=\mathbf{f}-A\tilde{\Phi}\) at each node as
After computing \(\mathbf{r}\), we restrict the values down to the next coarse level to form the right-hand side of
Eqn. \(A\mathbf{e}=\mathbf{r}\). For a weighted restriction, we will include information from all of the fine nodes that surround a
particular coarse mesh node. To accomplish this, we will set the residual at a coarse mesh node to be the
sum of a contribution from the coincident node (1/4 of the value), the nodes that are part of the stencil in
the north, south, east, and west directions (1/8 of the value), and the diagonal neighbors, i.e., north-east,
north-west, south-east, and south-west (1/16 of the value).
A prolongation operation is one that transfers the correction from a coarse mesh to a fine mesh. Similar to
the weighted method for restriction, we will perform a weighted prolongation by setting the correction at a
fine mesh node to be the value of the correction at a coincident coarse node, if applicable, or as the sum of
a contribution from the coarse nodes that are nearest in the north, south, east, and west directions (1/2 of
the value) and the nearest diagonal neighbors, i.e., north-east, north-west, south-east, and south-west (1/4
of the value).
The second class of intergrid transfer operations involves moving vectors from
a fine grid to a coarse grid. They are generally called restriction operators and are
denoted by \(I_{h}^{2h}\). The most obvious restriction operator is injection. It is defined by
\(I_{h}^{2h}\mathbf{v}^{h}=\mathbf{v}^{2h}\), where
In other words, with injection, the coarse-grid vector simply takes its value directly
from the corresponding fine-grid point.
An alternate restriction operator, called full weighting, is defined by \(I_{h}^{2h}\mathbf{v}^{h}=\mathbf{v}^{2h}\), where
The fact that the interpolation operator and the full weighting operator are transposes of each other up to a constant is called a variational property and will soon
be of importance.
For the sake of completeness, we give the full weighting operator in two dimensions. It is just an averaging of the fine-grid nearest neighbors. Letting \(I_{h}^{2h}\mathbf{v}^{h}=\mathbf{v}^{2h}\),
we have that
In our discussion of intergrid transfers, we consider only the case in which the
coarse grid has twice the grid spacing of the next finest grid. This is a nearly
universal practice, because there is usually no advantage in using grid spacings
with ratios other than 2. Think for a moment about the step in the correction
scheme that requires transferring the error approximation \(\mathbf{e}^{2h}\) from the coarse grid
\(\Omega^{2h}\) to the fine grid \(\Omega^{h}\). This is a common procedure in numerical analysis and is
generally called interpolation or prolongation. Many interpolation methods could
be used. Fortunately, for most multigrid purposes, the simplest of these is quite
effective. For this reason, we consider only linear interpolation.
The linear interpolation operator will be denoted \(I_{2h}^{h}\). It takes coarse-grid vectors and produces fine-grid vectors according to the rule \(I_{2h}^{h}\mathbf{v}^{2h}=\mathbf{v}^{h}\), where
For two-dimensional problems, the interpolation operator may be defined in a
similar way. If we let \(I_{2h}^{h}\mathbf{v}^{2h}=\mathbf{v}^{h}\), then the components of \(\mathbf{v}^{h}\) are given by
