Page 298 - Numerical Methods for Chemical Engineering

Discretized PDEs with more than two spatial dimensions 287

It can also be used without even having to store A in memory if we supply a rou-
tine that returns the vector Av for an input vector v. Such use of pcg is demonstrated in
BVP 2D Poisson FD cg.m.

The generalized minimum residual (GMRES) Krylov subspace method

Theefficiencyoftheconjugategradientmethodleadsustoconsidertheexistenceofmethods
that solve linear systems in a similar manner but without requiring positive-definiteness.
Such methods exist and are invoked in MATLAB by the keywords bicg, bicgstab, and gmres.
We present here a brief discussion of gmres.
We wish to solve Ax = b, but, as A is no longer required to be positive-definite, we
introduce the term residual for the vector that we have previously called the steepest descent
[0]
direction, r = Ax − b. Given an initial guess x , the initial residual is
r [0] = b − Ax [0] (6.145)
[1]
We wish to generate a sequence x , x [2] , . . . such that the sequence of residual norms
[1]
[2]
r , r ,...converges to zero. The GMRES method looks for new solution estimates
[0]
that minimize the residual among the members of a Krylov subspace of r , that at order
k + 1is
[0] [0] [0] 2 [0] k [0]
K k+1 r ; A ≡ span r , Ar , A r ,..., A r (6.146)
This subspace is sufficiently large to include all searches based on taking some step a [k] in
the residual direction at each iteration,
[k] [k]
x [k+1] = x [k] + a r (6.147)
This follows from the relation among residuals at successive iterations,

b − Ax [k+1] = b − Ax [k] − a [k] Ar [k]
(6.148)
[k]
[k]
[k+1]

r = I − a A r
from which we obtain
k−1

r [k] = I − a [ j] A r [ j] (6.149)
j=0
[0]
We thus find that r [k] ∈ K k+1 (r ; A). It also follows that
[0]

x [k] = x [0] + p [k] p [k] ∈ K k r ; A (6.150)
In GMRES, we compute successive estimates
x [1] = x [0] + p [1] p [1] = α 10 r [0]
x [2] = x [0] + p [2] p [2] = α 20 r [0] + α 21 Ar [0] (6.151)
2 [0]
x [3] = x [0] + p [3] p [3] = α 30 r [0] + α 31 Ar [0] + α 32 A r
[0]
by finding the p [k] ∈ K k (r ; A) that minimizes the 2-norm of the residual,
& & 2 & [0] & 2 & [0] & 2 [0]
r = b − A(x + p ) ≤ b − A(x + p) ∀p ∈ K k r ; A (6.152)
[k] &
& [k] & & & &
2 2 2