Multigrid Method

We usually use iterative methods for solving linear system when the matrix arising from problem is sparse. A clear advantage of using iterative methods is that they require far less computational effort. The principle of iterative method is beginning with an initial guess, and improve approximation successively until it is as accurate as desired. Most of iterative methods can reduce efficiently oscillatory components of the error, but it is much less effective as smooth components remained (sometimes relaxation is also called smoother).

The notion of smooth and oscillatory components are relative to the grid size which the solution is defined. The smooth error on the fine grid is more oscillatory when projected on the coarse grid. Therefore, we might relax on the fine grid to reduce oscillatory components of the error and then relax on the coarse grid when smooth components remained.

What is the relationship between the error and the approximation? An important scheme: residual correction gives a appropriate answer. Suppose that the system Au = f has a unique solution and that v is a computed approximation to u. It is easy to compute the residual r = f − Av. The error e = u − v satisfies the residual equation: Ae = r.

To improve the approximation v, we might solve the residual equation for e, and then compute a new approximation using the definition of the error u = v + e.

Now, we can combine these two ideas to produce two-grid correction scheme, as out-lined below [1].

• Interpolate the coarse-grid error to the fine grid by e^h = I_2h^h e^2h and correct the fine-grid approximation by v^h ←− v^h+ e^h.

• Relax ν₂ times on A^hu^h = f^h on Ω^h with initial guess v^h.

Ω^2h denotes coarse-grid has twice the grid spacing of the fine grid Ω^h. The procedure transferring the vectors from a fine-grid to a coarse-grid is called restriction, and this operator is denoted by I_h^2h; the procedure transferring the vectors from a coarser-grid to a fine-grid is called interpolation or prolongation, and this operator is denoted by I_2h^h . The arrow notation stands for replacement or overwriting. The integers ν₁ and ν₂ are parameters in the scheme that control the number of relaxation sweeps before and after

2

visiting the coarse grid. They are usually fixed at the start, based on either theoretical considerations or on past experimental results, and they are often small.

In general, injection and full weighting operators are more common restriction we used;

linear and cubic interpolations are popular methods. The issue of intergrid transfers is discussed at some books about multigrid method [1][2]. And some basic iterative method like Jacobi, Gauss-Seidel (GS) or red-black Gauss-Seidel (RBGS) can be found in elementary numerical books. Therefore, we do not mention all of them here. We only give definition of full weighting restriction, linear interpolation and red-black Gauss-Seidel method for two-dimensional problem since we will use these tools in the later numerical experiment.

In order to illustrate them more conveniently, we consider a simple second-order bound-ary value problem

½ −∆u = f (x, y), 0 < x < 1, 0 < y < 1

u = 0, on the boundary

on rectangle. The domain of the problem {(x, y) : 0 ≤ x ≤ 1, 0 ≤ y ≤ 1} is partitioned into n × n subrectangle by introducing the grid points x_i = ih, y_j = jh, where h = 1/n is the constant width of the each edge of the subrectangle. We also introduce v_ij as an approximation to the exact solution u(xi, yj) and introduce fij as the value of f (x, y) at (x_i, y_j).

•

Linear interpolation

If we let I_2h^h v^2h = v^h, then the components of v^h are given by v^h_2i,2j = v_ij^2h,

v^h_2i+1,2j = 1

2(v_ij^2h+ v_i+1,j^2h ), v^h_2i,2j+1= 1

2(v_ij^2h+ v_i,j+1^2h ), v_2i+1,2j+1^h = 1

4(v_ij^2h+ v_i+1,j^2h + v_i,j+1^2h + v_i+1,j+1^2h ), 0 ≤ i, j ≤ n 2 − 1.

Linear interpolation is effective when the vector is smooth.

•

Full weighting restriction

If we let I_h^2hv^h = v^2h, then the components of v^2h are given by v_ij^2h= 1

16[v_{2i−1,2j−1}^h + v_2i−1,2j+1^h + v^h_2i+1,2j−1+ v_2i+1,2j+1^h + 2(v^h_2i,2j−1+ v_2i,2j+1^h + v^h_2i−1,2j+ v_2i+1,2j^h ) + 4v^h_2i,2j], 1 ≤ i, j ≤ n

2 − 1.

The values of the coarse-grid vector are weighted averages of values at neighboring fine-grid points.

•

Red-black Gauss-Seidel relaxation

This method may be expressed in component form as below procedure. First, update red points v_ij by

v_ij ←− 1

4(v_i−1,j + v_i+1,j+ v_i,j−1+ v_i,j+1+ h²f_ij), where sum of index i and j is even. Then update black points v_ij by

v_ij ←− 1

4(v_i−1,j + v_i+1,j+ v_i,j−1+ v_i,j+1+ h²f_ij),

where sum of index i and j is odd. Fig. 1 shows red and black points. Red-black Gauss-Seidel has a clear advantage in terms of parallel computation.

4

Figure 1: A one-dimensional grid (top) and a two-dimensional grid (bottom) showing the red points ◦ and the black points • for red-black relaxation.

Two-grid correction is the basis of multigrid. Of course, we can repeat this process on successively coarser grids until arriving the coarsest grid. This is so-called V-cycle. To obtain more benefit from the coarser grids, where computations are cheaper, the W-cycle zigzags among the lower-level grids before moving back up to the finest grid. We can also join nested iteration idea which uses V-cycle on coarse grids to obtain improved initial guesses for fine grids. This is so-called full multigrid V-cycle (FMG).

The schedule of grids for V-cycle, W-cycle, and FMG show in Fig. 2. The complete algorithms of them are given in Algorithm 1,2,3 respectively. To describe Algorithm conveniently, we change some notations. We call the right-side vector of the residual equation f, rather than r, because it is just another right-side vector. Instead of calling the solution of the residual equation e, we use u because it is just a solution vector. We can then use v to denote approximations to u.

Algorithm 1 V^h(v^h, f^h) V-cycle Method if (h=coarsest) then

v^h = (A^h)⁻¹f^h {solve A^hu^h = f^h directly}

else

v^h=Relax(v^h, f^h, ν₁) {pre-relax ν₁ times with initial v^h } f^2h= I_h^2h(f^h− A^hv^h) {compute f^2h by restricting residual}

v^2h= 0

v^2h= V^2h(v^2h, f^2h) {V-cycle on the coarser grid } v^h = v^h+ I_2h^h v^2h {correct v^h by interpolating v^2h} v^h=Relax(v^h, f^h, ν₂) {post-relax ν₂ times with modified v^h } end if

Algorithm 2 W^h(v^h, f^h) W-cycle Method if (h=coarsest) then

v^h = (A^h)⁻¹f^h {solve A^hu^h = f^h directly}

else

v^h=Relax(v^h, f^h, ν₁) {pre-relax ν₁ times with initial v^h } f^2h = I_h^2h(f^h− A^hv^h) {compute f^2h by restricting residual}

v^2h = 0

v^2h = V^2h(v^2h, f^2h) 2 times {V-cycle on the coarser grid 2 times}

v^h = v^h+ I_2h^h v^2h {correct v^h by interpolating v^2h} v^h=Relax(v^h, f^h, ν₂) {post-relax ν₂ times with modified v^h } end if

Algorithm 3 F MG^h(f^h) Full Multigrid V-cycle Method if (h=coarsest) then

v^h = (A^h)⁻¹f^h {solve A^hu^h = f^h directly}

else

f^2h= I_h^2h(f^h) {compute f^2h by restricting f^h} v^2h= F MG^2h(f^2h) {FMG on the coarser grid }

v^h = I_2h^h v^2h {compute v^h by interpolating v^2h} v^h = V^h(v^h, f^h) {V-cycle with modified v^h }

end if

In this section, we have studied the elements of multigrid method. Although we don’t introduce all the details of multigrid method, we have enough ability to solve Helmholtz equations. [1] has complete introduction of multigrid method.

In brief, multigrid method integrate some easy ideas and schemes. In fact, these ideas may have some individual defects. Multigrid method arrange them skillfully such that they can work together, and result in a very useful numerical method.

6

Figure 2: Schedule of grids for V-cycle, W-cycle and FMG method from top to bottom.