Journal of Applied MathematicsVolume 2012 (2012), Article ID 365124, 9 pagesdoi:10.1155/2012/365124

Research Article

A Modified SSOR Preconditioning Strategy for Helmholtz Equations

Shi-Liang Wu¹ and Cui-Xia Li²

¹College of Mathematics, Chengdu University of Information Technology, Chengdu 610225, China²School of Mathematics and Statistics, Anyang Normal University, Anyang 455002, China

Received 22 August 2011; Revised 7 November 2011; Accepted 18 November 2011

Academic Editor: Kok Kwang Phoon

Copyright © 2012 Shi-Liang Wu and Cui-Xia Li. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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: 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: Making use of the above stencil, the following linear system is obtained: 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

with

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: 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: 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 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 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: where 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 with 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 , that is, 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);(3);(4)for (5);(6);(7); if then quit (failure);(8);(9)(10)if is accurate enough, then quit (convergence);(11);(12); if too small, then quit (failure);(13);(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 is too small, but exists (line 12), algorithm PCOCG breaks down and the other is that when the search direction , but (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 and preconditioner for solving the Helmholtz equation.

Example 3.1 (see [25]). Consider the following complex Helmholtz equation: where and and 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 .

All tests are started from the zero vector, preformed in MATLAB with machine precision . The COCG iteration terminates if the relative residual error satisfies 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 and 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 meshes, using COCG with the preconditioner and for solving the complex symmetric indefinite linear systems.

(p, q)

h
Preconditioner
(800, 10)
(800, 20)
(800, 30)
(800, 40)

IT
246
242
239
214

CPU(s)
0.8438
0.8125
0.7813
0.7344

IT
138
133
126
117

CPU(s)
0.25
0.2188
0.2031
0.1875

Table 2: Iteration numbers for the two-dimensional Helmholtz equations for meshes, using COCG with the preconditioner and for solving the complex symmetric indefinite linear systems.

(p, q)

h
Preconditioner
(1400, 40)
(1500, 40)
(1600, 40)
(1700, 40)

IT
230
274
310
336

CPU(s)
3.2813
3.9375
4.8438
4.7188

IT
214
228
237
249

CPU(s)
1.4375
1.5156
1.6563
1.6250

Table 3: Iteration numbers for the two-dimensional Helmholtz equations for and , using COCG with the preconditioner and .

q

h
Preconditioner
100
120
150
160

IT
471
456
374
349

CPU(s)
32.9344
30.9688
27.9844
26.0156

IT
373
366
347
326

CPU(s)
12.4688
12.2500
11.6094
10.9063

Table 4: Iteration numbers for the two-dimensional Helmholtz equations for and , using COCG with the preconditioner and .

p

h
Preconditioner
15000
15500
16000
16500

IT
148
153
147
165

CPU(s)
55.3594
57.4531
57.8594
52.4844

IT
134
137
139
144

CPU(s)
22.625
22.9063
23.2813
23.7344

From Tables 1–4, it is not difficult to find that when the COCG method preconditioned by and is used to solve the complex symmetric indefinite linear systems, the convergence rate of the preconditioner is more efficient than that of the preconditioner by the iteration numbers and CPU time. That is, the preconditioner outperforms the preconditioner under certain conditions. Compared with the preconditioner , the preconditioner 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).

References

C.-H. Guo, “Incomplete block factorization preconditioning for linear systems arising in the numerical solution of the Helmholtz equation,” Applied Numerical Mathematics, vol. 19, no. 4, pp. 495–508, 1996.

Y. A. Erlangga, C. Vuik, and C. W. Oosterlee, “On a class of preconditioners for solving the Helmholtz equation,” Applied Numerical Mathematics, vol. 50, no. 3-4, pp. 409–425, 2004.

A. Bayliss, C. I. Goldstein, and E. Turkel, “An iterative method for the Helmholtz equation,” Journal of Computational Physics, vol. 49, no. 3, pp. 443–457, 1983.

G. Bao and W. W. Sun, “A fast algorithm for the electromagnetic scattering from a large cavity,” SIAM Journal on Scientific Computing, vol. 27, no. 2, pp. 553–574, 2005.

Y. Wang, K. Du, and W. W. Sun, “Preconditioning iterative algorithm for the electromagnetic scattering from a large cavity,” Numerical Linear Algebra with Applications, vol. 16, no. 5, pp. 345–363, 2009.

Y. Saad, Iterative Methods for Sparse Linear Systems, PWS, Boston, Mass, USA, 1996.

J. Gozani, A. Nachshon, and E. Turkel, “Conjugate gradient comipled with multi-grid for an indefinite problem,” in Advances in Computer Methods for Partial Differential Equations, R. Vichnevestsky and R. S. Tepelman, Eds., vol. 5, pp. 425–427, IMACS, New Brunswick, NJ, USA, 1984.

A. L. Laird, “Preconditioned iterative solution of the 2D Helmholtz equation,” First Year's Report 02/12, Hugh's College, Oxford, UK, 2002.

Y. A. Erlangga, C. Vuik, and C. W. Oosterlee, “Comparison of multigrid and incomplete LU shifted-Laplace preconditioners for the inhomogeneous Helmholtz equation,” Applied Numerical Mathematics, vol. 56, no. 5, pp. 648–666, 2006.

M. Gander, “AILU for Helmholtz problems: a new preconditioner based on the analytic parabolic factorization,” Journal of Computational Acoustics, vol. 9, no. 4, pp. 1499–1506, 2001.

R. E. Plessix and W. A. Mulder, “Separation-of-variables as a preconditioner for an iterative Helmholtz solver,” Applied Numerical Mathematics, vol. 44, no. 3, pp. 385–400, 2003.

M. M. M. Made, R. Beauwens, and G. Warzée, “Preconditioning of discrete Helmholtz operators perturbed by a diagonal complex matrix,” Communications in Numerical Methods in Engineering, vol. 16, no. 11, pp. 801–817, 2000.

M. Benzi, “Preconditioning techniques for large linear systems: a survey,” Journal of Computational Physics, vol. 182, no. 2, pp. 418–477, 2002.

F. Mazzia and R. Alan McCoy, “Numerical experimetns with a shifted SSOR preconditioner for symmetric matrices,” type TR/PA/98/12, CERFACS, 1998.

Z.-Z. Bai, “Mofidied block SSOR preconditioners for symmetric positive definite linear systems,” Annals of Operations Research, vol. 103, pp. 263–282, 2001.

X. Chen, K. C. Toh, and K. K. Phoon, “A modified SSOR preconditioner for sparse symmetric indefinite linear systems of equations,” International Journal for Numerical Methods in Engineering, vol. 65, no. 6, pp. 785–807, 2006.

Y.-H. Ran and L. Yuan, “On modified block SSOR iteration methods for linear systems from steady incompressible viscous flow problems,” Applied Mathematics and Computation, vol. 217, no. 7, pp. 3050–3068, 2010.

O. Axelsson, “A generalized SSOR method,” BIT, vol. 18, pp. 443–467, 1972.

A. Hadjidimos, “Successive overrelaxation (SOR) and related methods,” Journal of Computational and Applied Mathematics, vol. 123, no. 1-2, pp. 177–199, 2000.

O. Axelsson, Iterative Solution Methods, Cambridge University Press, New York, NY, USA, 1995.

M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” Journal of Research of the National Bureau of Standards, vol. 49, pp. 409–436, 1952.

H. A. van der Vorst and J. B. M. Melissen, “A Petrov-Galerkin type method for solving , where A is symmetric complex,” IEEE Transactions on Magnetics, vol. 26, pp. 706–708, 1990.

R. W. Freund, “Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices,” SIAM Journal on Scientific Computing, vol. 13, no. 1, pp. 425–448, 1992.

Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM Journal on Scientific Computing, vol. 7, no. 3, pp. 856–869, 1986.

S.-L. Wu, T.-Z. Huang, L. Li, and L.-L. Xiong, “Positive stable preconditioners for symmetric indefinite linear systems arising from Helmholtz equations,” Physics Letters A, vol. 373, no. 29, pp. 2401–2407, 2009.