The conjugate gradient method is one of the popular methods for solving linear system Ax = b. It is very suitable for large-scale sparse matrices arising form FD or FE approx-imation of boundary-value problems. When the n × n symmetric positive definite matrix has been preconditioned to make the calculations more effective, good results are obtain in only about √
n steps.
In this section, we introduce roughly the conjugate gradient method, and its complete derivation can be found in [5].
First, we define the quadratic function φ(x) = 1
2hAx, xi − hb, xi = 1
2xTAx − xTb, and let
h(α) = φ(x + αs).
Figure 3: Convergence of steepest descent.
also find that h(α) has a minimal value when
α = hs, b − Axi hs, Asi
Now, we denote x an approximate solution to Ax∗ = b, and vector s 6= 0 gives a search direction to improve the approximation. Let r = b − Ax be the residual vector associated with x and
α = hs, b − Axi
hs, Asi = hs, ri hs, Asi.
If r 6= 0 and if s and r are not orthogonal, then φ(x + αs) is smaller than φ(x) and x + αs is presumably closer to x∗ than x. This suggest the following method.
Let x0 be an initial approximation to x∗, and let s0 6= 0 be an initial search direction.
For k = 0, 1, 2, . . . , we compute
αk = hsk, b − Axki hsk, Aski , xk+1 = xk+ αksk
and choose a new search direction sk+1. The object is to make this selection so that the sequence of approximations {xk} converges rapidly to x∗.
Our direct idea is using −∇φ(x) as a search direction because it is the direction of greatest decrease in the value of φ(x). And it is just the direction of the residual r. The method that choose
sk+1 = rk+1 = b − Axk+1
is called the method of steepest descent. Unfortunately, the convergence rate of steepest descent is often very poor owing to repeated searches in the same directions (see Fig. 3).
Therefore, the search direction requires some modifications not negative gradient direction any more. The conjugate gradient method chooses search directions {s0, . . . , sn−1} so that they are A-orthogonal set; that is,
hsi, Asji = 0, if i 6= j.
8
Now, we use rk+1 to generate sk+1 by setting
sk+1 = rk+1+ βk+1sk. We want to choose βk+1 so that
hsk, Ask+1i = 0.
We can obtain
βk+1 = −hsk, Ark+1i hsk, Aski .
It can be shown that with this choice of βk+1, {s0, . . . , sk+1} is an A-orthogonal set. Then we can simplify
αk = hrk, rki hsk, Aski. Thus,
xk+1 = xk+ αksk.
To compute rk+1, we multiply by A and subtract b to obtain rk+1 = rk− αkAsk.
Then we can change
βk+1 = hrk+1, rk+1i hrk, rki .
Above derivation is the process of conjugate gradient method. The algorithm of conjugate gradient method is given in algorithm 4.
Algorithm 4 Conjugate Gradient Method x0 =initial guess
r0 = b − Ax0 s0 = r0
for k = 0, 1, 2, . . .
αk = rTkrk/sTkAsk {compute search parameter}
xk+1= xk+ αksk {update solution}
rk+1 = rk− αkAsk {compute new residual}
βk+1 = rTk+1rk+1/rTkrk
sk+1 = rk+1+ βk+1sk {compute new search direction}
end
conjugate gradient method can be accelerated by preconditioning. Preconditioning means that choose a matrix M for which systems of the form Mz = y are easily solved. And M−1 ≈ A−1 so that M−1A is relatively well-conditioned. In fact, we should apply con-jugate gradient to L−1AL−T instead of M−1A to preserve symmetry of matrix, where M = LLT. However, the algorithm can be suitably rearrange so that only M is used and the corresponding matrix L is not required explicitly. The algorithm of preconditioned conjugate gradient method is given in algorithm 5.
Algorithm 5 Preconditioned Conjugate Gradient Method x0 =initial guess
Choosing an appropriate preconditioner is very important. The choice of precondi-tioner depends on the usual trade-off between the gain in the convergence rate and the increased cost per iteration that results from applying the preconditioner. Many types of preconditioner can be find in [6]. We only introduce some of them.
• Jacobi : M is taken to be a diagonal matrix with diagonal entries equal to those of A. This is the easiest preconditioner since M−1 is also a diagonal matrix (can be regarded as vector when we compute) whose entries are reciprocal of the diag-onal entries of A. Although it only increases less storage and cost per iteration, it need more number of necessary iterations than the following incomplete Cholesky factorization.
• Incomplete Cholesky factorization: Ideally, we can factor A into LLTby Cholesky factorization, but this may incur unacceptable fill. We may instead compute an ap-proximate factorization A ≈ ˆL ˆLT that allows little or no fill (e. g. restricting the nonzero entries of ˆL to be in the same positions as those of the lower triangle of A) to save storage and CPU time, then use M = ˆL ˆLT as a preconditioner.
• Multigrid preconditioner : We choose directly M = A and find M−1r by multi-grid method. This is so-called multimulti-grid conjugate gradient method (MGCG). This
10
method is also researched widely in many papers. We will use MGCG in the later numerical experiments.
Up to now, We have introduced MG and PCG. These methods belong to iterative method.
In fact, fast direct methods have been developed for solving general poisson equation [9].
The following Concus and Golub method make a little modification to Helmholtz equation so that fast direct method can be used.
4 Concus and Golub Method
Concus and Golub propose an iterative scheme which uses fast direct solvers for the repeated solution of a Helmholtz problem [7]. Suppose original equation is
Ψu = −4u + g(x, y)u = f (x, y).
Concus and Golub provide a approach to utilize a modified form of the iterative procedure
−∆un+1= −∆un− τ (Ψun− f ), where τ is a parameter. We use the shifted iteration
(−4 + K)un+1= (−4 + K)un− τ (Ψun− f ).
The discrete form is given by
(−4h+ KI)Vn+1= (−4h+ KI)Vn− τ [(−4h+ G)Vn− F )],
where K is a parameter, V is approximation to exact u, −4his a matrix from operator −4 discretization with mesh space h, G is a diagonal matrix with elements Gij = g(ih, jh), F is a vector with elements Fij = f (ih, jh), and I is the identity matrix.
Now, we want to find spectral radius ρ for above iteration. We denote M ≡ −4h+ G operator, Φ is a vector, and νm,νM are the minimum and maximum eigenvalues of eigenvalue problem MΦ = ν(−4h+ KI)Φ. To obtain it,
ρ(I − τ [−4h+ KI]−1M) = max{|1 − τ νm|, |1 − τ νM|}.
To estimate νm and νM, we use the Rayleigh quotient for ν, ΦTMΦ
ΦT(−4h+ KI)Φ = 1 + ΦT(G − KI)Φ ΦT(−4h+ KI)Φ then
β − K β − K B − K B − K
where β and B are the minimal and maximal of function g(x, y),λm and λM are the smallest and largest eigenvalue of −∆h, and K is between β and B.
In this, the Lemma in [7] help us finding the optimal choice of τ
τ = 2
We will use this formula to observe the convergence rate of Concus and Golub method in the later numerical examples.
In an attempt to make the operator −4 + K on the left-hand side agree closely with Ψ, we usually set K the so-called min-max value,
1
2(min(g(x, y)) + max(g(x, y)))
so the optimal τ = 1. And then, the shifted iteration formula can be simplified (−4h+ KI)Vn+1= (KI − G)Vn+ F.
According to above discrete form, we give an initial guess first, and compute right-hand side. Now, we can solve “Poisson equation” by FPS to gain the approximation and then update right-hand side again. Repeat this action until it reaches stopping criterion.
In fact, K can also be optimized for higher rates of convergence. The efficiency of Concus and Golub’s method can be increased by extending its formulation to accommo-date the use of a parameter K which is a one-dimensional function instead of a constant [7]. Unfortunately, for a special solution, that would have required a prohibitively large computational time since there is no fast Helmholtz solver available for variable coefficient.
In the past it has been demonstrated that the number of necessary iterations can vary dramatically, depending on the function g(x, y) which has a critical role on the rate of convergence. The smoother g(x, y) is, the faster rate of convergence. We will give a example in the next section.
5 Numerical Examples
We will solve Helmholtz equations with Dirichlet boundary condition
−4u(x, y) + g(x, y)u = f (x, y), (x, y) ∈ Ω (3) u(x, y) = p(x, y), (x, y) ∈ ∂Ω
on domain Ω = [−1, 1] × [−1, 1].
12
First, we use five-point finite difference formula to discretize e.q. (3). Assume domain is partitioned into n × n subdomain by introducing the grid points
xi = ih, yj = jh, 0 ≤ i, j ≤ n
Then, we can obtain discrete component form
−1
h2 (vi−1,j + vi+1,j+ vi,j−1+ vi,j+1− 4vij) + Gijvij = Fij, and we can use iterative method to improve approximation vij according to
vij = 1
We may also represent system in matrix form as Av = b.
A is an (n − 1) × (n − 1) block matrix, T is an (n − 1) × (n − 1) matrix and I is the identity matrix of order n − 1.
A = 1
And the unknown vector v is defined by
And the right-hand side vector b is defined by
b = Our stopping criterion is krkh
kbkh < 10−8, where r = b − Av is residual, and k • kh means discrete L2 norm.
We use the following methods for solving later examples.
• Concus and Golub method: Choose K equal to min-max value, and use “GEN-BUN” routine which is in “FISHPACK” as our fast poisson solver.
• V-cycle method: Use RBGS relaxation with 2 sweeps before the correction step and 1 sweep after the correction step. And use full weighting restriction and linear interpolation operators.
• ICCG method: Use Incomplete Cholesky factorization as a preconditioner of con-jugate gradient method.
• MGCG method: Use V-cycle method as a preconditioner of conjugate gradient method.
14
We choose two coefficient functions contrary, go(x, y) is a large amplitude function. The maximum and minimum of go(x, y) are 160 and −32 respectively.
We also choose two functions
us(x, y) = sin(πx) + cos(πy) and
uo(x, y) = (y2− 1)(y − 1)(e5ysin(4πx) + e−5ycos(4πx))
as our exact solutions. The function us(x, y) is smoother than uo(x, y). This also means that the approximation of us(x, y) is more accurate than the approximation of uo(x, y) under the same conditions.
Now, we use these functions gs(x, y), go(x, y), us(x, y) and uo(x, y) to result four examples, and show numerical results including CPU time (sec.) and necessary iterations to achieve stopping criterion on different grid numbers N = 1282, 2562, 5122, 10242.
Example 1.1: gs(x, y) uo(x, y)
Iterations Concus and Golub ICCG MGCG V-cycle
N = 1282 5 117 7 7
N = 2562 5 218 7 7
N = 5122 5 428 7 6
N = 10242 5 796 7 6
CPU time Concus and Golub ICCG MGCG V-cycle
N = 1282 0.3 0.7 0.6 0.6
N = 2562 1.1 4.6 2.3 2.2
N = 5122 4.8 34.0 9.2 7.8
N = 10242 23.0 264.0 36.6 31.0
In this example, we can find that ICCG is not a good method since it need longer CPU time to achieve stopping criterion. So incomplete Cholesky factorization is not a good
Cholesky factorization (MIC) for construction of the preconditioning matrix, see [6]. We also find that the number of iterations of ICCG increases when N becomes larger.
The other methods like Concus and Golub, MGCG, V-cycle converge with very few iterations. And their number of iterations is almost independent of N. Therefore, the CPU time needed for these methods reach stopping criterion is much less than ICCG when N is very large.
Up to now, Concus and Golub method seem to converge rapidly. It is even better than multigrid method. Let us continue seeing the next example.
Example 1.2: go(x, y) uo(x, y)
Iterations Concus and Golub ICCG MGCG V-cycle
N = 1282 271 144 9 11
N = 2562 267 282 9 10
N = 5122 266 533 9 9
N = 10242 266 1006 9 8
CPU time Concus and Golub ICCG MGCG V-cycle
N = 1282 9.0 0.73 0.7 0.8
N = 2562 36.0 5.2 2.8 3.0
N = 5122 158.0 39.0 11.6 10.8
N = 10242 784.0 305.0 46.0 40.0
This example uses a high amplitude function go(x, y) as coefficient. Clearly, the iterations and CPU time of all methods increase. But we note that Concus and Golub method results in very very slow convergence. Its CPU time needed is even longer than ICCG.
This phenomenon can be explained by below inequality introduced in the last section.
ρ ≤ M − m
2λ + M + m < 1,
where M and m denote the maximum and minimum of g(x, y) respectively. λ is the smallest eigenvalue of the discrete matrix of Laplace operator, and ρ is the spectral radius of iterative matrix of Concus and Golub method. In example 1.1
M − m
2λ + M + m ≈ 0.08, but in example 1.2
M − m
2λ + M + m ≈ 0.93.
That is why Concus and Golub method need so many iterations to reach stopping criterion here. In next two examples 2.1,2.2 we replace uo(x, y) in examples 1.1,1.2 with us(x, y).
Example 2.1: gs(x, y) us(x, y)
16
Iterations Concus and Golub ICCG MGCG V-cycle
N = 1282 9 124 7 7
N = 2562 9 232 7 6
N = 5122 9 442 7 6
N = 10242 9 869 6 6
CPU time Concus and Golub ICCG MGCG V-cycle
N = 1282 0.3 0.6 0.5 0.5
N = 2562 1.3 4.1 1.9 1.6
N = 5122 5.6 30.6 7.8 6.4
N = 10242 28.3 257.0 27.3 25.6
Example 2.2: go(x, y) us(x, y)
Iterations Concus and Golub ICCG MGCG V-cycle
N = 1282 394 144 9 12
N = 2562 393 282 9 10
N = 5122 392 554 9 10
N = 10242 392 1078 9 8
CPU time Concus and Golub ICCG MGCG V-cycle
N = 1282 12.0 0.63 0.6 0.8
N = 2562 53.0 4.9 2.6 2.6
N = 5122 228.0 38.0 10.2 10.6
N = 10242 1153.0 317.0 41.0 35.0
After observing results of these two examples, we get the same conclusion discussed before. But we can find another phenomenon by comparing with example 1.1 and example 1.2. The number of iterations and CPU time needed of Concus and Golub method and ICCG increase, but MGCG and V-cycle almost remain the same number of iterations even if us(x, y) may cause smoother error in these examples. That is because multigrid method can reduce smooth error efficiently.
In this thesis, we do not compare MGCG and V-cycle (ordinary multigrid) methods.
Which is the better method? What is the difference of behaviors between the MCCG and multigrid methods? We can find answer in [10]. In [10], it gives a special example. The matrix produced by the discretization of this example has scattered eigenvalues. This scattered eigenvalues prevent the ordinary multigrid method from converging rapidly.
We can observe the eigenvalue’s distribution after multigrid preconditioning from [10].
Almost all eigenvalues are clustered, and a few eigenvalues are scattered. The conjugate gradient method hides the defect of the multigrid method. Therefore, the MGCG method
6 Conclusion
After testing above examples, we know that multigrid method suits to be used for solving general nonseparable elliptic equations. Since no matter g(x, y) is smooth or oscillatory and smooth or oscillatory error resulted, multigrid method always converges in very few iterations. This means that multigrid method convergence stably. We also know an important advantage of multigrid method from our nu- merical examples. That is the number of necessary iterations is independent of mesh size.
There is another advantage of multigrid method. It is easy to be parallelized. We do not introduce the parallelism of multigrid method in this thesis. The content of parallel multigrid can be studied from [11]. These above reasons explain why we use multigrid method for solving Helmholtz equation here.
18
References
[1] William L. Briggs Van Emden Henson Steve F. McCormick, A multigrid toturial 2nd ed. 2000.
[2] A. BRANDT, Multigrid Techniques: 1984 Guide with Applications to Fluid Dynam-ics, GMD-Studien Nr. 85, Gesellschaft f¨ur Mathematik und Datenverarbeitung, St.
Augustin, Bonn, 1984.
[3] Michael T. Heath, Scientific computing :an introductory survey 2nd ed. 2002.
[4] W. M. Pickering and P. J. Harley , Numerical solution of non-separable elliptic equa-tions by iterative application of FFT methods. Int. Jnl. Computer Math. 55, 211-222, 1995.
[5] Richard L. Burden, J. Douglas Faires, Numerical analysis 7th ed. 2001.
[6] Yousef Saad, Iterative methods for sparse linear systems, 1996.
[7] Paul Concus and Gene H. Golub, Use of fast direct met hods for the efficient nu-merical solution of nonseparable elliptic equations, SIAM J. Numer. Anal. 10 pp.
1103-1120, 1973.
[8] Dimitropoulos, C.D. and Beris, A.N. An efficient and robust spectral solver for non-separable elliptic equations J. Comp. Phys. 133 pp. 186-191, 1997.
[9] B. Buzbee, G. Golub, and C. Nielson, On direct methods for solving Poissons equa-tion, SIAM J. Numer. Anal. 7 pp. 627-656, 1970.
[10] Osamu Tatebe, The multigrid preconditioned conjugate gradient method [11] Jim E. Jones, A Parallel Multigrid tutorial.