110
5. Solution of Linear Equation Systems
Golub and van Loan, 1990). In this description, pk is the residual at the kth
iteration, pk is the kth search direction, z k is an auxiliary vector and a k and
Dk are parameters used in constructing the new solution, residual, and search
direction. The algorithm can be summarized as follows:
0 Initialize by setting: k = 0, 4' = 4in, p0 = Q - A4in, p0 = 0, so = 10 30
0 Advance the counter: k = k + 1
0 Solve the system: M Z ~
= pk-'
0 Calculate: sk = pk--' - z k
p = sklsk-1
p k = zz" + p"k-1
ak = sk/(pk ' Apk)
4 k = 4 k - 1 + a k p k
pk = & I - ak Apk
0 Repeat until convergence.
This algorithm involves solving a system of linear equations at the first step.
The matrix involved is M = C-' where C is the pre-conditioning matrix,
which is in fact never actually constructed. For the method to be efficient,
M must be easy to invert. The choice of M used most often is incomplete
Cholesky factorization of A but in tests it was found that if M = LU where
L and U are the factors used in Stone's SIP method, faster convergence is
obtained. Examples will be presented below.
5.3.7 Biconjugate Gradients and CGSTAB
The conjugate gradient method given above is applicable only to symmetric systems; the matrices obtained by discretizing the Poisson equation are
often symmetric (examples are the heat conduction equation and pressure
or pressure-correction equations to be introduced in Chap. 7). To apply the
method to systems of equations that are not necessary symmetric (for example, any convection/diffusion equation), we need to convert an asymmetric
problem to a symmetric one. There are a couple of ways of doing this of which
the following is perhaps the simplest. Consider the system:
This system can be decomposed into two subsystems. The first is the original system; the second involves the transpose matrix and is irrelevant. (If
there were a need to do so, one could solve a system of equations involving
the transpose matrix at little extra cost.) When the pre-conditioned conjugate gradient method is applied to this system, the following method, called
biconjugate gradients, results:
0 Initialize by setting: k = 0, 4' = 4 i n l p0 = Q - A4in, = Q - AT4in,
= 8 = 0, so = 1o30
5. Solution of Linear Equation Systems
Golub and van Loan, 1990). In this description, pk is the residual at the kth
iteration, pk is the kth search direction, z k is an auxiliary vector and a k and
Dk are parameters used in constructing the new solution, residual, and search
direction. The algorithm can be summarized as follows:
0 Initialize by setting: k = 0, 4' = 4in, p0 = Q - A4in, p0 = 0, so = 10 30
0 Advance the counter: k = k + 1
0 Solve the system: M Z ~
= pk-'
0 Calculate: sk = pk--' - z k
p = sklsk-1
p k = zz" + p"k-1
ak = sk/(pk ' Apk)
4 k = 4 k - 1 + a k p k
pk = & I - ak Apk
0 Repeat until convergence.
This algorithm involves solving a system of linear equations at the first step.
The matrix involved is M = C-' where C is the pre-conditioning matrix,
which is in fact never actually constructed. For the method to be efficient,
M must be easy to invert. The choice of M used most often is incomplete
Cholesky factorization of A but in tests it was found that if M = LU where
L and U are the factors used in Stone's SIP method, faster convergence is
obtained. Examples will be presented below.
5.3.7 Biconjugate Gradients and CGSTAB
The conjugate gradient method given above is applicable only to symmetric systems; the matrices obtained by discretizing the Poisson equation are
often symmetric (examples are the heat conduction equation and pressure
or pressure-correction equations to be introduced in Chap. 7). To apply the
method to systems of equations that are not necessary symmetric (for example, any convection/diffusion equation), we need to convert an asymmetric
problem to a symmetric one. There are a couple of ways of doing this of which
the following is perhaps the simplest. Consider the system:
This system can be decomposed into two subsystems. The first is the original system; the second involves the transpose matrix and is irrelevant. (If
there were a need to do so, one could solve a system of equations involving
the transpose matrix at little extra cost.) When the pre-conditioned conjugate gradient method is applied to this system, the following method, called
biconjugate gradients, results:
0 Initialize by setting: k = 0, 4' = 4 i n l p0 = Q - A4in, = Q - AT4in,
= 8 = 0, so = 1o30