A Modified SSOR Preconditioning Strategy for Helmholtz Equations

Abstract

The finite difference method discretization of Helmholtz equations usually leads to the large spare linear systems. Since the coefficient matrix is frequently indefinite, it is difficult to solve iteratively. In this paper, a modified symmetric successive overrelaxation (MSSOR) preconditioning strategy is constructed based on the coefficient matrix and employed to speed up the convergence rate of iterative methods. The idea is to increase the values of diagonal elements of the coefficient matrix to obtain better preconditioners for the original linear systems. Compared with SSOR preconditioner, MSSOR preconditioner has no additional computational cost to improve the convergence rate of iterative methods. Numerical results demonstrate that this method can reduce both the number of iterations and the computational time significantly with low cost for construction and implementation of preconditioners.

1. Introduction

The finite difference method is one of the most effective and popular techniques in computational electromagnetics and seismology, such as time-harmonic wave propagations, scattering phenomena arising in acoustic and optical problems, and electromagnetics scattering from a large cavity. More information about applications of this method in electromagnetics can be found in [1–5].

In this paper, we focus on the following form of the complex Helmholtz equation: −Δ𝑢−𝑝𝑢+𝑖𝑞𝑢=𝑓inΩ,𝑢=𝑔on𝜕Ω. (1.1) Here Ω is a bounded region in ℝ². 𝑝 and 𝑞 are real continuous coefficient functions on Ω, while 𝑓 and 𝑔 are given continuous functions on Ω and 𝜕Ω, respectively.

To conveniently find numerical solutions of (1.1), the Laplace operator is approximated by using the second-order accurate 5-point difference stencil:𝐿_ℎ1=̇ℎ²111−41_ℎ. (1.2) Making use of the above stencil, the following linear system is obtained: 𝒜𝑥=𝐴−ℎ²𝐷+𝑖ℎ²𝐸𝑥=𝑏, (1.3) where 𝐷(𝐸) is a diagonal matrix whose diagonal elements are just the values of 𝑝(𝑞) at the mesh points and 𝐴 is the symmetric positive definite 𝑀-matrix arising from the discrete Laplace operator and is of the block tridiagonal form𝐺𝐴=¹𝐹₂𝐸₂𝐺₂⋱⋱⋱𝐹_𝑚𝐸_𝑚𝐺_𝑚, (1.4)

with𝐺_𝑘=4−1−14⋱⋱⋱−1−14_{𝑛𝑘×𝑛𝑘}𝐸,𝑘=1,2,…,𝑚,_𝑘=𝐹^𝑇_𝑘=𝑂_{𝑛𝑘×𝑝𝑘}−𝐼_𝑛𝑘𝑂_{𝑛𝑘×𝑞𝑘},𝑝_𝑘,𝑞_𝑘𝑝≥0,_𝑘+𝑞_𝑘=𝑛_𝑘−1−𝑛_𝑘if𝑛_𝑘≤𝑛_𝑘−1,𝑂_{𝑛𝑘−1×𝑠𝑘}−𝐼_𝑛𝑘−1𝑂_{𝑛𝑘−1×𝑡𝑘}^𝑇,𝑠_𝑘,𝑡_𝑘𝑠≥0,_𝑘+𝑡_𝑘=𝑛_𝑘−𝑛_𝑘−1if𝑛_𝑘>𝑛_𝑘−1,𝑘=2,3,…,𝑚. (1.5)

Obviously, from the linear systems (1.3), it is not difficult to find that the matrix 𝒜 is a complex symmetric coefficient matrix. Matrix 𝒜 becomes highly indefinite and ill-conditioned as 𝑝 is a sufficiently large positive number. So, large amount of computation times and memory are needed in order to solve the linear systems (1.3) efficiently.

As is well known, direct methods and iterative methods can be employed to solve the linear systems (1.3). The former is widely employed when the order of the coefficient matrix 𝒜 is not too large and is usually regarded as robust methods. The memory and the computational requirements for solving the large sparse linear systems may seriously challenge the most efficient direct solution method available today. Currently, the latter employed to solve the large sparse linear systems is popular. The reason is that iterative methods are easier to implement efficiently on high performance computers than direct methods. In practice, a natural choice is that we make use of iterative methods instead of direct methods to solve the large sparse linear systems.

At present, Krylov subspace methods are considered as one kind of the important and efficient iterative techniques for solving the large sparse linear systems because these methods are cheap to be implemented and are able to fully exploit the sparsity of the coefficient matrix. It is well known that the convergence speed of Krylov subspace methods depends on the distribution of the eigenvalues of the coefficient matrix [6]. When the coefficient matrix is typically extremely ill-conditioned and highly indefinite, the convergence of Krylov subspace methods can be unacceptably slow. In this case, Krylov subspace methods are not competitive without a good preconditioner. That is, preconditioning technique is a key ingredient for the success of Krylov subspace methods in applications. The idea of preconditioning technique is based on consideration of the linear system with the same solution as the original equation. The problem is that each preconditioning technique is suited for a different type of problem. Until current days, no robust preconditioning technique appears for all or at least much types of problems. Finding a good preconditioner to solve a given large sparse linear systems is often viewed as a combination of art and science.

In recent years, a great deal of effort has been invested in solving indefinite linear systems from the discrete Helmholtz equations. Most of the work has been aimed at developing effective preconditioning techniques. In general, there exist two classes of preconditioners for Helmholtz equations: the “operator-based" preconditioning technique and the “matrix-based" preconditioning technique.

The former is built based on an operator, such as the Laplace preconditioner [2, 3, 7–9], Analytic ILU [10], the Separation-of-variables [11]. The purpose of this class of preconditioners is that the spectrum of the corresponding preconditioned matrix is favorably clustered. Its advantage is that this operator does not have to be a representation of the inverse of the Helmholtz operator.

The latter is established based on an approximation of the inverse of the coefficient matrix. For this class, one of the natural and simplest ways of structuring a preconditioner is to employ a diagonal or block diagonal of the coefficient matrix as a preconditioner [12]. The above two diagonal preconditioners have no remarkable reduce with respect to the iterative number and CPU time. Another one of the simplest preconditioners is to perform an incomplete factorization (ILU) of the coefficient matrix [1]. The main idea of ILU factorizations depends on the implementation of Gaussian elimination which is used, see the survey [13] and the related references therein.

When the coefficient matrix of the linear systems (1.1) is complex symmetric and indefinite, it is difficult to solve iteratively. Using the symmetric successive overrelaxation (SSOR) as a preconditioner preserves the symmetry of the iterative matrix and also taking little initialization cost, which in some cases makes it preferable over other factorization methods such as ILU. So far, some variant SSOR preconditioning techniques have been proposed to improve the convergence rate of the corresponding iterative method for solving the linear systems. Mazzia and Alan [14] introduced a shift parameter to develop a shifted SSOR preconditioner for solving the linear systems from an electromagnetics application. Bai in [15] used a (block) diagonal matrix instead of the diagonal matrix of the coefficient matrix to establish a modified (block) SSOR preconditioner for the second-order selfadjoint elliptic boundary value problems. Chen et al. [16] used a diagonal matrix with a relaxation parameter instead of the diagonal matrix of the coefficient matrix to establish a modified block SSOR preconditioner for the Biot’s consolidation equations. Ran and Yuan in [17] also discussed a class of modified (block) SSOR preconditioners for linear systems from steady incompressible viscous flow problems. We refer the reader to [18, 19] for a general discussion.

Although the SSOR preconditioner with Krylov subspace methods can improve convergence improvement, the disadvantage of the SSOR preconditioner is that the convergence rate may still remain unsatisfactorily slow in some cases, especially in indefinite linear systems. In this paper, a modified symmetric successive overrelaxation (MSSOR) preconditioning strategy is presented, which can significantly improve the convergence speed and CPU time. Our motivation for this method arises from the solution of complex symmetric and indefinite linear systems from the discrete Helmholtz equations. The idea is to increase the values of diagonal elements of the coefficient matrix to obtain better preconditioners for the original linear systems, which is different from [14–17]. This modification does not require any significant computational cost as compared with the original SSOR preconditioner and also requires no additional storage cost.

The remainder of this paper is organized as follows. In Section 2, the MSSOR preconditioner for solving the resulting linear system is presented. In Section 3, numerical experiments are given to illustrate the efficiency of the presented preconditioner. Finally, in Section 4 some conclusions are drawn.

2. Modified SSOR Preconditioner

To improve the convergence rate of iterative methods, an appropriate preconditioner should be incorporated. That is, it is often preferable to solve the preconditioned linear system as follows:𝑃⁻¹𝒜𝑥=𝑃⁻¹𝑏, (2.1) where 𝑃, called the preconditioner, is a nonsingular matrix. The choice of the preconditioner 𝑃 plays an important role in actual implements. In general, the preconditioner 𝑃 is chosen such that the condition number of the preconditioned matrix 𝑃⁻¹𝒜 is less than that of the original matrix 𝒜. Based on the excellent survey of [13] by Benzi, a good preconditioner should meet the following requirements:(1)the preconditioned system should be easy to solve;(2)the preconditioner should be cheap to construct and apply.

Certainly, the best choice for 𝑃⁻¹ is the inverse of 𝒜. However, it is unpractical in actual implements because the cost of the computation of 𝒜⁻¹ may be high. If 𝒜 is a symmetric positive definite matrix, the approximation of 𝒜⁻¹ can be replaced by SSOR or multigrid. However, in fact, the Helmholtz equation often results in an indefinite linear system, for which SSOR or multi-grid may be not guaranteed to converge.

To introduce the modified SSOR preconditioner, a brief review of the classical and well-known SSOR preconditioner is needed. The SSOR preconditioner is established by the SSOR iterative method, which is a symmetric version of the well-known SOR iterative method. Based on matrix splitting, the coefficient matrix 𝒜 is split as follows: 𝒜=𝒟+ℒ+ℒ^𝑇, (2.2) where 𝒟 and ℒ are the diagonal parts and strictly lower triangular of 𝒜. According to the foregoing matrix splitting (2.2), the standard SSOR preconditioner [6, 20] is defined by 𝑃_SSOR=(𝒟+ℒ)𝒟⁻¹𝒟+ℒ^𝑇. (2.3) It is difficult to show theoretically the behavior of a preconditioner when the coefficient matrix 𝒜 is a large, sparse, and symmetric indefinite. The SSOR iterative method is not convergent, but 𝑃_SSOR may be still used as a preconditioner. By a simple modification on the original indefinite linear systems (1.3), we establish the following coefficient matrix: 𝒜=𝐴+ℎ²𝐷+ℎ²𝐸𝑖=𝒟+ℒ+ℒ^𝑇, (2.4) where 𝒟=diag(𝐴)+ℎ²𝐷+ℎ²𝐸𝑖. (2.5) Obviously, 𝒜 is a symmetric and positive stable 𝐻-matrix. To increase the values of diagonal elements of the coefficient matrix to obtain better preconditioners for the original linear systems and reduce computation times and amount of memory, based on (2.4), the MSSOR preconditioner is defined by 𝑃_MSSOR=𝒟+ℒ𝒟⁻¹𝒟+ℒ^𝑇, (2.6) with𝒟=diag𝒟,𝒟=𝑑_𝑖𝑖(𝑖=1,2,…,𝑛). (2.7) This idea is based on an absolute diagonal scaling technique, which is cheap and easy to implement.

Since the coefficient matrix of the linear systems (1.3) is neither positive definite nor Hermitian with 𝑝 being a sufficiently large positive number, the Conjugate Gradient (CG) method [21] may breakdown. To solve the complex symmetric linear systems, van der Vorst and Melissen [22] proposed the conjugate orthogonal conjugate gradient (COCG) method, which is regarded as an extension of CG method.

To solve the linear systems (1.3) efficiently, (1.3) is transformed into the following form with the preconditioner 𝑃_MSSOR, that is, 𝑃_MSSOR𝐴−ℎ²𝐷+𝑖ℎ²𝐸𝑥=𝑃_MSSOR𝑏. (2.8) Then the MSSOR preconditioned COCG (PCOCG) method can be employed to solve the preconditioned linear systems (2.8).

In the following, we give the MSSOR preconditioned COCG method for solving the linear systems (2.8). The MSSOR preconditioned COCG (PCOCG) algorithm is described as follows

Algorithm PCOCG [22]: given an initial guess 𝑥₀(1)𝑣₀=𝑏−𝒜𝑥₀;(2)𝑝₋₁=0;𝛽₋₁=0;(3)𝑤₀=𝑃⁻¹_MSSOR𝑣₀;𝜌₀=𝑣^𝑇₀𝑤₀;(4)for 𝑖=0,1,2,…(5)𝑝_𝑖=𝑤_𝑖+𝛽_𝑖−1𝑝_𝑖−1;(6)𝑢_𝑖=𝒜𝑝_𝑖;(7)𝜇_𝑖=𝑢^𝑇_𝑗𝑝_𝑗; if 𝜇_𝑖=0 then quit (failure);(8)𝛼=𝜌_𝑖/𝜇_𝑖;(9)𝑥_𝑖+1=𝑥_𝑖+𝛼𝑝_𝑖;𝑣_𝑖+1=𝑣_𝑖−𝛼𝑢_𝑖(10)if 𝑥_𝑖+1 is accurate enough, then quit (convergence);(11)𝑤_𝑖+1=𝑃⁻¹_MSSOR𝑣_𝑖+1;(12)𝜌_𝑖+1=𝑣^𝑇_𝑖+1𝑤_𝑖+1; if |𝜌_𝑖+1| too small, then quit (failure);(13)𝛽_𝑖=𝜌_𝑖+1/𝜌_𝑖;(14)End for 𝑖.

It is not difficult to find that the main computation of algorithm PCOCG involves one matrix-vector multiplication and two triangular linear systems. These computations are very easy to implement. The main advantage is no extra computational cost in construction of MSSOR preconditioner.

Note that the transpose in all dot products in this algorithm is essential [23]. Meanwhile, note that two different breakdowns of this algorithm may occur: one is that if 𝑣^𝑇_𝑖+1𝑃⁻¹_MSSOR𝑣_𝑖+1 is too small, but 𝑣_𝑖+1 exists (line 12), algorithm PCOCG breaks down and the other is that when the search direction 𝑝_𝑖≠0, but 𝑝^𝑇_𝑖𝒜𝑝_𝑖=0 (line 7), algorithm PCOCG breaks down. The breakdown can be fixed to some extent by restarting the process [22], such as the restarted process in GMRES [24]. However, breakdown scarcely happens in the actual computation of the Helmholtz equation.

3. Numerical Experiments

In this section, some numerical experiments are given to demonstrate the performance of both preconditioner 𝑃_SSOR and preconditioner 𝑃_MSSOR for solving the Helmholtz equation.

Example 3.1 (see [25]). Consider the following complex Helmholtz equation: −Δ𝑢−𝑝𝑢+𝑖𝑞𝑢=𝑓inΩ,𝑢=𝑔on𝜕Ω, (3.1) where Ω=[0,1]×[0,1] and 𝑝≥0 and 𝑞≥0 are real constants. Discretizing (3.1) with the approach above in introduction, we obtain the complex symmetric indefinite linear systems 𝒜𝑥=(𝐴−ℎ²𝐷+𝑖ℎ²𝐸)𝑥=𝑏, and 𝑓 and 𝑔 are adjusted such that 𝑏=𝒜𝑒(𝑒=(1,1,…,1)^𝑇).

All tests are started from the zero vector, preformed in MATLAB with machine precision 10⁻¹⁶. The COCG iteration terminates if the relative residual error satisfies ‖𝑟^(𝑘)‖₂/‖𝑟⁽⁰⁾‖₂<10⁻⁶ or the iteration number is more than 500.

In Tables 1, 2, 3, and 4, we present some iteration results to illustrate the convergence behaviors of the COCG method preconditioned by 𝑃_SSOR and 𝑃_MSSOR to solve the complex symmetric indefinite linear systems with the different values of 𝑝 and 𝑞. In Tables 1–4, (⋅,⋅) denotes the values of 𝑝 and 𝑞. “CPU(s)’’ denotes the time (in seconds) required to solve a problem. “IT’’ denotes the number of iteration.

Table 1 

Iteration numbers for the two-dimensional Helmholtz equations for ℎ=1/19 meshes, using COCG with the preconditioner 𝑃_SSOR and 𝑃_MSSOR for solving the complex symmetric indefinite linear systems.

Table 2 

Iteration numbers for the two-dimensional Helmholtz equations for ℎ=1/34 meshes, using COCG with the preconditioner 𝑃_SSOR and 𝑃_MSSOR for solving the complex symmetric indefinite linear systems.

Table 3 

Iteration numbers for the two-dimensional Helmholtz equations for ℎ=1/64 and 𝑝=4100, using COCG with the preconditioner 𝑃_SSOR and 𝑃_MSSOR.

Table 4 

Iteration numbers for the two-dimensional Helmholtz equations for ℎ=1/119 and 𝑞=2000, using COCG with the preconditioner 𝑃_SSOR and 𝑃_MSSOR.

From Tables 1–4, it is not difficult to find that when the COCG method preconditioned by 𝑃_SSOR and 𝑃_MSSOR is used to solve the complex symmetric indefinite linear systems, the convergence rate of the preconditioner 𝑃_MSSOR is more efficient than that of the preconditioner 𝑃_SSOR by the iteration numbers and CPU time. That is, the preconditioner 𝑃_MSSOR outperforms the preconditioner 𝑃_SSOR under certain conditions. Compared with the preconditioner 𝑃_SSOR, the preconditioner 𝑃_MSSOR may be the “preferential” choice under certain conditions.

4. Conclusions

In this paper, MSSOR preconditioned COCG algorithm has been applied for solving the complex symmetric indefinite systems arising from Helmholtz equations. Due to the reduction of the iteration numbers and CPU time, the MSSOR preconditioner presented is feasible and effective. Without extra costs, MSSOR preconditioner is more efficient than SSOR preconditioner.

Acknowledgments

The authors are grateful to the referees and the editors for their helpful suggestions to improve the quality of this paper. The research of this author was supported by NSFC Tianyuan Mathematics Youth Fund (11026040).