# New Preconditioners with Two Variable Relaxation Parameters for the Discretized Time-Harmonic Maxwell Equations in Mixed Form

Authors: Yuping Zeng; Chenliang Li
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

We provide new preconditioners with two variable relaxation parameters for the saddle point linear systems arising from finite element discretization of time-harmonic Maxwell equations in mixed form. The new preconditioners are of block-triangular forms and Schur complement-free. They are extensions of the results in Cheng et al., 2009, Grief and Schötzau, 2007, and Huang et al., 2009. Theoretical analysis shows that all eigenvalues of the preconditioned matrices are tightly clustered, and numerical tests confirm our analysis.


---

## Body

## 1. Introduction

We consider the preconditioning techniques for solving the saddle point linear systems arising from finite element discretization of the following time-harmonic Maxwell equations in mixed form [1–5]: find the vector field u and the Lagrangian multiplier p such that
(1.1)∇×∇×u-k2u+∇p=finΩ,∇·u=0inΩ,u×n→=0on∂Ω,p=0on∂Ω.
Here, Ω⊂ℝ2 is a simply connected polyhedron domain with a connected boundary ∂Ω, and n→ denotes the outward unit normal on ∂Ω. The datum f is a given source (not necessarily divergence-free), and the wave number k2=ω2ϵμ, where ω≥0 is the frequency, and ϵ and μ are positive permittivity and permeability parameters, respectively.In recent years, there have been many techniques for solving Maxwell equations, such as the geometry multigrid methods [6–8], algebraic multigrid methods [9], domain decomposition methods [4, 10–13], Nodal auxiliary space preconditioning methods [14], and the solution methods to the corresponding saddle-point linear systems [2, 3, 15]. We can also use Uzawa-type iterative methods [16, 17] and preconditioned Krylov subspace methods [18–24] to solve the saddle-point linear systems. Based on the previous works in [2, 3, 15], we will further study solution methods for the saddle-point linear systems in this paper.Using Nédélec elements of the first kind [25–27] for the approximation of the vector field and standard nodal elements for the Lagrangian multiplier yields the following saddle-point linear system:
(1.2)𝒜x≡(A-k2MBTB0)(up)=(g0)≡b,
where u∈Rn and p∈Rm are finite arrays, and g∈Rn is a load vector associated with f. The matrix A∈Rn×n is symmetric positive semidefinite with nullity m and corresponds to the curl-curl operator; B∈Rm×n is a discrete divergence operator with full-row rank, and M∈Rn×n is a vector mass matrix.For convenience, we denote the standard Euclidean inner product of vectors by〈·,·〉  and the null space of a matrix by null (·). For a given positive (semi)definite matrix W and a vector x, we define the (semi)norm:
(1.3)|x|W=〈Wx,x〉.The matricesA and B have the following stability properties [3]. Let 〈Au,u〉=|u|A2. Then there exists an α, 0<α<1, such that
(1.4)|u|A2≥α-|u|M2,u∈null(B),
where α-=α/(1-α). Matrix B satisfies the discrete inf-sup condition:
(1.5)inf0≠q∈Rmsup0≠v∈null(A)〈Bv,q〉|v|M|q|L≥β>0,
where the inf-sup constant β>0 is only dependent on the domain Ω.If the wave numberk2>0, then the (1,1) block of (1.2) is indefinite. For difficulty and corresponding solution methods of this problem, we refer to [18, 28]. Recently, by using the spectral equivalent properties similar to [4], Grief and Schötzau [3] construct the block-diagonal preconditioner:
(1.6)ℳk=(A-k2M+BTL-1B00L),
where L∈Rm×m is the discrete Laplace operator introduced in [3], k2<1, and ℳk is a symmetric positive definite block-diagonal matrix. As L is augmentation-free and Schur complement-free, this approach overcomes the difficulty in forming the Schur complement in general. However, the computational work of BTL-1B may be too large. Using the fact that the matrices A+BTL-1B and A+M are spectrally equivalent, [3] considers the following preconditioner:
(1.7)ℳ^k=(A+(1-k2)M00L),
and shows that the eigenvalues of the preconditioned matrix are tightly clustered.Based on the work of Grief and Schötzau [3], [2] gives block-triangular Schur complement-free preconditioners for the linear system (1.2). And it is shown that all eigenvalues of the proposed block-triangular preconditioning matrices are more tightly clustered. Compared with the restriction k2<1 in [3], [2] considers the general case k2∈R+. Furthermore, [15] provides block-triangular preconditioners when k2=0 with two variable relaxation parameters.Based on the previous work [2, 3, 15], mentioned above, in this paper we are devoted to give new preconditioners with two scaling parameters. The new block triangular preconditioners in the general case k2∈R+ contain the preconditioners discussed in [2]. Theoretical analysis shows that all eigenvalues of the preconditioned matrices are tightly clustered. Numerical experiments demonstrate efficiency of the new method and show that preconditioner ℳ^k,η,ε is more efficient than ℳ^k,t.The remainder of the paper is as follows. In Section2, we establish new block-triangular preconditioners for the linear systems (1.2) in the general case k2∈R+, and then the corresponding spectral analysis is presented. In Section 3, we provide numerical examples to examine our analysis. Finally, some conclusions are drawn in Section 4.


## 2. New Block-Triangular Preconditioners for Anyk2

We consider the saddle-point linear system (1.2) arising from the discretized time-harmonic Maxwell equations in mixed form (1.1) and assume that k2 is not an eigenvalue and k2∈R+.Grief and Schötzau [3] provide the block-diagonal Schur complement-free preconditioner ℳk as in (1.6). Using the fact that the matrices A+BTL-1B and A+M are spectrally equivalent, the argumentation-free and Schur complement-free preconditioner ℳ^k is defined in (1.7). Spectral analysis shows that the eigenvalues of the preconditioned saddle-point matrices ℳk-1𝒜 and ℳ^k-1𝒜 are strongly clustered when k2 is small.Reference [2] provides the block-triangular Schur complement-free preconditioners for the linear system (1.2). In particular, they considered preconditioning matrices for the general case k2∈R+ with
(2.1)ℳk,t=(A-k2M+tBTL-1BtBT01-ttL),ℳ^k,t=(A+(t-k2)MtBT01-ttL),
where 1≠t>k2.Fork2=0, [15] provides the block-triangular preconditioner for linear system (1.2):
(2.2)𝒫η,ε=(A+ηBTL-1B(1-ηε)BT0εL),
where η>0 and ε>0.Based on the works in [2, 3, 15], we provide the following new block-triangular Schur complement-free preconditioners ℳk,η,ε and ℳ^k,η,ε:
(2.3)ℳk,η,ε=(A-k2M+ηBTL-1B(1-ηε)BT0εL),(2.4)ℳ^k,η,ε=(A+(η-k2)M(1-ηε)BT0εL),
where η>k2 and ε≠0 are scaling parameters. It is interesting to note that when parameters η=t and ε=(1-t)/t, ℳk,η,ε and ℳ^k,η,ε apparently reduce to ℳk,t and ℳ^k,t, respectively. We also see that when k2=0, the preconditioner ℳk,η,ε in (2.3) (ε≠0) is different from 𝒫η,ε in (2.2) (ε>0).We stress thatℳk,η,ε is not the preconditioner we eventually use in actual computation. It is only introduced to lay theoretical basis and motivation for the preconditioner ℳ^k,η,ε in (2.4), which we will use in practice. We note that the (1,1) block A-k2M+ηBTL-1B in ℳk,η,ε is symmetric positive definite for k is sufficiently small [3]. But this is not true when k is large enough. However, this situation may not appear in the actual preconditioner ℳ^k,η,ε. The (1,1) block A+(η-k2)M in ℳ^k,η,ε is always symmetric positive definite when η>k2. In this paper, we will apply the BiCGSTAB with the preconditioner ℳ^k,η,ε as an outer solver for the saddle-point system (1.2). Then, the overall computational cost of solution procedure relies on how to efficiently solve the linear systems A+τM(τ=η-k2) and L, which are called by inner solvers. For the linear system L arising from a standard scalar elliptic problem, many efficient solution methods exist. On the other hand, for solving the linear system A+τM, we refer to [6, 8, 9, 14], and some detailed numerical examples are provided in [29].For the spectral analysis, we recall some results which are contained in the following lemma.Lemma 2.1 (see [3]).
The following relations hold:(i)
ℛn=null(A)⊕null(B);(ii)
〈MuA,uB〉=0foranyuA∈null(A)andanyuB∈null(B);(iii)
〈BTL-1BuA,uA〉=〈MuA,uA〉foranyuA∈null(A).Theorem 2.2.
Let𝒜 be the saddle-point matrix in (1.2). Then the matrix ℳk,η,ε-1𝒜 has two distinct eigenvalues, which are given by
(2.5)λ1=1,λ2=-1ε(η-k2),
with the algebraic multiplicities n and m, respectively.Proof.
Suppose thatλ is an eigenvalue of ℳk,η,ε-1𝒜, whose eigenvector is (vq). Then the corresponding eigenvalue problem is
(2.6)(A-k2MBTB0)(vq)=λ(A-k2M+ηBTL-1B(1-ηε)BT0εL)(vq).
From the second row we can obtainq=(1/λε)L-1Bv. By substituting it into the first row we have
(2.7)(1-λ)[λ(A-k2M)v+(1ε+λη)BTL-1Bv]=0.
It is straightforward to see that any vectorv∈ℛn satisfies (2.7) with λ=1, so λ=1 is an eigenvalue of ℳk,η,ε-1𝒜. By the similar technique of linear independence considerations from [3], we can demonstrate that the eigenvalue λ=1 has algebraic multiplicity n.
Since there arem linearly independent null vectors of A, by Lemma 2.1,
(2.8)v=vA+vB(vA≠0),vA∈null(A),vB∈null(B).
By Lemma 2.1 (ii) and (iii), and using the inner product in (2.7) with vA, we have
(2.9)(1-λ)[1ε+λ(η-k2)]|vA|M2=0.
Since vA≠0, from (2.9) we can obtain that λ=-1/ε(η-k2) is another eigenvalue of ℳk,η,ε-1𝒜, and we claim that the eigenvalue λ=-1/ε(η-k2) has algebraic multiplicity m.Corollary 2.3.
Let-1/ε=η-k2. Then the corresponding preconditioned matrix ℳk,η,ε-1𝒜 has only one eigenvalue λ=1 of algebraic multiplicity n+m.Proof.
From Theorem2.2, we can easily obtain the corresponding conclusion.Remark 2.4.
From Theorem2.2, we demonstrate that the preconditioned matrix ℳk,η,ε-1𝒜 has precisely two distinct eigenvalues. Then if Krylov subspace methods are used to solve (1.2) with ℳk,η,ε as a preconditioner, the iteration will require merely two steps if round-off errors are ignored [30]. And from Corollary 2.3, for any η, we can find a number ε which makes the preconditioned matrix ℳk,η,ε-1𝒜 have only one eigenvalue. Therefore, we can demonstrate that our preconditioners are more efficient than the block-triangular preconditioner proposed in [2].Remark 2.5.
From (2.3) we know that if ηε=1, the new preconditioner reduces to the diagonal preconditioner ℳk,η:
(2.10)ℳk,η=(A-k2M+ηBTL-1B001ηL).
Then we can use MINRES to solve the linear system (1.2).Theorem 2.6.
Let𝒜 be the saddle-point matrix in (1.2). Then
(2.11)λ1=1,λ2=-1ε(η-k2)
are the eigenvalues of ℳ^k,η,ε-1𝒜, having algebraic multiplicity m. The rest of the eigenvalues satisfy
(2.12)α--k2α-+η-k2≤λ<1,
where α- is defined as in (1.4).Proof.
Suppose thatλ is an eigenvalue of ℳ^k,η,ε-1𝒜, whose eigenvector is (vq). Then the corresponding eigenvalue problem is
(2.13)(A-k2MBTB0)(vq)=λ(A+(η-k2)M(1-ηε)BT0εL)(vq).
From the second row we can obtain q=(1/ελ)(L-1Bv). By substituting it into the first row we have
(2.14)λ(A-k2M)v+1εBTL-1Bv=λ2(A+(η-k2)M)v+λε(1-ηε)BTL-1Bv.
Consider the m linearly independent null vectors of A, by Lemma 2.1 (i),
(2.15)v=vA+vB(vA≠0),
where vA∈ null(A) and vB∈ null(B). By Lemma 2.1 (ii) and (iii), and taking the inner product in (2.14) with vA, we obtain
(2.16)(1-λ)[1ε+λ(η-k2)]|vA|M2=0.
Since |vA|M≠0, λ1=1 and λ2=-1/(ε(η-k2)) are two eigenvalues of ℳ^k,η,ε-1𝒜 and by the similar technique of linear independence considerations from [3], we claim that each eigenvalue has algebraic multiplicity m.
For the rest of eigenvectors we havevB≠0. Noting that
(2.17)〈BTL-1BvA,vB〉=〈L-1BvA,BvB〉=0,
by Lemma 2.1 (ii) and by taking the inner product in (2.14) with vB and using (2.17), we obtain
(2.18)(1-λ)|vB|A2=(λ(η-k2)+k2)|vB|M2.
It is impossible to have λ=1, since (2.18) leads to |vB|M=0, which contradicts with vB≠0. We cannot have λ>1, since the left-hand side is negative but the right-hand side is positive (because we assume η>k2). Thus, we claim that λ<1.
From (1.4), we recall that for any u∈null(B), |u|A2≥α-|u|M2 with α-=α/(1-α). Applying this to (2.18), we have (λ(η-k2)+k2)/(1-λ)≥α-. Since η>k2>0, we have α-+η-k2>0 and
(2.19)α--k2α-+η-k2≤λ<1.Corollary 2.7.
Let-1/ε=η-k2. Then the corresponding preconditioned matrix ℳ^k,η,ε𝒜 has only one eigenvalue λ=1 with algebraic multiplicity 2m. The remaining eigenvalues satisfy (2.12).Remark 2.8.
From (2.4) we know that when ηε=1, the new preconditioner reduces to the diagonal preconditioner ℳ^k,η:
(2.20)ℳ^k,η=(A+(η-k2)M001ηL).
Then we can use MINRES to solve the linear system (1.2).Remark 2.9.
From the proof of Theorem2.6, we easily see that the new preconditioner ℳ^k,η,ε is also efficient for k2=0. Then from (2.12), we conclude that if -1/ε=η-k2 and k2=0, then the closer η is to 0 and the closer (α--k2)/(α-+η-k2)=α-/(α-+η) is to 1; that is, the preconditioned matrix has more tightly clustered eigenvalues. For a general case of k2∈ℝ+, we can only obtain similar results when α->k2. The following numerical experiments show that the closer η is to k2, the less iteration counts we have used for a fixed k2∈ℝ+. However, choosing η-k2 too small may result in too large ε, then result in ill-conditioning of ℳ^k,η,ε. So we choose η-k2 to be moderate size in practice.


## 3. Numerical Experiments

The test problem is a two-dimensional time-harmonic Maxwell equations in mixed form (1.1) in a square domain Ω=(0<x<1;0<y<1). We set the right-hand side function so that the exact solution is given by
(3.1)u(x,y)=(u1(x,y)u2(x,y))=(y(1-y)x(1-x))
and p≡0.We consider five uniformly refined meshes, which are constructed by subsequently splitting each triangle into four triangles by joining the midpoints of the edges of the triangle. Two of five mesh grids are depicted in Figures1 and 2. The lowest order elements are used to discretize equations. The matrix sizes on different meshes are given in Table 1.Table 1
Values of matrix size of the linear system for five meshes.
Mesh
n
m
n + m
8
×
8
176
49
225
16
×
16
736
225
961
32
×
32
3008
961
3969
64
×
64
12160
3969
16129
128
×
128
48896
16129
65025Figure 1
8
×
8 mesh.Figure 2
16
×
16 mesh.Our numerical experiments were performed using MATLAB. The machine is a PC-Intel (R), Pentium(R)Dual CPU E2200 2.20 GHz, 1.00 G of RAM. The purpose of our experiments is to investigate the convergence behavior of preconditioned BiCGSTAB by choosing different parametersη and and ε in the preconditioner ℳ^k,η,ε. Thus, we apply exact inner solver, and the outer iteration is used as a zero initial guess and stopped when ∥r(k)(=b-𝒜x(k))∥2/∥r(0)(=b)∥2≤5×10-10.From Theorem2.6 and Corollary 2.7 we know that the preconditioned matrix ℳ^k,η,ε-1𝒜 has one eigenvalue λ=1, and the remaining eigenvalues are satisfying (2.12). Figure 3 depicts the eigenvalues of the preconditioned matrix ℳ^k,η,ε-1𝒜 with k2=2, where η and ε satisfy -1/ε=η-k2, and m+n=961. From it we observe that the eigenvalue distribution of preconditioned matrix ℳ^k,η,ε-1𝒜 with η=k2+0.1=2.1 denoted by solid line is more tightly clustered than with η=k2+30=32 denoted by dotted line. From Remark 2.9 we know that the new preconditioner ℳ^k,η,ε is also efficient for k2=0. Figures 4, 5, 6, and 7 show the eigenvalue distribution of the preconditioned matrix ℳ^k,η,ε-1𝒜 for different η with k2=0, where η and ε satisfy -1/ε=η-k2, and m+n=961. From Figures 4–7 we know that the closer the parameter η is to 0, the more tightly clustered the eigenvalues of the preconditioned matrix will be.Figure 3
The eigenvalue distribution of the preconditioned matrixℳ^k,η,ε-1𝒜 with k2=2, η=2.1 and η=32, respectively, and m+n=961. The case η=2.1 is indicated by the solid line while the case η=32 is indicated by the dotted line.Figure 4
The eigenvalue distribution of the preconditioned matrixℳ^k,η,ε-1𝒜 with k2=0, η=5 and m+n=961.Figure 5
The eigenvalue distribution of the preconditioned matrixℳ^k,η,ε-1𝒜 with k2=0, η=0.1 and m+n=961.Figure 6
The eigenvalue distribution of the preconditioned matrixℳ^k,η,ε-1𝒜 with k2=0, η=0.01 and m+n=961.Figure 7
The eigenvalue distribution of the preconditioned matrixℳ^k,η,ε-1𝒜 with k2=0, η=0.001 and m+n=961.Table2 shows the outer iteration counts for different k2 and η, applying BiCGSTAB with the block-triangular preconditioner, where η and ε satisfy -1/ε=η-k2, and m+n=16129. The iteration counts are denoted by Iter. We observe that for a fixed k2, the closer η is to k2, the less iteration counts are produced. For comparison, we also give the outer iteration counts for η=t*=1+(k2+1+k4)/2. We refer the definition of t* to [2]. It shows that the preconditioner ℳ^k,η,ε is more efficient than ℳ^k,t.Table 2
Iteration counts for differentk2 and η, using BiCGSTAB with the preconditioner M^k,η,ε, and m+n=16129. The iteration was stopped once ∥r(k)∥/∥r(0)∥≤5×10-10.
η
k
2
0.1
0.2
0.5
1
2
3
4
5
k
2
+
0.001
Iter
2
2
2.5
2.5
3
3.5
4
4.5
k
2
+
0.01
Iter
2
2
2.5
2.5
3
3.5
4
4.5
k
2
+
0.1
Iter
2
2
2.5
2.5
3
3.5
4
4.5
k
2
+
0.25
Iter
2.5
2.5
2.5
2.5
3.5
3.5
4
4.5
k
2
+
0.5
Iter
2.5
2.5
2.5
2.5
3.5
4
4
4.5
k
2
+
1
Iter
2.5
2.5
2.5
3
3.5
4
4
4.5
k
2
+
1.5
Iter
2.5
2.5
3
3
3.5
4
4.5
4.5
k
2
+
2
Iter
3
3
3
3.5
3.5
4
4.5
4.5
k
2
+
5
Iter
3.5
3.5
3.5
4
4.5
4.5
4.5
5
k
2
+
15
Iter
4.5
5
5
5
5
5.5
5.5
5.5
k
2
+
30
Iter
5.5
5.5
5.5
6
6
6
6
6
t
*
=
1
+
(
(
k
2
+
1
+
k
4
)
/
2
)
Iter
3
3
3
3
3.5
4
4
4.5Tables3 and 4 show the outer iteration numbers for different meshes, applying BiCGSTAB with the preconditioner ℳ^k,η,ε, where η are set to be η-k2=0.1 and η-k2=6, and -1/ε=η-k2. We observe that the outer iteration numbers of the preconditioned BiCGSTAB are hardly sensitive to the changes in the mesh size.Table 3
Iteration counts for different meshes, using BiCGSTAB with the preconditionerM^k,η,ε satisfying η-k2=0.1 and -1/ε=η-k2. The iteration was stopped once ∥r(k)∥/∥r(0)∥≤5×10-10.
Mesh
k
2
=
0
k
2
=
0.25
k
2
=
0.5
k
2
=
1
k
2
=
3
k
2
=
4
k
2
=
6
k
2
=
10
8
×
8
2
2.5
2.5
2.5
3.5
4
4.5
5.5
16
×
16
2
2.5
2.5
2.5
3.5
4
4.5
5
32
×
32
2
2.5
2.5
2.5
3.5
4
4.5
5.5
64
×
64
2
2.5
2.5
2.5
3.5
4
5
5.5
128
×
128
2
2.5
3
3
4
4.5
5
6Table 4
Iteration counts for different meshes, using BiCGSTAB with the preconditionerM^k,η,ε satisfying η-k2=6 and -1/ε=η-k2. The iteration was stopped once ∥r(k)∥/∥r(0)∥≤5×10-10.
Mesh
k
2
=
0
k
2
=
0.25
k
2
=
0.5
k
2
=
1
k
2
=
3
k
2
=
4
k
2
=
6
k
2
=
10
8
×
8
3.5
3.5
3.5
4
4.5
4.5
5
6.5
16
×
16
3.5
4
4
4
4.5
4.5
5.5
6.5
32
×
32
3.5
4
4
4
4.5
4.5
5.5
6
64
×
64
3.5
4
4
4
4.5
5
5.5
6
128
×
128
3.5
4
4.5
4.5
5
5
6
6.5


## 4.   Conclusions

We have investigated the use of new block-triangular preconditioners with two variable relaxation parameters for solving the mixed formulation of the time-harmonic Maxwell equations. Our results are extensions of the work in [2, 3, 15]. The preconditioned matrices are demonstrated to have clustering eigenvalues. We have shown experimentally that the outer iteration numbers of BiCGSTAB with the new preconditioner are hardly any sensitive to the changes in the mesh size.



---

*Source: 617390-2012-08-21.xml*