Page 301 - Numerical Methods for Chemical Engineering

290 6 Boundary value problems

because, being generated by elimination, fill-in causes them to have many more nonzero
elements than the original matrix A. However, it still may be possible to find approximate
factors PA ≈ M 1 M 2 such that M −1 PAM −1 has a condition number closer to 1 than the
1 2
original system PA. Then, the conjugate gradient or GMRES search directions point more
closely towards the solution and require fewer iterations to reduce the residual norm to near
zero. PA ≈ M, M = M 1 M 2 is then said to serve as a preconditioner for A. The update to
the original system at each iteration is obtained by backward substitution of
[k+1] [k] [k]
M 2 x = ˆx + ˆp (6.166)
T
For a positive-definite matrix A, the existence of the Cholesky factorization A = R R,
T
allows us to use a preconditioner A ≈ M M 2 , such that the transformed system
2
T(−1) −1 T
M 2 AM 2 ˆ x = c M c = b (6.167)
2
is also positive-definite, and thus can be solved by conjugate gradients.

How should we choose the preconditioner?

There are many possible approaches, and sometimes one uses a preconditioner specifically
tailored to the PDE being solved. Here, we consider only general approaches. The simplest
preconditioner is the Jacobi preconditioner, for which M is merely the diagonal part of A,
M = diag(diag(A)). More effective preconditioners are obtained from incomplete factor-
izations of A, using an incomplete Cholesky factorization if A is positive-definite (and we
are using pcg)oran incomplete LU factorization if A is not (and we are using gmres or
bicgstab).

What do we mean by an incomplete factorization?
T
Consider the Cholesky factorization of a positive-definite matrix, A = R R, where R is
upper-triangular. While this decomposition is exact, it is inefficient to use because it fills in
zero elements of A with nonzero values. However, we can generate an incomplete Cholesky
T
factorization A ≈ M, M = M M 2 if we execute the Cholesky algorithm, but throw out
2
some fraction of the new nonzero values introduced by fill-in to control the number of
nonzero values to a manageable amount. Alternatively, in a modified incomplete Cholesky
factorization,wesubsumethediscardednonzerovaluesintothediagonalelements.Asimilar
procedure for Gaussian elimination yields an incomplete LU factorization, A ≈ M, M =
LU, that can be applied if A is not positive-definite. Using either of these factorizations,
A ≈ M, M = M 1 M 2 we supply the preconditioner with the syntax
x = pcg(A,b,tol,maxit,M); x = pcg(A,b,tol,maxit,M1,M2);
x = gmres (A,b,restart,tol,maxit,M1,M2);
tol sets the desired level of accuracy, maxit is the allowable number of iterations, and
restart asks GMRES to rebuild the Krylov subspace every restart iterations. The pre-
conditioner is supplied as M or as M1, M2. To specify the initial guess of the solution,