Multilevel Preconditioner with Stable Coarse Grid Corrections for the Helmholtz Equation

Let $\mathcal{T}_{h}$ be a conforming quasi-uniform triangulation of $\Omega$, and denote the collection of edges/faces by $\mathcal{E}_{h}$, while the set of interior edges/faces by $\mathcal{E}_{h}^{I}$ and the set of boundary edges/faces by $\mathcal{E}_{h}^{B}$. For any $T\in{\mathcal{T}}_{h}$, we define $h_{T}:={\rm diam}(T)$. Similarly, for $e\in{\mathcal{E}}_{h}$, define $h_{e}:={\rm diam}(e)$. Let $h:=\max_{T\in{\mathcal{T}}_{h}}h_{T}$. Throughout this paper we use the standard notations and definitions for Sobolev spaces (cf. [1]). In particular, $(\cdot,\cdot)_{Q}$ and $\left\langle\cdot,\cdot\right\rangle_{\Sigma}$ for $\Sigma\subset\partial Q$ denote the $L^{2}$-inner product on complex-valued $L^{2}(Q)$ and $L^{2}(\Sigma)$ spaces, respectively. Denote by $(\cdot,\cdot):=(\cdot,\cdot)_{\Omega}$ and $\left\langle\cdot,\cdot\right\rangle:=\left\langle\cdot,\cdot\right\rangle_{\partial\Omega}$.

Now we define the energy space $V:=H^{1}(\Omega)\cap\prod_{T\in{\mathcal{T}}_{h}}H^{2}(T)$. For any $v\in V$ and an interior edge/face $e=T_{1}\cap T_{2}$, where $T_{1}$ and $T_{2}$ are two distinct elements of ${\mathcal{T}}_{h}$ with respective outer normals $n_{1}$ and $n_{2}$, we introduce the jump $\Lbrack\nabla v\cdot n\Rbrack|_{e}=\nabla v\cdot n_{1}|_{T_{1}}+\nabla v\cdot n_{2}|_{T_{2}}$. Define the sesquilinear form $b_{h}(\cdot,\cdot)$ on $V\times V$ as follows:

$$b_{h}(u,v):=(\nabla u,\nabla v)+J(u,v),\qquad u,v\in V,$$

where

$\displaystyle J(u,v):=\sum_{e\in{\mathcal{E}}_{h}^{I}}{\bf i}\gamma_{e}h_{e}\left\langle\Lbrack\nabla u\cdot n\Rbrack,\Lbrack\nabla v\cdot n\Rbrack\right\rangle_{e},$

(2.1)

where ${\bf i}\gamma_{e}$ is a complex number with positive imaginary part. The terms in $J(u,v)$ are so-called penalty terms and ${\bf i}\gamma_{e}$ are penalty parameters (cf. [19, 21]).

Clearly, $J(u,v)=0$ if $u\in H^{2}(\Omega)$ and $v\in V$. Thus, if $u\in H^{2}(\Omega)$ is the solution of (1.1)-(1.2), then there holds

$$a_{h}(u,v):=b_{h}(u,v)-\kappa^{2}(u,v)+{\bf i}\kappa\langle u,v\rangle=(f,v)+\langle g,v\rangle,\qquad v\in V.$$

(2.2)

We define the CIP approximation space $V_{h}$ as the standard finite element space of order $p$, i.e.,

$$V_{h}:=\big{\{}v_{h}\in H^{1}(\Omega):v_{h}|_{T}\in\mathcal{P}_{p}(T),\,T\in{\mathcal{T}}_{h}\big{\}},$$

where $\mathcal{P}_{p}(T)$ is the space of polynomials of degree at most $p$ on $T$. The CIP finite element approximation is to find $u_{h}\in V_{h}$ such that

$$a_{h}(u_{h},v_{h})=(f,v_{h})+\langle g,v_{h}\rangle,\qquad v_{h}\in V_{h}.$$

(2.3)

It is clear that if the parameters $\gamma_{e}\equiv 0$, then the CIP-FEM reduces to the standard FEM. It has been proved that the CIP-FEM is stable for any $\kappa,h,p>0$ [19, 21]. The penalty parameters may be tuned to reduce the pollution errors. The numerical results in [19] show that using about the same total degrees of freedom (DOFs), the CIP-FEM yields the least phase error comparing to the standard FEM and IPDG method [12].

Let $\{\mathcal{T}_{l}\}_{l=0}^{L}$ be a shape regular family of nested geometrically conforming simplicial triangulations of the computational domain $\Omega$ obtained by successive quasi-uniform refinement of an intentionally chosen coarse grid $\mathcal{T}_{0}$. We denote by $V_{l}$ the CIP approximation space on ${\mathcal{T}}_{l}$. It is easy to see that the spaces $\{V_{l}\}_{l=0}^{L}$ are nested, i.e., $V_{0}\subset V_{1}\subset\cdots\subset V_{L}$. But the bilinear forms $\{a_{l}(\cdot,\cdot)\}_{l=0}^{L}$ defined as in (2.3) on each $V_{l}$ are nonnested. Thus in this work we consider the following multilevel methods for CIP-FEM discretizations, which has also been used for solving the linear systems from nonconforming P1 finite element approximations (cf. [17]).

For brevity, we denote by $a_{l}(\cdot,\cdot)$ the bilinear form $a_{h_{l}}(\cdot,\cdot)$ on $V_{l}$, where $h_{l}$ is the mesh size of ${\mathcal{T}}_{l}$. Define projections $P^{C}_{l}$, $Q_{l}$ : $V_{L}\rightarrow V_{l}$ according to

$$a_{l}(P^{C}_{l}v,w)=a_{L}(v,w),\quad(Q_{l}v,w)=(v,w),\quad\forall v\in V_{L},\,w\in V_{l}.$$

The existence and uniqueness of the discrete problem (2.3) imply the well-posedness of each $P^{C}_{l}$. For $0\leq l\leq L$, we also define $A^{C}_{l}:V_{l}\rightarrow V_{l}$ by means of

$\displaystyle(A^{C}_{l}v,w)=a_{l}(v,w),\quad\forall v,w\in V_{l}.$

(2.4)

Define $F_{l}\in V_{l}$ by

$$(F_{l},v)=(f,v)+\langle g,v\rangle\quad\forall v\in V_{l}.$$

Then the CIP-FEM on level $l$ (cf. (2.3)) can be rewritten as:

$$A^{C}_{l}u^{C}_{l}=F_{l}.$$

(2.5)

For the smoothing strategy, in fact, when $\kappa h_{l}/p$ is small enough, either of weighted Jacobi and Gauss-Seidel relaxation can be applied. Otherwise, GMRES relaxation can be used as a smoother and this will be explained in the following sections. To be precise, we describe the smoothing strategy as follows:

We remark that we choose $\alpha=0.5$ in this paper. This choice is motivated by the efficient relaxation of Jacobi and Gauss-Seidel smoothers, and it ensures that the amplification factor for the Jacobi and Gauss-Seidel smoothers will not become too large (cf. section 4.3 in the following and section 2.1 in [5]). When $l\in S_{L}$, we also define the operator $(R_{l}^{C})^{t}$ by

$\displaystyle\big{(}(R_{l}^{C})^{t}\bar{w},\bar{v}\big{)}=\big{(}R_{l}^{C}v,w\big{)},\quad\forall v,w\in V_{l}.$

(2.6)

We can now state the multilevel method for solving the CIP-FEM discretization system on level $L$ which is non-recursive version of multilevel method.

Here we note that the GMRES relaxation is not linear in the starting value, and the error operator of the above algorithm can not be written directly in a product form which holds only for the particular case $G_{L}=\emptyset$. When this particular case is considered reasonably, we set

$\displaystyle T_{l}:=\mu_{l}R_{l}^{C}A^{C}_{l}P^{C}_{l}\text{ and }T_{l}^{*}:=\mu_{l}(R_{l}^{C})^{t}A^{C}_{l}P^{C}_{l},\quad l=0,1,\ldots,L,\ G_{L}=\emptyset,$

Then the error operator of Algorithm 2.1 for the case $G_{L}=\emptyset$ can be derived as

$\displaystyle E_{M}^{*}E_{M},\text{ where }E_{M}:=(I-T_{L})\cdots(I-T_{1})(I-T_{0}),E_{M}^{*}:=(I-T_{0}^{*})(I-T_{1}^{*})\cdots(I-T_{L}^{*}),$

(2.7)

where $I$ is the identity operator in $V_{L}$.

The LFA is based on certain idealized assumptions and simplifications: the boundary conditions are neglected and the problem is considered on regular indefinite grids ${\mathcal{G}}_{h}=\{x:x=x_{j}=jh,j\in\mathbb{Z}\}$. Although the Robin boundary condition (1.2) and other absorbing boundary conditions are often applied in realistic Helmholtz problem, the neglect of boundary condition does usually not affect the validity of LFA (cf. [18]). The LFA in this section can be considered as simplification for the analysis of multilevel method for Helmholtz problem (1.1-1.2) in one dimension.

Let ${\boldsymbol{L}}_{h}$ be a discrete operator with a stencil representation ${\boldsymbol{L}}_{h}=[l_{j}]_{h},j\in\mathbb{Z}$. For any $u_{h}$ defined on ${\mathcal{G}}_{h}$ and a fixed point $x\in{\mathcal{G}}_{h}$, ${\boldsymbol{L}}_{h}u_{h}$ reads in stencil notation as

$\displaystyle{\boldsymbol{L}}_{h}u_{h}(x)=\sum_{j\in J}l_{j}u_{h}(x+jh),$

where $J\subset\mathbb{Z}$ is a certain finite index defining the so-called stencil. The basic quantities in the LFA are the Fourier modes $\varphi_{h}(\theta,x)=e^{{\bf i}\theta x/h}$ with $x\in{\mathcal{G}}_{h}$ and Fourier frequency $\theta\in\mathbb{R}$. In fact, the frequency $\theta$ can be restricted to the interval $(-\pi,\pi]$ as a fact that $\varphi_{h}(\theta+2\pi,x)=\varphi_{h}(\theta,x)$. It is easy to see that the Fourier modes are all the formal eigenfunctions of ${\boldsymbol{L}}_{h}$:

$\displaystyle{\boldsymbol{L}}_{h}\varphi_{h}(\theta,x)=\widetilde{{\boldsymbol{L}}}_{h}(\theta)\varphi_{h}(\theta,x),\qquad x\in{\mathcal{G}}_{h},\theta\in(-\pi,\pi],$

(4.1)

where the eigenvalues of ${\boldsymbol{L}}_{h}$ can be presented as $\widetilde{{\boldsymbol{L}}}_{h}(\theta)=\sum_{j\in J}l_{j}e^{{\bf i}\theta j}$, which is called the Fourier symbol of the operator ${\boldsymbol{L}}_{h}$. Given a so-called low frequency $\theta^{0}\in\Theta_{\rm low}=(-\pi/2,\pi/2]$, its complementary frequency $\theta^{1}$ is defined as

$\displaystyle\theta^{1}=\theta^{0}-{\rm sign}(\theta^{0})\pi.$

(4.2)

Interpreting the Fourier modes as coarse grid functions yields

$$\varphi_{h}(\theta^{0},x)=\varphi_{2h}(2\theta^{0},x)=\varphi_{2h}(2\theta^{1},x)=\varphi_{h}(\theta^{1},x),\qquad\theta^{0}\in\Theta_{\rm low},\ x\in{\mathcal{G}}_{2h}.$$

The Fourier modes $\varphi_{h}(\theta^{0},x)$ and $\varphi_{h}(\theta^{1},x)$ are called $2h$-harmonics. These Fourier modes coincide on the coarse grid with mesh size $H=2h$, and they can be represented by a single coarse grid mode $\varphi_{2h}(2\theta^{0},x)$. Hence, each low frequency mode is associated with a high frequency mode. For a given $\theta^{0}\in\Theta_{\rm low}$, define the two dimensional subspace of $2h$-harmonics by

$\displaystyle E_{2h}^{\theta^{0}}:={\rm span}\{\varphi_{h}(\theta^{0},x),\varphi_{h}(\theta^{1},x)\},$

(4.3)

where $\theta^{1}$ is defined as in (4.2). A crucial observation is that the space $E_{2h}^{\theta^{0}}$ is invariant under both smoothing operators and correction schemes for general cases by two-level method. The invariance property holds for many well-known smoothing methods (cf. [18]), such as Jacobi relaxation, lexicographical Gauss-Seidel relaxation, et al.

Let ${\boldsymbol{M}}_{h}$ be a discrete two-level operator. In the following, we will show that a block-diagonal representation for ${\boldsymbol{M}}_{h}$ consists of $2\times 2$ blocks $\widetilde{{\boldsymbol{M}}}_{h}(\theta)$ (cf. [18]), which denotes the representation of ${\boldsymbol{M}}_{h}$ on $E_{2h}^{\theta^{0}}$. Then the convergence factor of ${\boldsymbol{M}}_{h}$ by the LFA is defined as follows:

$$\rho({\boldsymbol{M}}_{h})={\rm sup}\{\rho(\widetilde{{\boldsymbol{M}}}_{h}(\theta)):\ \theta\in\Theta_{\rm low}\},$$

where $\rho(\widetilde{{\boldsymbol{M}}}_{h}(\theta))$ is the spectral radius of the matrix $\widetilde{{\boldsymbol{M}}}_{h}(\theta)$. We can refer to [18] for generalizations to $k$-level analysis.

Now we give the Fourier symbols of different operators in multilevel method for CIP-P1 discretization of one dimensional Helmholtz equation (1.1). We always assume $\gamma_{e}\equiv\gamma$ for some constant $\gamma$ in (2.1). Denote by $t=\kappa h$, $R=-1-{\bf i}4\gamma-t^{2}/6$, $S=1+{\bf i}3\gamma-t^{2}/3$. Since the boundary condition is neglected in the LFA, the stencil presentation of discretization operator ${\boldsymbol{A}}^{C}_{h}$ from (2.4) can be derived as (cf. [20])

$\displaystyle{\boldsymbol{A}}^{C}_{h}=\frac{1}{h}[{\bf i}\gamma\ \ R\ \ 2S\ \ R\ \ {\bf i}\gamma]_{h}.$

Combining the above expression and (4.1) yields the Fourier symbol of ${\boldsymbol{A}}^{C}_{h}$ as

$\displaystyle\widetilde{{\boldsymbol{A}}}^{C}_{h}(\theta)=\frac{1}{h}({\bf i}2\gamma\cos{2\theta}+2R\cos{\theta}+2S).$

(4.4)

Obviously, the Fourier symbol of standard FEM with piecewise P1 approximation (FEM-P1) is $\widetilde{{\boldsymbol{A}}}^{F}_{h}(\theta)=\widetilde{{\boldsymbol{A}}}^{C}_{h}(\theta)|_{\gamma=0}$.

For simplicity, we use standard weighted Jacobi ($\omega$-JAC) and lexicographical Gauss-Seidel (GS-LEX) relaxations as the smoothers in the LFA. It is easy to derive the weighted Jacobi relaxation matrix as ${\boldsymbol{S}}^{J}_{h}={\boldsymbol{I}}_{h}-\omega{\boldsymbol{D}}^{-1}_{h}{\boldsymbol{A}}^{C}_{h}$, where ${\boldsymbol{I}}_{h}$ is indentity matrix, ${\boldsymbol{D}}_{h}$ consists of the diagonal of ${\boldsymbol{A}}^{C}_{h}$ and $\omega$ is a weighted parameter. Due to the fact that ${\boldsymbol{D}}_{h}=\frac{2S}{h}{\boldsymbol{I}}_{h}$, one can easily deduce the Fourier symbol of weighted Jacobi relaxation as follows:

$\displaystyle\widetilde{{\boldsymbol{S}}}^{J}_{h}(\theta)=1-\frac{\omega}{S}({\bf i}\gamma\cos{2\theta}+R\cos{\theta}+S).$

(4.5)

The GS-LEX relaxation matrix is ${\boldsymbol{S}}^{GS}_{h}=({\boldsymbol{D}}_{h}-{\boldsymbol{L}}_{h})^{-1}{\boldsymbol{U}}_{h}$, where $-{\boldsymbol{L}}_{h}$ is the strictly lower triangular part of ${\boldsymbol{A}}^{C}_{h}$ and $-{\boldsymbol{U}}_{h}$ is the strictly upper triangular part of ${\boldsymbol{A}}^{C}_{h}$. The Fourier symbol of ${\boldsymbol{S}}^{GS}_{h}$ can also be directly derived that

$\displaystyle\widetilde{{\boldsymbol{S}}}^{GS}_{h}(\theta)=-\frac{Re^{{\bf i}\theta}+{\bf i}\gamma e^{{\bf i}2\theta}}{Re^{{\bf-i}\theta}+{\bf i}\gamma e^{{\bf-i}2\theta}+2S}.$

(4.6)

Note that for the restriction matrix ${\boldsymbol{I}}^{2h}_{h}=[r_{j}]^{2h}_{h}$ and $x\in{\mathcal{G}}_{2h}$, there holds

$$({\boldsymbol{I}}^{2h}_{h}\varphi_{h}(\theta^{\alpha},\cdot))(x)=\sum_{j\in J}r_{j}e^{{\bf i}j\theta^{\alpha}}\varphi_{h}(\theta^{\alpha},x)=\sum_{j\in J}r_{j}e^{{\bf i}j\theta^{\alpha}}\varphi_{2h}(2\theta^{0},x),\qquad\alpha=0,1.$$

By an analogous stencil argument, the stencil presentation of full weighting restriction matrix can be derived to be ${\boldsymbol{I}}^{2h}_{h}=[1/4,1/2,1/4]_{h}^{2h}$. Then one can obtain the Fourier symbol of ${\boldsymbol{I}}^{2h}_{h}$ is

$$\widetilde{{\boldsymbol{I}}}^{2h}_{h}(\theta)=\frac{1}{2}(1+\cos{\theta}).$$

For the linear prolongation matrix ${\boldsymbol{I}}^{h}_{2h}$, one can also derive its Fourier symbol as follows (cf. [18]):

$$\widetilde{{\boldsymbol{I}}}^{h}_{2h}(\theta)=\frac{1}{2}(1+\cos{\theta}).$$

Weighted Jacobi and lexicographical Gauss-Seidel relaxations are the general smoothing operators. It is well-known that such two smoothers are unstable especially for linear systems from indefinite Helmholtz equation by standard FEM approximations. This is caused by negative eigenvalues of the associated linear system and divergence occurs under such two smoothers. To make up the problem, an improvement is introduced by the shifted Laplacian preconditioner [8, 9, 10], which preconditions the Helmholtz equation with a complex shifted operator

$\displaystyle-\Delta-(1+{\bf i}\beta)k^{2},$

(4.7)

where $\beta$ are free parameter. The use of a shift is important to guarantee multilevel convergence, whereas the multilevel method for the shifted operator only converges for a sufficiently large shift. But this contradicts the fact that the outer Krylov acceleration prefers a small shift. Based on the idea of adding a shift to the original problem, we consider the stable CIP-FEM which is consistent with original equation and adds a complex shift in the bilinear form.

Recalling that every two dimensional subspace of $2h$-harmonics $E^{\theta^{0}}_{2h}$ with $\theta^{0}\in\Theta_{\rm low}$ is left invariant under the $\omega$-JAC and GS-LEX relaxations, then the Fourier representation of smoother ${\boldsymbol{S}}_{h}={\boldsymbol{S}}^{J}_{h}$ or ${\boldsymbol{S}}^{GS}_{h}$ with respect to $E^{\theta^{0}}_{2h}$ can be written as

$$\left[\begin{array}[]{cc}\widetilde{{\boldsymbol{S}}}_{h}(\theta^{0})&0\\ 0&\widetilde{{\boldsymbol{S}}}_{h}(\theta^{1})\\ \end{array}\right],$$

(4.8)

where $\widetilde{{\boldsymbol{S}}}_{h}(\theta)$ is the smoother symbol derived in (4.5) and (4.6). The spectral radius of the smoother operator can be easily calculated since the above matrix is diagonal.

For simplicity, we concern on the weighted Jacobi relaxation. The frequency $\theta$ maximizing $|\widetilde{{\boldsymbol{S}}}^{J}_{h}(\theta)|$ over $(-\pi,\pi]$ can be calculated by its first and second derivatives, which reveal $\theta=0$ or $\theta=\pi$ maximizing $|\widetilde{{\boldsymbol{S}}}^{J}_{h}(\theta)|$. Hence, the spectral radius of ${\boldsymbol{S}}^{J}_{h}$ can be deduced as follows: $\rho({\boldsymbol{S}}^{J}_{h})=\max\{|\widetilde{{\boldsymbol{S}}}^{J}_{h}(0)|,|\widetilde{{\boldsymbol{S}}}^{J}_{h}(\pi)\}|,$ where

$\displaystyle|\widetilde{{\boldsymbol{S}}}^{J}_{h}(0)|=|1-\omega-\frac{\omega}{S}({\bf i}\gamma+R)|,\quad|\widetilde{{\boldsymbol{S}}}^{J}_{h}(\pi)|=|1-\omega-\frac{\omega}{S}({\bf i}\gamma-R)|.$

Since the parameter $\gamma$ in CIP-FEM may influence the pollution error, it is critical to make a suitable choice. For the case $t=\kappa h\leq 1$, it has been derived in [20] that for one dimensional problem there exists an optimal choice ${\bf i}\gamma_{o}=\frac{6\cos{t}-6+t^{2}\cos{t}+2t^{2}}{12(1-\cos{t})^{2}}$ such that the pollution error vanishes when the penalty parameter is chosen as $|{\bf i}\gamma-{\bf i}\gamma_{o}|\leq\frac{C}{\kappa^{2}h}$, where the constant $C$ is independent of $\kappa$ and $h$. The left graph of Figure 1 shows the Fourier symbols $\widetilde{{\boldsymbol{S}}}^{J}_{h}(\theta)$ for $\omega$-JAC smoother with $\gamma=\gamma_{o}$ and $\omega=0.6$. We find that $\widetilde{{\boldsymbol{S}}}^{J}_{h}(\theta)\geq 1$ always occur at the low frequencies, and small $t$ leads to a better relaxation. Similar phenomenon is observed for GS-LEX smoother in the right graph of Figure 1. Thus, for fixed wave number $\kappa$, both $\omega$-JAC and GS-LEX relaxations can be used as smoother on fine grid. Actually, the relaxation properties of these two smoothers are similar when ${\bf i}\gamma$ is a complex number.

When $\omega$-JAC smoother is utilized in the standard P1 FEM approximation, then

$\displaystyle\rho({\boldsymbol{S}}^{J}_{h})|_{\gamma=0}=\max\big{\{}\mid 1-\frac{\omega(12-t^{2})}{6-2t^{2}}\mid,\mid 1+\frac{3\omega t^{2}}{6-2t^{2}}\mid\big{\}}.$

(4.9)

Thus, for $t=\kappa h$ near $\sqrt{3}$, the error is amplified by the smoothing. For CIP-FEM, there is a (complex) shift in the bilinear form, therefore the amplification factor is reduced. This may permit the use of $\omega$-JAC smoother again. In fact, in our modified multilevel method, GMRES iteration is used on the intermediate grids (cf. [5]). In contrast with $\omega$-JAC and GS-LEX relaxations, the Fourier symbol can not be derived for GMRES smoothing. In the following, we will give some explanations for the performance of GMRES smoothing from the numerical approach.

For simplification, we consider the one dimensional Helmholtz equation on an interval $(0,10)$ with homogeneous Dirichlet boundary conditions. The piecewise P1 CIP-FEM discretization leads to a linear system with $N\times N$ coefficient matrix (cf. [20])

$\displaystyle{\boldsymbol{A}}^{C}_{h}=\frac{1}{h}\left(\begin{array}[]{ccccccc}2S-{\bf i}\gamma&R&{\bf i}\gamma&&\\ R&2S&R&{\bf i}\gamma&\\ {\bf i}\gamma&R&2S&R&{\bf i}\gamma\\ &\ddots&\ddots&\ddots&\ddots&\ddots&\\ &&{\bf i}\gamma&R&2S&R&{\bf i}\gamma\\ &&&{\bf i}\gamma&R&2S&R\\ &&&&{\bf i}\gamma&R&2S-{\bf i}\gamma\\ \end{array}\right)_{N\times N},$

(4.17)

where $N=10/h-1$. For $\kappa=200$, we apply the grid with mesh size $h=0.004$, i.e., $t=0.8$, and choose $\gamma=\gamma_{o}$. Let the vector $u_{0}=e^{{\bf i}\theta{\boldsymbol{x}}/h}$ be an initial choice for smoothing, where ${\boldsymbol{x}}=[x_{1},\cdots,x_{N-1}],x_{k}=x_{k-1}+h,x_{0}=0,k=1,\cdots,N-1$, $\theta\in(-\pi,\pi]$. We assume that ${\boldsymbol{S}}_{h}$ is the relaxation iteration matrix and $u_{1}={\boldsymbol{S}}_{h}u_{0}$ is the new vector after one step of smoothing. Then for fixed $\theta$, we obtain the amplification factor for one step of smoothing $\rho_{s}(\theta)=\|u_{1}\|/\|u_{0}\|$, where $\|\cdot\|$ stands for the Euclidean norm. Figure 2 shows the amplification factors of three different kinds of smoothing strategies: GS, Jacobi and GMRES relaxations. We can see that for high and low frequencies $\theta$, the amplification factor of GMRES relaxation is always smaller than that of GS and Jacobi relaxations. The GMRES relaxation can still be convergent even when the other two smoothers lead to a divergent relaxation. The similar phenomenon can also be observed for other choices of $t$ and $\gamma$. Therefore, when standard Jacobi or GS relaxation fails as a smoother, we can replace this with the GMRES smoothing.

From the above analysis, we can see that for the relaxation of linear system for Helmholtz problem, the standard Jacobi or GS relaxation can be performed on the fine grids, but they fail on the coarse grids. In particular, the GMRES relaxation is efficient for smoothing on the intermediate coarse grids. In the next section, the correction scheme will be taken into account to obtain a more realistic convergence analysis of the multilevel method.

For more details about convergence property of multilevel method in Algorithm 3.1, we will focus on the LFA of two-level method. The LFA of three-level method will also be mentioned. For simplicity, we consider the case without postsmoothing. Since the Fourier symbol can not be obtained for GMRES smoothing, we only consider the two and three level methods with $\omega$-JAC or GS-LEX relaxation. Then the iteration operator of Algorithm 3.1 in this case is the nonsymmetric version and can be derived as in (2.7). Three cases of iteration operator will be analyzed in the following: multilevel methods by shifted Laplace approach, multilevel methods with stable CIP-FEM coarse grid corrections, and multilevel methods with stable CIP-FEM for fine grid smoothing and coarse grid corrections.

The iteration operator of two-level method can be written as $(I-T_{1})(I-T_{0})$, and the operator $T_{l}$ is chosen for the corresponding algorithm. Since standard FEM is applied on the finest grid, we have $T_{l}=\mu_{l}R_{l}A^{F}_{l}P^{F}_{l}$, where $R_{l}$ is smoother determined by the algorithm. When CIP-FEM is used for smoothing on the grids ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{0}$ respectively, the iteration matrix is given by

$${\boldsymbol{M}}^{C}_{2}=\Big{(}{\boldsymbol{I}}_{1}-\mu_{1}({\boldsymbol{I}}_{1}-{\boldsymbol{S}}^{C}_{1})({\boldsymbol{A}}^{C}_{1})^{-1}{\boldsymbol{A}}^{F}_{1}\Big{)}\Big{(}{\boldsymbol{I}}_{1}-\mu_{0}{\boldsymbol{I}}^{1}_{0}({\boldsymbol{A}}^{C}_{0})^{-1}{\boldsymbol{I}}^{0}_{1}{\boldsymbol{A}}^{F}_{1}\Big{)}.$$

Here, for $l\geq 0$, ${\boldsymbol{S}}^{C}_{l}={\boldsymbol{I}}_{l}-{\boldsymbol{R}}^{C}_{l}{\boldsymbol{A}}^{C}_{l}$ is smoothing relaxation matrix, ${\boldsymbol{I}}_{l}$ with the same size as ${\boldsymbol{A}}^{F}_{l}$ is identity matrix, ${\boldsymbol{I}}^{s}_{l}(s>l)$ is prolongation matrix from level $l$ to $s$, ${\boldsymbol{I}}^{s}_{l}(s<l)$ is restriction matrix from level $l$ to $s$, and ${\boldsymbol{R}}^{C}_{l}$ stands for matrix representation of smoother $R_{l}^{C}$ with CIP-FEM approximation. Let ${\boldsymbol{M}}^{SL}_{2}$ (cf. [4]) stand for two-level iteration matrix, taking the smoothing based on shifted Laplace operator (4.7) with standard FEM approximation. Similarly, when applying standard FEM for smoothing on ${\mathcal{T}}_{1}$ and CIP-FEM for correction on ${\mathcal{T}}_{0}$, the associated iteration matrix is derived as

$${\boldsymbol{M}}^{FC}_{2}=\Big{(}{\boldsymbol{I}}_{1}-\mu_{1}({\boldsymbol{I}}_{1}-{\boldsymbol{S}}^{F}_{1})\Big{)}\Big{(}{\boldsymbol{I}}_{1}-\mu_{0}{\boldsymbol{I}}^{1}_{0}({\boldsymbol{A}}^{C}_{0})^{-1}{\boldsymbol{I}}^{0}_{1}{\boldsymbol{A}}^{F}_{1}\Big{)},$$

where ${\boldsymbol{S}}^{F}_{1}={\boldsymbol{I}}_{1}-{\boldsymbol{R}}^{F}_{1}{\boldsymbol{A}}^{F}_{1}$, and ${\boldsymbol{R}}^{F}_{1}$ stands for smoothing matrix representation of $R_{1}$ with standard FEM approximation.

Since every two dimensional subspace (4.3) of $2h$-harmonics $E^{\theta^{0}}_{2h_{1}}$ with $\theta^{0}\in(-\pi/2,\pi/2]$ is left invariant under $\omega$-JAC or GS-LEX smoothing operator and correction operator, the representation of two-level iteration matrix of ${\boldsymbol{M}}^{C}_{2}$ on $E^{\theta^{0}}_{2h_{1}}$ is given by a $2\times 2$ matrix as follows:

	$\displaystyle\widetilde{{\boldsymbol{M}}}^{C}_{2}$	$\displaystyle=$	$\displaystyle\left[\widetilde{{\boldsymbol{I}}}_{1}-\mu_{1}\left(\widetilde{{\boldsymbol{I}}}_{1}-\left[\begin{array}[]{c}\widetilde{{\boldsymbol{S}}}_{1}^{C}(\theta^{0})\\ \widetilde{{\boldsymbol{S}}}_{1}^{C}(\theta^{1})\\ \end{array}\right]_{D}\right)\left[\begin{array}[]{cc}\widetilde{{\boldsymbol{A}}}_{1}^{C}(\theta^{0})\\ \widetilde{{\boldsymbol{A}}}_{1}^{C}(\theta^{1})\\ \end{array}\right]_{D}^{-1}\left[\begin{array}[]{cc}\widetilde{{\boldsymbol{A}}}_{1}^{F}(\theta^{0})\\ \widetilde{{\boldsymbol{A}}}_{1}^{F}(\theta^{1})\\ \end{array}\right]_{D}\right]$		(4.31)
			$\displaystyle\cdot\left[\widetilde{{\boldsymbol{I}}}_{1}-\mu_{0}\left[\begin{array}[]{c}\widetilde{{\boldsymbol{I}}}_{0}^{1}(\theta^{0})\\ \widetilde{{\boldsymbol{I}}}_{0}^{1}(\theta^{1})\\ \end{array}\right]\widetilde{{\boldsymbol{A}}}^{C}_{0}(2\theta^{0})^{-1}\left[\begin{array}[]{c}\widetilde{{\boldsymbol{I}}}_{1}^{0}(\theta^{0})\\ \widetilde{{\boldsymbol{I}}}_{1}^{0}(\theta^{1})\\ \end{array}\right]^{t}\left[\begin{array}[]{c}\widetilde{{\boldsymbol{A}}}_{1}^{F}(\theta^{0})\\ \widetilde{{\boldsymbol{A}}}_{1}^{F}(\theta^{1})\\ \end{array}\right]_{D}\right],$

where $\widetilde{{\boldsymbol{I}}}_{1}$ is $2\times 2$ identity matrix and the subscript-_D denotes the transformation of a vector into a diagonal matrix. Similarly, the representation $\widetilde{{\boldsymbol{M}}}^{SL}_{2}$ of two-level iteration matrix based on shifted Laplace operator, can be easily obtained (cf. [4]). For the iteration operator ${\boldsymbol{M}}^{FC}_{2}$, its representation on $E^{\theta^{0}}_{2h_{1}}$ is given by

	$\displaystyle\widetilde{{\boldsymbol{M}}}^{FC}_{2}$	$\displaystyle=$	$\displaystyle\left[\widetilde{{\boldsymbol{I}}}_{1}-\mu_{1}\left(\widetilde{{\boldsymbol{I}}}_{1}-\left[\begin{array}[]{c}\widetilde{{\boldsymbol{S}}}_{1}^{F}(\theta^{0})\\ \widetilde{{\boldsymbol{S}}}_{1}^{F}(\theta^{1})\\ \end{array}\right]_{D}\right)\right]$		(4.41)
			$\displaystyle\cdot\left[\widetilde{{\boldsymbol{I}}}_{1}-\mu_{0}\left[\begin{array}[]{c}\widetilde{{\boldsymbol{I}}}_{0}^{1}(\theta^{0})\\ \widetilde{{\boldsymbol{I}}}_{0}^{1}(\theta^{1})\\ \end{array}\right]\widetilde{{\boldsymbol{A}}}^{C}_{0}(2\theta^{0})^{-1}\left[\begin{array}[]{c}\widetilde{{\boldsymbol{I}}}_{1}^{0}(\theta^{0})\\ \widetilde{{\boldsymbol{I}}}_{1}^{0}(\theta^{1})\\ \end{array}\right]^{t}\left[\begin{array}[]{c}\widetilde{{\boldsymbol{A}}}_{1}^{F}(\theta^{0})\\ \widetilde{{\boldsymbol{A}}}_{1}^{F}(\theta^{1})\\ \end{array}\right]_{D}\right].$

Then the spectral radius of $\widetilde{{\boldsymbol{M}}}^{SL}_{2},\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ and $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ for different $\theta^{0}\in\Theta_{\rm low}$ can be obtained analytically and numerically.

The left graph of Figure 3 shows the spectral radius of $\widetilde{{\boldsymbol{M}}}^{SL}_{2}$ with $\beta=0.5,0.3$ and $0.1$ over $\theta^{0}\in\Theta_{\rm low}$, and $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$, $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ under $\omega$-JAC relaxation for $t=0.8$ on the finest grid. In order to make some comparisons between different algorithms, the penalty parameter ${\bf i}\gamma$ in the CIP-FEM is always chosen as a complex number in the following if there is no special notation. If no confusion is possible, we will always denote $t={\it constant}$ for $t=\kappa h_{L}$ on the finest grid. The parameter in $\omega$-JAC relaxation is always chosen as $\omega=0.6$ in this section. We observe that for most of the frequencies $\theta^{0}\in\Theta_{\rm low}$ the spectral radius is smaller than one. The spectral radius tends to larger than one only under a few frequencies. Actually, the appearance of such a resonance is caused by the coarse grid correction and originates from the inversion of the coarse grid discretization symbol, $\widetilde{{\boldsymbol{A}}}^{C}_{0}(2\theta^{0})$ in (4.31) and (4.41) for instance. The minimum of $|\widetilde{{\boldsymbol{A}}}^{C}_{0}(2\theta^{0})|$ maximizes the spectral radius of $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ and $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ for fixed $t$ and ${\bf i}\gamma$.

Moreover, for different choices of parameter $\beta$ in shifted Laplace operator, the properties of $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ and $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ with ${\bf i}\gamma={\bf i}\gamma_{o}+0.01{\bf i}$ on the finest grid always perform better than that of $\widetilde{{\boldsymbol{M}}}^{SL}_{2}$. Although the spectral radius of $\widetilde{{\boldsymbol{M}}}^{SL}_{2}$ with $\beta=0.1$ is almost equivalent to that of $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ over the frequencies spreading near zero, there is a relatively large resonance causing by the coarse grid correction. Small $\beta$ in (4.7) deteriorates the coarse grid correction. However, large $\beta$ amplifies the corresponding spectral radius over the frequencies spreading near zero. We can refer to a recent work [4] about the choice of minimal complex shift parameter in shifted Laplace preconditioner. The choice of ${\bf i}\gamma$ in $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ or $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ would also influence its spectral radius. The right graph of Figure 3 shows the spectral radius of $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ with different ${\bf i}\gamma$ under $\omega$-JAC relaxation for $t=0.8$. The influence of ${\bf i}\gamma$ in $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ is similar to $\beta$ in shifted Laplace preconditioner. For other small $t<1$, the properties of $\widetilde{{\boldsymbol{M}}}^{SL}_{2}$, $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ and $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ can also be observed as above. Moreover, when $t=\kappa h$ on the finest grid is small enough, the spectral radiuses of $\widetilde{{\boldsymbol{M}}}^{SL}_{2}$, $\widetilde{{\boldsymbol{M}}}^{FC}_{2}$ and $\widetilde{{\boldsymbol{M}}}^{C}_{2}$ can all be smaller than one over all $\theta^{0}\in\Theta_{\rm low}$. Due to the fact that our aim is to study the influences of smoothing corrections in two- and three-level methods, in the following we do not always assume $t$ to be sufficiently small.