[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

δ<

2γ

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-