A Cascadic Multigrid Algorithm on the Shishkin Mesh for a Singularly Perturbed Elliptic Problem with regular layers

A two-dimensional linear elliptic equation with regular boundary layers is considered in the unit square. It is solved by using an upwind difference scheme on the Shishkin mesh which converges uniformly with respect to a small perturbation parameter. The scheme is resolved based on an iterative method. It is known that the application of multigrid methods leads to essential reduction of arithmetical operations amount. Earlier we investigated the cascadic two-grid method with the application of Richardson extrapolation to increase accuracy of the difference scheme by an order uniform with respect to a perturbation parameter, using an interpolation formula uniform with respect to a perturbation parameter. In this paper a cascadic multigrid algorithm of the same structure is studied. We construct an extrapolation of initial guess using numerical solutions on two coarse meshes to reduce the arithmetical operations amount. The application of the Richardson extrapolation method based on numerical solutions on the last three meshes leads to increase accuracy of the difference scheme by two orders uniformly with respect to a perturbation parameter. We compare the proposed cascadic multigrid method with a multigrid method with V-cycle. The results of some numerical experiments are discussed.


Introduction
The two-dimensional linear elliptic problem with regular boundary layers in the unit square is considered. It is well known that the application of classical difference schemes for a singularly perturbed problem leads to large errors for small values of the perturbation parameter [1][2][3][4][5][6][7][8][9][10]. The uniform convergence of a difference scheme for such problem can be provided by fitting the scheme to a boundary layer component [1,2] or by using a mesh which is dense in a boundary layer [3,5].
The multigrid method [11][12][13][14][15] and the two-grid method [16,17] as a special case of multigrid method lead to essential reduction of the number of arithmetical operations. These methods can be effective applied for singularly perturbed problems, see [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] and the references therein. To increase the accuracy of the difference scheme in multigrid method the Richardson extrapolation [7,34] can be applied. The two-grid method with the application Richardson extrapolation to increase the ε-uniform accuracy of the difference scheme on the Shishkin mesh is investigated in [25][26][27][28][29][30]33]. In [29,30,33] the multigrid algorithm of the same structure based on three meshes is considered and to reduce the number of iteration we apply the idea of the extrapolation of numerical solutions on coarse meshes [35,36].  [37][38][39] and the references therein. It is shown numerically that the cascadic two-level method for solving semilinear elliptic boundary value problem base on radial basis function-generated finite difference method is computationally effective and has simplicity of programm, may be apply in the cases of the square and L-shape domains for 2D and 3D [38]. It is shown numerically that the extrapolation cascadic multigrid (EXCMG) method combined with the fourth-order compact difference scheme to solve the Dirichlet boundary value problem of the 3D Poisson equation in rectangular domains is much more efficient comparing to the classical V-cycle and W-cycle multigrid method and it is particularly suitable for solving large scale problems [37]. It is shown numerically that the EXFMG method with a fourth-order compact difference scheme, a Lagrange interpolation and a simple extrapolation technique to solve the 2D convectiondiffusion equation improves the solution accuracy and keeps less computational costs, compared to the classical multigrid methods and can also be extended to solve other kinds of partial differential equations [39].
For reaction-diffusion problems a boundary-layer preconditioning approach is proposed and analyzed in [22]. It is shown in [18] that standard multigrid methods are not convergent for singularly perturbed elliptic problem and a multigrid method with V-cycle and a modified restriction operator is proposed. The modified multigrid algorithm has similar properties that standard multigrid for classical elliptic problems.
The aim of this work is to study a cascadic multigrid method using Richardson extrapolation for a singularly perturbed elliptic problem based on the ε-uniform difference scheme on the Shishkin mesh and to compare with V-cycle multigrid method form [18].
Notation: Here C, sometimes subscripted, denotes a generic positive constant that is independent of the perturbation parameter ε and the step size of the mesh.

Preliminaries
We consider the following boundary value problem: where the coefficients a, b, c are bounded and satisfy the condition: the perturbation parameter ε takes arbitrary values in the open-closed interval (0, 1]. The coefficients, the right-hand side f and the boundary function g are sufficiently smooth. We also assume that the sufficient compatibility conditions are satisfied. It is known [5][6][7] that the solution of problem (1) under condition (2) is uniformly bounded and has two regular boundary layers near x = 0 and y = 0.
Define the piecewise-uniform mesh [5] in the domain Ω: where α and β are given in (2).  3 We consider the upwind difference scheme for the problem (1) on the mesh Ω N : Let [u] Ω be the projection of a function u(x, y) on a mesh Ω. According to [5,6] the difference scheme (4) on the mesh (3) converges ε-uniformly and the following accuracy estimate is satisfied: where C is positive constant that is independent of the perturbation parameter ε and the parameters of the mesh (3).
Notice that the matrix of this system (4) is M -matrix and a number of iterative methods for its resolving converge [40,41].

Cascadic multigrid method
We study a cascadic multigrid method with structure like in the papers [25,[27][28][29][30]33], and with the number of nodes of the coarsest mesh 4. To improve the accuracy of the difference scheme (4) in the multigrid method we apply Richardson extrapolation [7,34] using three last numerical solutions.
Let u N be the solution of the difference scheme (4) on the mesh Ω N . For the application of Richardson extrapolation we use the solutions of the difference scheme (4) on the meshes Ω N/2 and Ω N/4 which ones have the same value of the parameters σ x and σ y as the refined mesh Ω N . Thus these meshes are nested that is Let us define the solution of the difference scheme (4) on the mesh Ω N/2 as u N/2 and on the mesh Ω N/4 as u N/4 . According to the Richardson extrapolation we introduce the function u nN on the mesh Ω N . At first we define the function u nN at the nodes of the coarse mesh (X l , Y m ) ∈ Ω N/4 the formula: Now let us specify u nN (x i , y j ) at the nodes of refined mesh Ω N \ Ω N/4 using the interpolation function.
To construct a better initial guess on refined mesh Ω N we use the same idea of articles [35,36] and the proposed formulas in [29]: where We can get the values at other nodes on the mesh Ω N \ Ω N/2 , using an appropriate interpolation.
It is stated in the following algorithm. Step We proposed to use a residual stop criterion for the coarsest mesh k = 1 and a corrected residual stop criterion for other meshes k > 1 Note that we can use the stop criterion (10) only since we have sufficiently good guess. In [18] the multigrid method with special restriction operator is proposed and it is shown that this algorithm has good properties multigrid method. A single iteration of the multigrid with V-cycle U (1) = M G(U (0) , F, N ).
Step 2.1. Apply relaxation to AU = F with initial guess U (0) , producing U (r) .
Step 2.4. Compute U (c) = U (r) + P U c for interpolation matrix P .
Step 2.5. Apply relaxation AU = F to with initial guess U (c) , producing U (r) .
Here the relaxation (the smoother) is a line Gauss-Seidel of alternating symmetric type, the interpolation (prolongation) operator P is the bilinear interpolation and the restriction operator R from [18].

Results of numerical experiments
We consider the following boundary layer problem: where f corresponds to the exact solution: The solution of the problem (11) is computed based on the difference scheme (4). We define the initial guess as The difference scheme (4) can be written as the five-point scheme: Then vector-matrix form of the Gauss-Seidel method can be presented as: Notice that the efficiency of Gauss-Seidel method in the case of the problem (1) depends on the ordering of equations and unknowns [40,41,47]. Therefore we study a line Gauss-Seidel of alternating symmetric type, the lines are processed in x-and y-directions in forward and backward lexicographical order. Table 1 contains the error norm for a multigrid method with V-cycle (left table) and for a cascadic multigrid method (right table) for various values of ε and N . Note that the error norm for a one-grid method and multigrid method with V-cycle is the same.  Table 2 contains the time of computing for a one-grid method for various values of N and ε.   Table 3 contains the time of computing for a multigrid method with V-cycle for various values of N and ε.
It follows from table 2, 3 that the time of computing for one-grid method depends from a small parameter ε and does not increment linearly and it does not similar to properties multigrid method with V-cycle. But for the size of mesh less than 128 × 128 one-grid method is slightly more computationally efficient. Table 4 contains the time of computing for a cascadic multigrid method for various values of N and ε. It follows from table 3, 4 that the time of computing for cascadic multigrid method is more computationally efficient and does not depend from a small parameter ε but it does not increment linearly too. Thus, we find that the numerical solutions from the cascadic multigrid method are much more accurate than those from the multigrid method with V-cycle. This may be considered how, to achieve the same accuracy, fewer mesh points are needed for the cascadic multigrid method than for the multigrid method. For instance, the numerical solution by the cascadic multigrid method on the 64 × 64 mesh is comparable to the numerical solution obtained by the multigrid method with V-cycle on the 512 × 512 mesh in Table 1.
We note that a cascadic multigrid method can have good properties similar to the multigrid method with V-cycle. To achieve it we can use instead of the residual stop criterions (9), (10) in Step 1.3 on each level (the coarse meshes) a fixed number of iterations, for instance, N k = 10 times. Table 5 contains the time of computing for a cascadic multigrid method with fixed number of iterations N k = 10 for various values of N and ε.  Here it should be noted that such a stop criterion leads to a loss of the accuracy of the numerical solution in the Step 1.4 since the Richardson extrapolation does not work. Table 6 contains the error norm for a cascadic multigrid method with fixed number of iterations N k = 10 without the Richardson extrapolation (left table) and with the Richardson extrapolation (right table) for various values of N and ε. It follows from table 6 that the application of extrapolation (6) in the case of the use of a stop criterion with a fixed number of iteration is useless. However we have that the time of computing for the cascadic multigrid method is much better than those of the multigrid method with V-cycle but it should be remembered that the number of iteration N k = 10 is a numerical result.