<<

. 3
( 7)



>>

3.6.1. Bi-CG. The earliest such method Bi-CG (Biconjugate gradient)
[122], [76], does not enforce a minimization principle; instead, the kth residual
must satisfy the bi-orthogonality condition
T
(3.16) rk w = 0 for all w ∈ Kk ,

where
Kk = span(ˆ0 , AT r0 , . . . , (AT )k’1 r0 )
r ˆ ˆ
is the Krylov space for AT and the vector r0 . r0 is a user-supplied vector and is
ˆˆ
often set to r0 . The algorithm gets its name because it produces sequences of
residuals {rk }, {ˆk } and search directions {pk }, {ˆk } such that bi-orthogonality
r p
T r = 0 if k = l and the search directions {p } and {ˆ } satisfy
holds, i. e. rk l
ˆ pk
k
the bi-conjugacy property

pT Apl = 0 if k = l.
ˆk

In the symmetric positive de¬nite case (3.16) is equivalent to the minimization
principle (2.2) for CG [89].
Using the notation of Chapter 2 and [191] we give an implementation of Bi-
CG making the choice r0 = r0 . This algorithmic description is explained and
ˆ
motivated in more detail in [122], [76], [78], and [191]. We do not recommend
use of Bi-CG and present this algorithm only as a basis for discussion of some
of the properties of this class of algorithms.
Algorithm 3.6.1. bicg(x, b, A, , kmax)
1. r = b ’ Ax, r = r, ρ0 = 1, p = p = 0, k = 0
ˆ ˆ

2. Do While r > b and k < kmax
2 2

(a) k = k + 1
(b) ρk = rT r, β = ρk /ρk’1
ˆ
(c) p = r + βp, p = r + β p
ˆˆ ˆ
(d) v = Ap
(e) ± = ρk /(ˆT v)
p
(f) x = x + ±p
(g) r = r ’ ±v; r = r ’ ±AT p
ˆˆ ˆ

Note that if A = AT is spd, Bi-CG produces the same iterations as CG
(but computes everything except x twice). If, however, A is not spd, there is
no guarantee that ρk in step 2b or pT Ap, the denominator in step 2e, will not
ˆ
T Ap = 0 we say that a breakdown has taken
vanish. If either ρk’1 = 0 or p ˆ
place. Even if these quantities are nonzero but very small, the algorithm can
become unstable or produce inaccurate results. While there are approaches



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




48 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

[80], [78], [148], [79], [81], to avoid some breakdowns in many variations of
Bi-CG, there are no methods that both limit storage and completely eliminate
the possibility of breakdown [74]. All of the methods we present in this section
can break down and should be implemented with that possibility in mind.
Once breakdown has occurred, one can restart the Bi-CG iteration or pass
the computation to another algorithm, such as GMRES. The possibility of
breakdown is small and certainly should not eliminate the algorithms discussed
below from consideration if there is not enough storage for GMRES to perform
well.
Breakdowns aside, there are other problems with Bi-CG. A transpose-
vector product is needed, which at the least will require additional program-
ming and may not be possible at all. The performance of the algorithm can
be erratic or even unstable with residuals increasing by several orders of mag-
nitude from one iteration to the next. Moreover, the e¬ort in computing r ˆ
at each iteration is wasted in that r makes no contribution to x. However,
ˆ
Bi-CG sometimes performs extremely well and the remaining algorithms in
this section represent attempts to capture this good performance and damp
the erratic behavior when it occurs.
We can compare the best-case performance of Bi-CG with that of GMRES.
Note that there is pk ∈ Pk such that both
¯

rk = pk (A)r0 and rk = pk (AT )ˆ0 .
(3.17) ¯ ˆ ¯ r

Hence, by the minimization property for GMRES
Bi’CG
rk RES
GM
¤ rk 2,
2

reminding us that GMRES, if su¬cient storage is available, will always
reduce the residual more rapidly than Bi-CG (in terms of iterations, but not
necessarily in terms of computational work). One should also keep in mind that
a single Bi-CG iteration requires two matrix-vector products and a GMRES
iterate only one, but that the cost of the GMRES iteration increases (in terms
of ¬‚oating-point operations) as the iteration progresses.

3.6.2. CGS. A remedy for one of the problems with Bi-CG is the Conjugate
Gradient Squared (CGS) algorithm [180]. The algorithm takes advantage of
the fact that (3.17) implies that the scalar product rT r in Bi-CG (step 2b of
ˆ
Algorithm bicg) can be represented without using AT as

rk rk = (¯k (A)r0 )T (¯k (AT )ˆ0 ) = (¯k (A)2 r0 )T r0 .
T
ˆ p p r p ˆ

The other references to AT can be eliminated in a similar fashion and an
iteration satisfying
rk = pk (A)2 r0
(3.18) ¯
is produced, where pk is the same polynomial produced by Bi-CG and used in
¯
(3.17). This explains the name, Conjugate Gradient Squared.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




49
GMRES ITERATION

The work used in Bi-CG to compute r is now used to update x. CGS
ˆ
replaces the transpose-vector product with an additional matrix-vector product
and applies the square of the Bi-CG polynomial to r0 to produce rk . This may,
of course, change the convergence properties for the worse and either improve
good convergence or magnify erratic behavior [134], [191]. CGS has the same
potential for breakdown as Bi-CG. Since p2 ∈ P2k we have
¯k

r2k RES
GM CGS
(3.19) ¤ rk 2.
2

We present an algorithmic description of CGS that we will refer to in our
discussion of TFQMR in § 3.6.4. In our description of CGS we will index the
vectors and scalars that would be overwritten in an implementation so that we
can describe some of the ideas in TFQMR later on. This description is taken
from [77].
Algorithm 3.6.2. cgs(x, b, A, , kmax)
1. x0 = x; p0 = u0 = r0 = b ’ Ax

2. v0 = Ap0 ; r0 = r0
ˆ

ˆT
3. ρ0 = r0 r0 ; k = 0

4. Do While r > b and k < kmax
2 2

(a) k = k + 1
ˆT
(b) σk’1 = r0 vk’1
(c) ±k’1 = ρk’1 /σk’1
(d) qk = uk’1 ’ ±k’1 vk’1
(e) xk = xk’1 + ±k’1 (uk’1 + qk )
rk = rk’1 ’ ±k’1 A(uk’1 + qk )
ˆT
(f) ρk = r0 rk ; βk = ρk /ρk’1
uk = rk + βk qk
pk = uk + βk (qk + βk pk’1 )
vk = Apk

Breakdowns take place when either ρk’1 or σk’1 vanish. If the algorithm
does not break down, then ±k’1 = 0 for all k. One can show [180], [77], that
if pk is the residual polynomial from (3.17) and (3.18) then
¯

(3.20) pk (z) = pk’1 (z) ’ ±k’1 z qk’1 (z)
¯ ¯ ¯

where the auxiliary polynomials qk ∈ Pk are given by q0 = 1 and for k ≥ 1 by
¯ ¯

(3.21) qk (z) = pk (z) + βk qk’1 (z).
¯ ¯ ¯

We provide no MATLAB implementation of Bi-CG or CGS. Bi-CGSTAB
is the most e¬ective representative of this class of algorithms.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




50 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

3.6.3. Bi-CGSTAB. The Bi-conjugate gradient stabilized method (Bi-
CGSTAB) [191] attempts to smooth the convergence of CGS by replacing
(3.18) with
(3.22) rk = qk (A)pk (A)r0
where
k
qk (z) = (1 ’ ωi z).
i=1

The constant ωi is selected to minimize ri = qi (A)pi (A)r0 as a function of ωi .
The examples in [191] indicate the e¬ectiveness of this approach, which can
be thought of [12] as a blend of Bi-CG and GMRES(1). The performance in
cases where the spectrum has a signi¬cant imaginary part can be improved by
constructing qk to have complex conjugate pairs of roots, [98].
There is no convergence theory for Bi-CGSTAB and no estimate of the
residuals that is better than that for CGS, (3.19). In our description of the
implementation, which is described and motivated fully in [191], we use r0 = r0
ˆ
and follow the notation of [191].
Algorithm 3.6.3. bicgstab(x, b, A, , kmax)
ˆT
1. r = b ’ Ax, r0 = r = r, ρ0 = ± = ω = 1, v = p = 0, k = 0, ρ1 = r0 r
ˆ ˆ

2. Do While r > b and k < kmax
2 2

(a) k = k + 1
(b) β = (ρk /ρk’1 )(±/ω)
(c) p = r + β(p ’ ωv)
(d) v = Ap
rT
(e) ± = ρk /(ˆ0 v)
(f) s = r ’ ±v, t = As
(g) ω = tT s/ t 2 , ρk+1 = ’ωˆ0 t
rT
2
(h) x = x + ±p + ωs
(i) r = s ’ ωt

Note that the iteration can break down in steps 2b and 2e. We provide an
implementation of Algorithm bicgstab in the collection of MATLAB codes.
The cost in storage and in ¬‚oating-point operations per iteration re-
mains bounded for the entire iteration. One must store seven vectors
(x, b, r, p, v, r0 , t), letting s overwrite r when needed. A single iteration re-
ˆ
quires four scalar products. In a situation where many GMRES iterations are
needed and matrix-vector product is fast, Bi-CGSTAB can have a much lower
average cost per iterate than GMRES. The reason for this is that the cost of
orthogonalization in GMRES can be much more than that of a matrix-vector
product if the dimension of the Krylov space is large. We will present such a
case in § 3.8.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




51
GMRES ITERATION

3.6.4. TFQMR. Now we consider the quasi-minimal residual (QMR)
family of algorithms [80], [81], [77], [37], [202]. We will focus on the transpose-
free algorithm TFQMR proposed in [77], to illustrate the quasi-minimization
idea. All algorithms in this family minimize the norm of an easily computable
quantity q = f ’ Hz, a quasi-residual over z ∈ RN . The quasi-residual is
related to the true residual by a full-rank linear transformation r = Lq and
reduction of the residual can be approximately measured by reduction in q.
The speci¬c formulation of q and L depends on the algorithm.
Returning to Algorithm cgs we de¬ne sequences
±
 uk’1 if m = 2k ’ 1 is odd,

ym =

q if m = 2k is even,
k

and ±
 (¯k (A))2 r0
p if m = 2k ’ 1 is odd,
wm =

 p (A)¯
¯k pk’1 (A)r0 if m = 2k is even.
Here the sequences {¯k } and {¯k } are given by (3.20) and (3.21). Hence, the
p q
CGS residual satis¬es
CGS
rk = w2k+1 .
We assume that CGS does not break down and, therefore, ±k = 0 for all
k ≥ 0. If we let r be the nearest integer less than or equal to a real r we
have
(3.23) Aym = (wm ’ wm+1 )/± (m’1)/2 ,
where the denominator is nonzero by assumption. We express (3.23) in matrix
form as
(3.24) AYm = Wm+1 Bm ,
where Ym is the N — m matrix with columns {yj }m , Wm+1 the N — (m + 1)
j=1
m+1
matrix with columns {wj }j=1 , and Bm is the (m + 1) — m matrix given by
« 
1 0 ... 0
¬ .·
.. .·
¬ ’1 .
1 .
¬ ·
¬ ·
. .
Bm = ¬ 0 · diag(±0 , ±0 , . . . , ± (m’1)/2 )’1 .
. .
. .
0
¬ ·
¬. ·
..
¬. ·
. ’1 1
. 
0 . . . 0 ’1

Our assumptions that no breakdown takes place imply that Km is the span
of {yj }m and hence
j=1
(3.25) xm = x0 + Ym z
for some z ∈ Rm . Therefore,

(3.26) rm = r0 ’ AYm z = Wm+1 (e1 ’ Bm z),



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




52 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

where, as in our discussion of GMRES, e1 = (1, 0, . . . , 0)T ∈ Rm . In a sense,
the conditioning of the matrix Wm+1 can be improved by scaling. Let

(3.27) „¦ = diag(ω1 , ω2 , . . . , ωm )

and rewrite (3.26) as

rm = r0 ’ AYm z = Wm+1 „¦’1 (fm+1 ’ Hm z)
(3.28) m+1

where fm+1 = „¦m+1 e1 = ωm+1 e1 and Hm = „¦m+1 Bm is bidiagonal and hence
upper Hessenberg. The diagonal entries in the matrix „¦m+1 are called weights.
The weights make it possible to derive straightforward estimates of the residual
in terms of easily computed terms in the iteration.
If Wm+1 „¦’1 were an orthogonal matrix, we could minimize rm by solving
m+1
the least squares problem

(3.29) minimizez∈Rm fm+1 ’ Hm z 2 .

The quasi-minimization idea is to solve (3.29) (thereby quasi-minimizing the
residual) to obtain zm and then set

(3.30) xm = x0 + Ym zm .

Note that if k corresponds to the iteration counter for CGS, each k produces
two approximate solutions (corresponding to m = 2k ’ 1 and m = 2k).
The solution to the least squares problem (3.29) can be done by using
Givens rotations to update the factorization of the upper Hessenberg matrix
Hm [77] and this is re¬‚ected in the implementation description below. As
one can see from that description, approximate residuals are not computed in
Algorithm tfqmr. Instead we have the quasi-residual norm

„m = fm+1 ’ Hm zm 2 .

then the columns of Wm+1 „¦’1 have
Now if we pick the weights ωi = wi 2 m+1
norm one. Hence

Wm+1 „¦’1 (fm+1
(3.31) rm = ’ Hm z) ¤ „m m + 1.
2 2
m+1

We can, therefore, base our termination criterion on „m , and we do this in
Algorithm tfqmr.
One can [80], [78], also use the quasi-minimization condition to estimate
the TFQMR residual in terms of residual polynomials.
Theorem 3.6.1. Let A be an N — N matrix and let x0 , b ∈ RN be given.
Assume that the TFQMR iteration does not break down or terminate with
the exact solution for k iterations and let 1 ¤ m ¤ 2k. Let rm RES be the
GM

residual for the mth GMRES iteration. Let ξm be the smallest singular value
of Wm+1 „¦’1 . Then if ξm > 0
m+1

„m ¤ rm RES
GM
(3.32) 2 /ξm .




Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




53
GMRES ITERATION

Proof. We may write

rm RES = r0 + Ym zm RES
GM GM


and obtain

= Wm+1 „¦’1 (fm+1 ’Hm zm RES )
rm RES
GM GM
≥ ξm fm+1 ’Hm zm RES
GM
2.
2 2
m+1

Hence, by the quasi-minimization property

rm RES
GM
≥ ξm „m
2

as desired.
As a corollary we have a ¬nite termination result and a convergence
estimate for diagonalizable matrices similar to Theorem 3.1.3.
Corollary 3.6.1. Let A be an N — N matrix and let x0 , b ∈ RN be given.
Then within (N + 1)/2 iterations the TFQMR iteration will either break down
or terminate with the solution.
We apply Theorem 3.1.3 to (3.32) and use (3.31) to obtain the following
result.
Theorem 3.6.2. Let A = V ΛV ’1 be a nonsingular diagonalizable matrix.
Assume that TFQMR does not break down or terminate with the solution for k
iterations. For 1 ¤ m ¤ 2k let ξm be the smallest singular value of Wm+1 „¦’1
m+1
and let xm be given by (3.30). Then, if ξm > 0,

rm 2 ’1
¤ m + 1ξm κ(V ) max |φ(z)|
r0 z∈σ(A)
2

for all φ ∈ Pm .
The implementation below follows [77] with r0 = r0 used throughout.
ˆ
Algorithm 3.6.4. tfqmr(x, b, A, , kmax)
1. k = 0, w1 = y1 = r0 = b ’ Ax, u1 = v = Ay1 , d = 0
T
ρ0 = r0 r0 , „ = r 2 , θ = 0, · = 0

2. Do While k < kmax
(a) k = k + 1
T
(b) σk’1 = r0 v, ± = ρk’1 /σk’1 , y2 = y1 ’ ±k’1 v, u2 = Ay2
(c) For j = 1, 2 (m = 2k ’ 2 + j)
i. w = w ’ ±k’1 uj
d = yj + (θ2 ·/±k’1 )d
ii.

θ = w 2 /„ , c = 1/ 1 + θ2
iii.
„ = „ θc, · = c2 ±k’1
iv.
v. x = x + ·d

vi. If „ m + 1 ¤ b 2 terminate successfully
T
(d) ρk = r0 w, β = ρk /ρk’1



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




54 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

(e) y1 = w + βy2 , u1 = Ay1
(f) v = u1 + β(u2 + βv)

Note that y2 and u2 = Ay2 need only be computed if the loop in step 2c
does not terminate when j = 1. We take advantage of this in the MATLAB
code tfqmr to potentially reduce the matrix-vector product cost by one.

3.7. Examples for GMRES iteration
These examples use the code gmres from the collection of MATLAB codes. The
inputs, described in the comment lines, are the initial iterate x0 , the right-hand
side vector b, a MATLAB function for the matrix-vector product, and iteration
parameters to specify the maximum number of iterations and the termination
criterion. We do not pass the preconditioner to the GMRES code, rather
we pass the preconditioned problem. In all the GMRES examples, we limit
the number of GMRES iterations to 60. For methods like GMRES(m), Bi-
CGSTAB, CGNR, and TFQMR, whose storage requirements do not increase
with the number of iterations, we allow for more iterations.
We consider the discretization of the partial di¬erential equation

(Lu)(x, y) = ’(uxx (x, y) + uyy (x, y)) + a1 (x, y)ux (x, y)
(3.33)
+a2 (x, y)uy (x, y) + a3 (x, y)u(x, y) = f (x, y)

on 0 < x, y < 1 subject to homogeneous Dirichlet boundary conditions

u(x, 0) = u(x, 1) = u(0, y) = u(1, y) = 0, 0 < x, y < 1.

For general coe¬cients {aj } the operator L is not self-adjoint and its discretiza-
tion is not symmetric. As in § 2.7 we discretize with a ¬ve-point centered dif-
ference scheme with n2 points and mesh width h = 1/(n + 1). The unknowns
are
uij ≈ u(xi , xj )

where xi = ih for 1 ¤ i ¤ n. We compute Lu in a matrix-free manner
as we did in § 2.7. We let n = 31 to create a system with 961 unknowns.
As before, we expect second-order accuracy and terminate the iteration when
rk 2 ¤ h2 b 2 ≈ 9.8 — 10’4 b 2 .
As a preconditioner we use the fast Poisson solver fish2d described in § 2.7.
The motivation for this choice is similar to that for the conjugate gradient
computations and has been partially explained theoretically in [33], [34], [121],
and [139]. If we let Gu denote the action of the Poisson solver on u, the
preconditioned system is GLu = Gf .
For the computations reported in this section we took

a1 (x, y) = 1, a2 (x, y) = 20y, and a3 (x, y) = 1.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




55
GMRES ITERATION

u0 = 0 was the initial iterate. As in § 2.7 we took the right hand side so that
the exact solution was the discretization of

10xy(1 ’ x)(1 ’ y) exp(x4.5 ).

In Figs. 3.1 and 3.2 we plot iteration histories corresponding to precondi-
tioned or unpreconditioned GMRES. For the restarted iteration, we restarted
the algorithm after three iterations, GMRES(3) in the terminology of § 3.4.
Note that the plots are semilog plots of r 2 / b 2 and that the right-hand
sides in the preconditioned and unpreconditioned cases are not the same. In
both ¬gures the solid line is the history of the preconditioned iteration and the
dashed line that of the unpreconditioned. Note the importance of precondi-
tioning. The MATLAB flops command indicates that the unpreconditioned
iteration required 7.6 million ¬‚oating point operations and converged in 56
iterations. The preconditioned iteration required 1.3 million ¬‚oating point op-
erations and terminated after 8 iterations. Restarting after 3 iterations is more
costly in terms of the iteration count than not restarting. However, even in
the unpreconditioned case, the restarted iteration terminated successfully. In
this case the MATLAB flops command indicates that the unpreconditioned
iteration terminated after 223 iterations and 9.7 million ¬‚oating point opera-
tions and the preconditioned iteration required 13 iterations and 2.6 million
¬‚oating point operations.
The execution times in our environment indicate that the preconditioned
iterations are even faster than the di¬erence in operation counts indicates.
The reason for this is that the FFT routine is a MATLAB intrinsic and is not
interpreted MATLAB code. In other computing environments preconditioners
which vectorize or parallelize particularly well might perform in the same way.

3.8. Examples for CGNR, Bi-CGSTAB, and TFQMR iteration
When a transpose is available, CGNR (or CGNE) is an alternative to GMRES
that does not require storage of the iteration history. For elliptic problems like
(3.33), the transpose (adjoint) of L is given by

L— u = ’uxx ’ uyy ’ a1 ux ’ a2 uy + a3 u

as one can derive by integration by parts. To apply CGNR to the precondi-
tioned problem, we require the transpose of GL, which is given by

(GL)— u = L— G— u = L— Gu.

The cost of the application of the transpose of L or GL is the same as that
of the application of L or GL itself. Hence, in the case where the cost of the
matrix-vector product dominates the computation, a single CGNR iteration
costs roughly the same as two GMRES iterations.
However, CGNR has the e¬ect of squaring the condition number, this e¬ect
has serious consequences as Fig. 3.3 shows. The dashed line corresponds to



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




56 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS


0
10




-1
10




Relative Residual
-2
10




-3
10




-4
10
0 10 20 30 40 50 60
Iterations




Fig. 3.1. GMRES for 2-D elliptic equation.
0
10




-1
10
Relative Residual




-2
10




-3
10




-4
10
0 50 100 150 200 250
Iterations




Fig. 3.2. GMRES(3) for 2-D elliptic equation.


CGNR on the unpreconditioned problem, i.e., CG applied to L— Lu = L— f .
The matrix corresponding to the discretization of L— L has a large condition
number. The solid line corresponds to the CGNR applied to the preconditioned
problem, i.e., CG applied to L— G2 Lu = L— G2 f . We limited the number of
iterations to 310 = 10m and the unpreconditioned formulation of CGNR had
made very little progress after 310 iterations. This is a result of squaring
the large condition number. The preconditioned formulation did quite well,
converging in eight iterations. The MATLAB flops command indicates that
the unpreconditioned iteration required 13.7 million ¬‚oating-point operations



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




57
GMRES ITERATION

and the preconditioned iteration 2.4 million. However, as you are asked to
investigate in Exercise 3.9.8, the behavior of even the preconditioned iteration
can also be poor if the coe¬cients of the ¬rst derivative terms are too large.
CGNE has similar performance; see Exercise 3.9.9.

0
10




-1
10
Relative Residual




-2
10




-3
10




-4
10
0 50 100 150 200 250 300 350
Iterations




Fig. 3.3. CGNR for 2-D elliptic equation.

The examples for Bi-CGSTAB use the code bicgstab from the collection
of MATLAB codes. This code has the same input/output arguments as gmres.
We applied Bi-CGSTAB to both the preconditioned and unpreconditioned
forms of (3.33) with the same data and termination criterion as for the GMRES
example.
We plot the iteration histories in Fig. 3.4. As in the previous examples, the
solid line corresponds to the preconditioned iteration and the dashed line to
the unpreconditioned iteration. Neither convergence history is monotonic, with
the unpreconditioned iteration being very irregular, but not varying by many
orders of magnitude as a CGS or Bi-CG iteration might. The preconditioned
iteration terminated in six iterations (12 matrix-vector products) and required
roughly 1.7 million ¬‚oating-point operations. The unpreconditioned iteration
took 40 iterations and 2.2 million ¬‚oating-point operations. We see a
considerable improvement in cost over CGNR and, since our unpreconditioned
matrix-vector product is so inexpensive, a better cost/iterate than GMRES in
the unpreconditioned case.
We repeat the experiment with TFQMR as our linear solver. We use
the code tfqmr that is included in our collection. In Fig.√3.5 we plot the
convergence histories in terms of the approximate residual 2k + 1„2k given
by (3.31). We plot only full iterations (m = 2k) except, possibly, for the
¬nal iterate. The approximate residual is given in terms of the quasi-residual
norm „m rather than the true residual, and the graphs are monotonic. Our



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




58 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS


1
10



0
10




Relative Residual
-1
10



-2
10



-3
10



-4
10
0 5 10 15 20 25 30 35 40
Iterations




Fig. 3.4. Bi-CGSTAB for 2-D elliptic equation.
0
10



-1
10
Relative Approximate Residual




-2
10



-3
10



-4
10



-5
10
0 10 20 30 40 50 60 70
Iterations




Fig. 3.5. TFQMR for 2-D elliptic equation.


termination criterion is based on the quasi-residual and (3.31) as described in

§ 3.6.4, stopping the iteration when „m m + 1/ b 2 ¤ h’2 . This led to actual
relative residuals of 7 — 10’5 for the unpreconditioned iteration and 2 — 10’4
for the preconditioned iteration, indicating that it is reasonable to use the
quasi-residual estimate (3.31) to make termination decisions.
The unpreconditioned iteration required roughly 4.1 million ¬‚oating-point
operations and 67 iterations for termination. The preconditioned iteration
was much more e¬cient, taking 1.9 million ¬‚oating-point operations and 7
iterations. The plateaus in the graph of the convergence history for the



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




59
GMRES ITERATION

unpreconditioned iteration are related (see [196]) to spikes in the corresponding
graph for a CGS iteration.




Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




60 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

3.9. Exercises on GMRES
3.9.1. Give an example of a matrix that is not diagonalizable.

3.9.2. Prove Theorem 3.1.4 and Theorem 3.1.5.

3.9.3. Prove that the Arnoldi process (unless it terminates prematurely with a
solution) produces a basis for Kk+1 that satis¬es (3.8).

3.9.4. Duplicate Table 3.1 and add two more columns corresponding to classical
Gram“Schmidt orthogonalization with reorthogonalization at every step
and determine when the test in (3.10) detects loss of orthogonality. How
do your conclusions compare with those in [160]?

3.9.5. Let A = J» be an elementary Jordan block of size N corresponding to
an eigenvalue » = 0. This means that
« 
» 1 0
¬0 ·
» 1 0
¬ ·
¬ ·
.. .. .. ..
¬ ·
. . . .
J» = ¬ ·.
¬ ·
0 » 1 0·
¬
¬ ·
 1
0 »
0 »

Give an example of x0 , b ∈ RN , and » ∈ R such that the GMRES
residuals satisfy rk 2 = r0 2 for all k < N .

3.9.6. Consider the centered di¬erence discretization of

’u + u + u = 1, u(0) = u(1) = 0.

Solve this problem with GMRES (without preconditioning) and then
apply as a preconditioner the map M de¬ned by

’(M f ) = f, (M f )(0) = (M f )(1) = 0.

That is, precondition with a solver for the high-order term in the
di¬erential equation using the correct boundary conditions. Try this
for meshes with 50, 100, 200, . . . points. How does the performance of
the iteration depend on the mesh? Try other preconditioners such as
Gauss“Seidel and Jacobi. How do other methods such as CGNR, CGNE,
Bi-CGSTAB, and TFQMR perform?

3.9.7. Modify the MATLAB GMRES code to allow for restarts. One way to do
this is to call the gmres code with another code which has an outer loop
to control the restarting as in Algorithm gmresm. See [167] for another
example. How would you test for failure in a restarted algorithm? What
kinds of failures can occur? Repeat Exercise 6 with a limit of three
GMRES iterations before a restart.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




61
GMRES ITERATION

3.9.8. Duplicate the results in § 3.7 and § 3.8. Compare the execution times
in your computing environment. Experiment with di¬erent values of m.
Why might GMRES(m) be faster than GMRES? Vary the dimension and
the coe¬cient ay . Try ay = 5y. Why is the performance of the methods
sensitive to ay and ax ?

3.9.9. Apply CGNE to the PDE in § 3.7. How does the performance di¬er from
CGNR and GMRES?

3.9.10. Consider preconditioning from the right. For the PDE in § 3.7, apply
GMRES, CGNR, Bi-CGSTAB, TFQMR, and/or CGNE, to the equation
LGw = f and then recover the solution from u = Gw. Compare timings
and solution accuracy with those obtained by left preconditioning.

3.9.11. Are the solutions to (3.33) in § 3.7 and § 3.8 equally accurate? Compare
the ¬nal results with the known solution.

3.9.12. Implement Bi-CG and CGS and try them on some of the examples in
this chapter.




Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




62 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS




Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




Chapter 4

Basic Concepts and Fixed-Point Iteration




This part of the book is about numerical methods for the solution of systems
of nonlinear equations. Throughout this part, we will let · denote a norm
on RN as well as the induced matrix norm.
We begin by setting the notation. We seek to solve
(4.1) F (x) = 0.
Here F : RN ’ RN . We denote the ith component of F by fi . If the
components of F are di¬erentiable at x ∈ RN we de¬ne the Jacobian matrix
F (x) by
‚fi
F (x)ij = (x).
‚xj
The Jacobian matrix is the vector analog of the derivative. We may express
the fundamental theorem of calculus as follows.
Theorem 4.0.1. Let F be di¬erentiable in an open set „¦ ‚ RN and let
x— ∈ „¦. Then for all x ∈ „¦ su¬ciently near x—
1

F (x— + t(x ’ x— ))(x ’ x— ) dt.
F (x) ’ F (x ) =
0


4.1. Types of convergence
Iterative methods can be classi¬ed by the rate of convergence.
Definition 4.1.1. Let {xn } ‚ RN and x— ∈ RN . Then
• xn ’ x— q-quadratically if xn ’ x— and there is K > 0 such that
xn+1 ’ x— ¤ K xn ’ x— 2 .

• xn ’ x— q-superlinearly with q-order ± > 1 if xn ’ x— and there is
K > 0 such that
xn+1 ’ x— ¤ K xn ’ x— ± .

• xn ’ x— q-superlinearly if
xn+1 ’ x—
lim = 0.
xn ’ x—
n’∞

65

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




66 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

• xn ’ x— q-linearly with q-factor σ ∈ (0, 1) if
xn+1 ’ x— ¤ σ xn ’ x—
for n su¬ciently large.
Definition 4.1.2. An iterative method for computing x— is said to be
locally (q-quadratically, q-superlinearly, q-linearly, etc.) convergent if the
iterates converge to x— (q-quadratically, q-superlinearly, q-linearly, etc.) given
that the initial data for the iteration is su¬ciently good.
Note that a q-superlinearly convergent sequence is also q-linearly conver-
gent with q-factor σ for any σ > 0. A q-quadratically convergent sequence is
q-superlinearly convergent with q-order 2.
In general, a q-superlinearly convergent method is preferable to a q-linearly
convergent one if the cost of a single iterate is the same for both methods. We
will see that often a method that is more slowly convergent, in terms of the
convergence types de¬ned above, can have a cost/iterate so low that the slow
iteration is more e¬cient.
Sometimes errors are introduced into the iteration that are independent of
the errors in the approximate solution. An example of this would be a problem
in which the nonlinear function F is the output of another algorithm that
provides error control, such as adjustment of the mesh size in an approximation
of a di¬erential equation. One might only compute F to low accuracy in the
initial phases of the iteration to compute a solution and tighten the accuracy
as the iteration progresses. The r-type convergence classi¬cation enables us to
describe that idea.
Definition 4.1.3. Let {xn } ‚ RN and x— ∈ RN . Then {xn } converges
to x— r-(quadratically, superlinearly, linearly) if there is a sequence {ξn } ‚ R
converging q-(quadratically, superlinearly, linearly) to zero such that
xn ’ x— ¤ ξn .
We say that {xn } converges r-superlinearly with r-order ± > 1 if ξn ’ 0
q-superlinearly with q-order ±.

4.2. Fixed-point iteration
Many nonlinear equations are naturally formulated as ¬xed-point problems
(4.2) x = K(x)
where K, the ¬xed-point map, may be nonlinear. A solution x— of (4.2) is
called a ¬xed point of the map K. Such problems are nonlinear analogs of the
linear ¬xed-point problems considered in Chapter 1. In this section we analyze
convergence of ¬xed-point iteration, which is given by
(4.3) xn+1 = K(xn ).
This iterative method is also called nonlinear Richardson iteration, Picard
iteration, or the method of successive substitution.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




67
FIXED-POINT ITERATION

Before discussing convergence of ¬xed-point iteration we make two de¬ni-
tions.
Definition 4.2.1. Let „¦ ‚ RN and let G : „¦ ’ RM . G is Lipschitz
continuous on „¦ with Lipschitz constant γ if

G(x) ’ G(y) ¤ γ x ’ y

for all x, y ∈ „¦.
Definition 4.2.2. Let „¦ ‚ RN . K : „¦ ’ RN is a contraction mapping
on „¦ if K is Lipschitz continuous on „¦ with Lipschitz constant γ < 1.
The standard result for ¬xed-point iteration is the Contraction Mapping
Theorem [9]. Compare the proof to that of Lemma 1.2.1 and Corollary 1.2.1.
Theorem 4.2.1. Let „¦ be a closed subset of RN and let K be a contraction
mapping on „¦ with Lipschitz constant γ < 1 such that K(x) ∈ „¦ for all x ∈ „¦.
Then there is a unique ¬xed point of K, x— ∈ „¦, and the iteration de¬ned by
(4.3) converges q-linearly to x— with q-factor γ for all initial iterates x0 ∈ „¦.
Proof. Let x0 ∈ „¦. Note that {xn } ‚ „¦ because x0 ∈ „¦ and K(x) ∈ „¦
whenever x ∈ „¦. The sequence {xn } remains bounded since for all i ≥ 1

xi+1 ’ xi = K(xi ) ’ K(xi’1 ) ¤ γ xi ’ xi’1 . . . ¤ γ i x1 ’ x0 ,

and therefore
n’1
xn ’ x0 = i=0 xi+1 ’ xi

n’1 n’1 i
¤ xi+1 ’ xi ¤ x1 ’ x0 i=0 γ
i=0


¤ x1 ’ x0 /(1 ’ γ).

Now, for all n, k ≥ 0,

xn+k ’ xn = K(xn+k’1 ) ’ K(xn’1 )

¤ γ xn+k’1 ’ xn’1

¤ γ K(xn+k’2 ) ’ K(xn’2 )

¤ γ 2 xn+k’2 ’ xn’2 ¤ . . . ¤ γ n xk ’ x0

¤ γ n x1 ’ x0 /(1 ’ γ).

Hence
lim xn+k ’ xn = 0
n,k’∞

and therefore the sequence {xn } is a Cauchy sequence and has a limit x— [162].
If K has two ¬xed points x— and y — in „¦, then

x— ’ y — = K(x— ) ’ K(y — ) ¤ γ x— ’ y — .



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




68 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Since γ < 1 this can only hold if x— ’ y — = 0, i. e. x— = y — . Hence the ¬xed
point is unique.
Finally we note that

xn+1 ’ x— = K(xn ) ’ K(x— ) ¤ γ xn ’ x— ,

which completes the proof.
One simple application of this theorem is to numerical integration of
ordinary di¬erential equations by implicit methods. For example, consider
the initial value problem

y = f (y), y(t0 ) = y0

and its solution by the backward Euler method with time step h,

y n+1 ’ y n = hf (y n+1 ).
(4.4)

In (4.4) y k denotes the approximation of y(t0 + kh). Advancing in time from
y n to y n+1 requires solution of the ¬xed-point problem

y = K(y) = y n + hf (y).

For simplicity we assume that f is bounded and Lipschitz continuous, say

|f (y)| ¤ M and |f (x) ’ f (y)| ¤ M |x ’ y|

for all x, y.
If hM < 1, we may apply the contraction mapping theorem with „¦ = R
since for all x, y

K(x) ’ K(y) = h|f (x) ’ f (y)| ¤ hM |x ’ y|.

From this we conclude that for h su¬ciently small we can integrate forward
in time and that the time step h can be chosen independently of y. This time
step may well be too small to be practical, and the methods based on Newton™s
method which we present in the subsequent chapters are used much more often
[83]. This type of analysis was used by Picard to prove existence of solutions
of initial value problems [152]. Nonlinear variations of the classical stationary
iterative methods such as Gauss“Seidel have also been used; see [145] for a
more complete discussion of these algorithms.

4.3. The standard assumptions
We will make the standard assumptions on F .
Assumption 4.3.1.
1. Equation 4.1 has a solution x— .

2. F : „¦ ’ RN —N is Lipschitz continuous with Lipschitz constant γ.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




69
FIXED-POINT ITERATION

3. F (x— ) is nonsingular.
These assumptions can be weakened [108] without sacri¬cing convergence
of the methods we consider here. However the classical result on quadratic
convergence of Newton™s method requires them and we make them throughout.
Exercise 5.7.1 illustrates one way to weaken the standard assumptions.
Throughout this part, we will always denote a root of F by x— . We let
B(r) denote the ball of radius r about x—

B(r) = {x | e < r},

where
e = x ’ x— .
The notation introduced above will be used consistently. If xn is the nth iterate
of a sequence, en = xn ’ x— is the error in that iterate.
Lemma 4.3.1 is an important consequence of the standard assumptions.
Lemma 4.3.1. Assume that the standard assumptions hold. Then there is
δ > 0 so that for all x ∈ B(δ)

F (x) ¤ 2 F (x— ) ,
(4.5)

F (x)’1 ¤ 2 F (x— )’1 ,
(4.6)
and
F (x— )’1 ’1
e /2 ¤ F (x) ¤ 2 F (x— ) e .
(4.7)
Proof. Let δ be small enough so that B(δ) ‚ „¦. For all x ∈ B(δ) we have

F (x) ¤ F (x— ) + γ e .

Hence (4.5) holds if γδ < F (x— ) .
The next result (4.6) is a direct consequence of the Banach Lemma if

I ’ F (x— )’1 F (x) < 1/2.

This will follow from
F (x— )’1 ’1
δ<

since then
I ’ F (x— )’1 F (x) = F (x— )’1 (F (x— ) ’ F (x))
(4.8)
¤ γ F (x— )’1 e ¤ γδ F (x— )’1 < 1/2.

To prove the ¬nal inequality (4.7), we note that if x ∈ B(δ) then x— + te ∈
B(δ) for all 0 ¤ t ¤ 1. We use (4.5) and Theorem 4.0.1 to obtain
1
F (x— + te) e dt ¤ 2 F (x— ) e
F (x) ¤
0




Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




70 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

which is the right inequality in (4.7).
To prove the left inequality in (4.7) note that
1
(x— )’1 F (x) (x— )’1 F (x— + te)e dt
F =F
0

1
(I ’ F (x— )’1 F (x— + te))e dt,
=e’
0

and hence, by (4.8)
1
— ’1
I ’ F (x— )’1 F (x— + te)dt ) ≥ e /2.
F (x ) F (x) ≥ e (1 ’
0

Therefore
e /2 ¤ F (x— )’1 F (x) ¤ F (x— )’1 F (x) ,
which completes the proof.




Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




Chapter 5

Newton™s Method




5.1. Local convergence of Newton™s method
We will describe iterative methods for nonlinear equations in terms of the
transition from a current iterate xc to a new iterate x+ . In this language,
Newton™s method is
x+ = xc ’ F (xc )’1 F (xc ).
(5.1)
We may also view x+ as the root of the two-term Taylor expansion or linear
model of F about xc
Mc (x) = F (xc ) + F (xc )(x ’ xc ).
In the context of single systems this method appeared in the 17th century
[140], [156], [13], [137].
The convergence result on Newton™s method follows from Lemma 4.3.1.
Theorem 5.1.1. Let the standard assumptions hold. Then there are K > 0
and δ > 0 such that if xc ∈ B(δ) the Newton iterate from xc given by (5.1)
satis¬es
e+ ¤ K ec 2 .
(5.2)
Proof. Let δ be small enough so that the conclusions of Lemma 4.3.1 hold.
By Theorem 4.0.1
1
’1 ’1
(F (xc ) ’ F (x— + tec ))ec dt.
e+ = ec ’ F (xc ) F (xc ) = F (xc )
0

By Lemma 4.3.1 and the Lipschitz continuity of F
e+ ¤ (2 F (x— )’1 )γ ec 2 /2.
This completes the proof of (5.2) with K = γ F (x— )’1 .
The proof of convergence of the complete Newton iteration will be complete
if we reduce δ if needed so that Kδ < 1.
Theorem 5.1.2. Let the standard assumptions hold. Then there is δ such
that if x0 ∈ B(δ) the Newton iteration
xn+1 = xn ’ F (xn )’1 F (xn )
converges q-quadratically to x— .
71

Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




72 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Proof. Let δ be small enough so that the conclusions of Theorem 5.1.1
hold. Reduce δ if needed so that Kδ = · < 1. Then if n ≥ 0 and xn ∈ B(δ),
Theorem 5.1.1 implies that
2
(5.3) en+1 ¤ K en ¤ · en < en

and hence xn+1 ∈ B(·δ) ‚ B(δ). Therefore if xn ∈ B(δ) we may continue the
iteration. Since x0 ∈ B(δ) by assumption, the entire sequence {xn } ‚ B(δ).
(5.3) then implies that xn ’ x— q-quadratically.
The assumption that the initial iterate be “su¬ciently near” the solution
(x0 ∈ B(δ)) may seem arti¬cial at ¬rst look. There are, however, many
situations in which the initial iterate is very near the root. Two examples are
implicit integration of ordinary di¬erential equations and di¬erential algebraic
equations, [105], [83], [16], where the initial iterate is derived from the solution
at the previous time step, and solution of discretizations of partial di¬erential
equations, where the initial iterate is an interpolation of a solution from a
coarser computational mesh [99], [126]. Moreover, when the initial iterate is
far from the root and methods such as those discussed in Chapter 8 are needed,
local convergence results like those in this and the following two chapters
describe the terminal phase of the iteration.

5.2. Termination of the iteration
We can base our termination decision on Lemma 4.3.1. If x ∈ B(δ), where δ
is small enough so that the conclusions of Lemma 4.3.1 hold, then if F (x— ) is
well conditioned, we may terminate the iteration when the relative nonlinear
residual F (x) / F (x0 ) is small. We have, as a direct consequence of applying
Lemma 4.3.1 twice, a nonlinear form of Lemma 1.1.1.
Lemma 5.2.1. Assume that the standard assumptions hold. Let δ > 0 be
small enough so that the conclusions of Lemma 4.3.1 hold for x ∈ B(δ). Then
for all x ∈ B(δ)

4κ(F (x— )) e
e F (x)
¤ ¤ ,
4 e0 κ(F (x— )) F (x0 ) e0

where κ(F (x— )) = F (x— ) F (x— )’1 is the condition number of F (x— )
relative to the norm · .
From Lemma 5.2.1 we conclude that if F (x— ) is well conditioned, the size
of the relative nonlinear residual is a good indicator of the size of the error.
However, if there is error in the evaluation of F or the initial iterate is near a
root, a termination decision based on the relative residual may be made too
late in the iteration or, as we will see in the next section, the iteration may
not terminate at all. We raised this issue in the context of linear equations
in Chapter 1 when we compared (1.2) and (1.4). In the nonlinear case, there
is no “right-hand side” to use as a scaling factor and we must balance the
relative and absolute size of the nonlinear residuals in some other way. In all



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




73
NEWTON™S METHOD

of our MATLAB codes and algorithms for nonlinear equations our termination
criterion is to stop the iteration if

(5.4) F (x) ¤ „r F (x0 ) + „a ,

where the relative error tolerance „r and absolute error tolerance „a are input
to the algorithm. Combinations of relative and absolute error tolerances are
commonly used in numerical methods for ordinary di¬erential equations or
di¬erential algebraic equations [16], [21], [23], [151].
Another way to decide whether to terminate is to look at the Newton step

s = ’F (xc )’1 F (xc ) = x+ ’ xc ,

and terminate the iteration when s is su¬ciently small. This criterion is
based on Theorem 5.1.1, which implies that

ec = s + O( ec 2 ).

Hence, near the solution s and ec are essentially the same size.
For methods other than Newton™s method, the relation between the step
and the error is less clear. If the iteration is q-linearly convergent, say, then

e+ ¤ σ ec

implies that
(1 ’ σ) ec ¤ s ¤ (1 + σ) ec .
Hence the step is a reliable indicator of the error provided σ is not too near
1. The di¬erential equations codes discussed in [16], [21], [23], and [151], make
use of this in the special case of the chord method.
However, in order to estimate ec this way, one must do most of the work
needed to compute x+ , whereas by terminating on small relative residuals or
on a condition like (5.4) the decision to stop the iteration can be made before
computation and factorization of F (xc ). If computation and factorization of
the Jacobian are inexpensive, then termination on small steps becomes more
attractive.

5.3. Implementation of Newton™s method
In this section we assume that direct methods will be used to solve the linear
equation for the Newton step. Our examples and discussions of operation
counts make the implicit assumption that the Jacobian is dense. However,
the issues for sparse Jacobians when direct methods are used to solve linear
systems are very similar. In Chapter 6 we consider algorithms in which iterative
methods are used to solve the equation for the Newton step.
In order to compute the Newton iterate x+ from a current point xc one
must ¬rst evaluate F (xc ) and decide whether to terminate the iteration. If
one decides to continue, the Jacobian F (xc ) must be computed and factored.



Buy this book from SIAM at http://www.ec-securehost.com/SIAM/FR16.html.
Copyright ©1995 by the Society for Industrial and Applied Mathematics. This electronic version is for personal use and may not be duplicated or distributed.




74 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Then the step is computed as the solution of F (xc )s = ’F (xc ) and the iterate
is updated x+ = xc + s. Of these steps, the evaluation and factorization of F
are the most costly. Factorization of F in the dense case costs O(N 3 ) ¬‚oating-

<<

. 3
( 7)



>>