is used to solve exit distribution problems in radiative transfer.

We will discretize the equation with the composite midpoint rule. Here we

approximate integrals on [0, 1] by

N

1 1

f (µ) dµ ≈ f (µj ),

N

0 j=1

where µi = (i ’ 1/2)/N for 1 ¤ i ¤ N . The resulting discrete problem is

« ’1

N

c µi xj

F (x)i = xi ’ 1 ’

(5.22) .

2N µ + µj

j=1 i

It is known [132] that both (5.21) and the discrete analog (5.22) have two

solutions for c ∈ (0, 1). Only one of these solutions has physical meaning,

however. Newton™s method, with the 0 function or the constant function with

value 1 as initial iterate, will ¬nd the physically meaningful solution [110].

The standard assumptions hold near either solution [110] for c ∈ [0, 1). The

discrete problem has its own physical meaning [41] and we will attempt to

solve it to high accuracy. Because H has a singularity at µ = 0, the solution of

the discrete problem is not even a ¬rst-order accurate approximation to that

of the continuous problem.

In the computations reported here we set N = 100 and c = .9. We used the

function identically one as initial iterate. We computed all Jacobians with the

MATLAB code diffjac. We set the termination parameters „r = „a = 10’6 .

We begin with a comparison of Newton™s method and the chord method.

We present some iteration statistics in both tabular and graphical form. In

Table 5.1 we tabulate the iteration counter n, F (xn ) ∞ / F (x0 ) ∞ , and the

ratio

Rn = F (xn ) ∞ / F (xn’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.

88 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

Table 5.1

Comparison of Newton and chord iteration.

n F (xn ) ∞ / F (x0 ) Rn F (xn ) ∞ / F (x0 ) Rn

∞ ∞

0 1.000e+00 1.000e+00

1 1.480e-01 1.480e-01 1.480e-01 1.480e-01

2 2.698e-03 1.823e-02 3.074e-02 2.077e-01

3 7.729e-07 2.865e-04 6.511e-03 2.118e-01

4 1.388e-03 2.132e-01

5 2.965e-04 2.136e-01

6 6.334e-05 2.136e-01

7 1.353e-05 2.136e-01

8 2.891e-06 2.136e-01

for n ≥ 1 for both Newton™s method and the chord method. In Fig. 5.1

the solid curve is a plot of F (xn ) ∞ / F (x0 ) ∞ for Newton™s method and

the dashed curve a plot of F (xn ) ∞ / F (x0 ) ∞ for the chord method. The

concave iteration track is the signature of superlinear convergence, the linear

track indicates linear convergence. See Exercise 5.7.15 for a more careful

examination of this concavity. The linear convergence of the chord method

can also be seen in the convergence of the ratios Rn to a constant value.

For this example the MATLAB flops command returns an fairly accurate

indicator of the cost. The Newton iteration required over 8.8 million ¬‚oating-

point operations while the chord iteration required under 3.3 million. In

Exercise 5.7.18 you are asked to obtain timings in your environment and will

see how the ¬‚oating-point operation counts are related to the actual timings.

The good performance of the chord method is no accident. If we think

of the chord method as a preconditioned nonlinear Richardson iteration with

the initial Jacobian playing the role of preconditioner, then the chord method

should be much more e¬cient than Newton™s method if the equations for the

steps can be solved cheaply. In the case considered here, the cost of a single

evaluation and factorization of the Jacobian is far more than the cost of the

entire remainder of the chord iteration.

The ideas in Chapters 6 and 7 show how the low iteration cost of the

chord method can be combined with the small number of iterations of a

superlinearly convergent nonlinear iterative method. In the context of this

particular example, the cost of the Jacobian computation can be reduced by

analytic evaluation and reuse of data that have already been computed for the

evaluation of F . The reader is encouraged to do this in Exercise 5.7.21.

In the remaining examples, we plot, but do not tabulate, the iteration

statistics. In Fig. 5.2 we compare the Shamanskii method with m = 2 with

Newton™s method. The Shamanskii computation required under 6 million

¬‚oating-point operations, as compared with over 8.8 for Newton™s method.

However, for this problem, the simple chord method is the most e¬cient.

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.

89

NEWTON™S METHOD

0

10

-1

10

Relative Nonlinear Residual

-2

10

-3

10

-4

10

-5

10

-6

10

-7

10

0 1 2 3 4 5 6 7 8

Iterations

Fig. 5.1. Newton and chord methods for H-equation, c = .9.

0

10

-1

10

Relative Nonlinear Residual

-2

10

-3

10

-4

10

-5

10

-6

10

-7

10

0 0.5 1 1.5 2 2.5 3 3.5 4

Iterations

Fig. 5.2. Newton and Shamanskii method for H-equation, c = .9.

The hybrid approach in nsol with the default parameters would be the chord

method for this problem.

We now reconsider the H-equation with c = .9999. For this problem the

Jacobian at the solution is nearly singular [110]. Aside from c, all iteration

parameters are the same as in the previous example. While Newton™s method

converges in 7 iterations and requires 21 million ¬‚oating-point operations, the

chord method requires 188 iterations and 10.6 million ¬‚oating-point operations.

The q-factor for the chord method in this computation was over .96, indicating

very slow convergence. The hybrid method in nsol with the default parameters

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.

90 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

computes the Jacobian four times during the course of the iteration, converges

in 14 iterations, and requires 12 million ¬‚oating-point operations. In Fig. 5.3

we compare the Newton iteration with the hybrid.

0

10

-1

10

Relative Nonlinear Residual

-2

10

-3

10

-4

10

-5

10

-6

10

0 2 4 6 8 10 12 14

Iterations

Fig. 5.3. Newton and hybrid methods for H-equation c = .9999.

One should approach plots of convergence histories with some caution. The

theory for the chord method asserts that the error norms converge q-linearly to

zero. This implies only that the nonlinear residual norms { F (xn ) } converge

r-linearly to zero. While the plots indicate q-linear convergence of the nonlinear

residuals, the theory does not, in general, support this. In practice, however,

plots like those in this section are common.

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.

91

NEWTON™S METHOD

5.7. Exercises on Newton™s method

In some of the exercises here, and in the rest of this part of the book, you will

be asked to plot or tabulate iteration statistics. When asked to do this, for

each iteration tabulate the iteration counter, the norm of F , and the ratio

of F (xn ) / F (xn’1 ) for n ≥ 1. A better alternative in the MATLAB

environment is to use the semilogy command to plot the norms. When

one does this, one can visually inspect the plot to determine superlinear

convergence (concavity) without explicitly computing the ratios. Use the l∞

norm.

5.7.1. A function G is said to be H¨lder continuous with exponent ± in „¦ if

o

± for all x, y ∈ „¦. Show that if the Lipschitz

G(x) ’ G(y) ¤ K x ’ y

continuity condition on F in the standard assumptions is replaced by

H¨lder continuity with exponent ± > 0 that the Newton iterates converge

o

with q-order 1 + ±.

5.7.2. Can the performance of the Newton iteration be improved by a linear

change of variables? That is, for nonsingular N — N matrices A and

B, can the Newton iterates for F (x) = 0 and AF (Bx) = 0 show any

performance di¬erence when started at the same initial iterate? What

about the chord method?

5.7.3. Let ∇h F be given by De¬nition 5.4.1. Given the standard assumptions,

prove that the iteration given by

xn+1 = xn ’ (∇hn F (xn ))’1 F (xn )

is locally q-superlinearly convergent if hn ’ 0. When is it locally q-

quadratically convergent?

5.7.4. Suppose F (xn ) is replaced by ∇hn F (xn ) in the Shamanskii method.

Discuss how hn must be decreased as the iteration progresses to preserve

the q-order of m + 1 [174].

5.7.5. In what way is the convergence analysis of the secant method changed if

there are errors in the evaluation of f ?

5.7.6. Prove that the secant method converges r-superlinearly with r-order

√

( 5 + 1)/2. This is easy.

5.7.7. Show that the secant method converges q-superlinearly with q-order

√

( 5 + 1)/2.

5.7.8. The bisection method produces a sequence of intervals (xk , xk ) that r

l

contain a root of a function of one variable. Given (xk , xk ) with r

l

k )f (xk ) < 0 we replace one of the endpoints with y = (xk + xk )/2

f (xl r r

l

k+1 k+1 ) < 0. If f (y) = 0, we terminate the iteration.

to force f (xl )f (xr

Show that the sequence xk = (xk + xk )/2 is r-linearly convergent if f is

r

l

0 , x0 ) such that f (x0 )f (x0 ) < 0.

continuous and there exists (xl r r

l

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.

92 ITERATIVE METHODS FOR LINEAR AND NONLINEAR EQUATIONS

5.7.9. Write a program that solves single nonlinear equations with Newton™s

method, the chord method, and the secant method. For the secant

method, use x’1 = .99x0 . Apply your program to the following

function/initial iterate combinations, tabulate or plot iteration statistics,

and explain your results.

1. f (x) = cos(x) ’ x, x0 = .5

2. f (x) = tan’1 (x), x0 = 1

3. f (x) = sin(x), x0 = 3

4. f (x) = x2 , x0 = .5

5. f (x) = x2 + 1, x0 = 10

5.7.10. For scalar functions of one variable a quadratic model of f about a point

xc is

mc (x) = f (xc ) + f (xc )(x ’ xc ) + f (xc )(x ’ xc )2 /2.

Consider the iterative method that de¬nes x+ as a root of mc . What

problems might arise in such a method? Resolve these problems and

implement your algorithm. Test it on the problems in Exercise 5.7.9. In

what cases is this algorithm superior to Newton™s method? For hints and

generalizations of this idea to functions of several variables, see [170]. Can

the ideas in [170] be applied to some of the high-order methods discussed

in [190] for equations in one unknown?

5.7.11. Show that if f is a real function of one real variable, f is Lipschitz

continuous, and f (x— ) = f (x— ) = 0 but f (x— ) = 0 then the iteration

xn+1 = xn ’ 2f (xn )/f (xn )

converges locally q-quadratically to x— provided x0 is su¬ciently near to

x— , but not equal to x— [171].

5.7.12. Show that if f is a real function of one real variable, the standard

assumptions hold, f (x— ) = 0 and f (x) is Lipschitz continuous, then

the Newton iteration converges cubically (q-superlinear with q-order 3).

5.7.13. Show that if f is a real function of one real variable, the standard

assumptions hold, x0 is su¬ciently near x— , and f (x0 )f (x0 ) > 0, then

the Newton iteration converges monotonically to x— . What if f (x— ) = 0?

How can you extend this result if f (x— ) = f (x— ) = . . . f (k) (x— ) = 0 =

f (k+1) (x— )? See [131].

5.7.14. How would you estimate the q-order of a q-superlinearly convergent

sequence? Rigorously justify your answer.

5.7.15. If xn ’ x— q-superlinearly with q-order ± > 1, show that log( en ) is a

concave function of n in the sense that

log( en+1 ) ’ log( en )

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.

93

NEWTON™S METHOD

is a decreasing function of n for n su¬ciently large. Is this still the case

if the convergence is q-superlinear but there is no q-order? What if the

convergence is q-linear?

5.7.16. Show that the sequence on the right side of (5.20) is r-quadratically (but

not q-quadratically) convergent to 0 if β·γ < 1/2.

5.7.17. Typically (see [89], [184]) the computed LU factorization of F (x) is the

exact LU factorization of a nearby matrix.

LU = F (x) + E.

Use the theory in [89], [184] or another book on numerical linear algebra

and Theorem 5.4.1 to describe how this factorization error a¬ects the

convergence of the Newton iteration. How does the analysis for the QR

factorization di¬er?

5.7.18. Duplicate the results on the H-equation in § 5.6. Try di¬erent values of

N and c. How does the performance of the methods di¬er as N and c

change? Compare your results to the tabulated values of the H function

in [41] or [14]. Compare execution times in your environment. Tabulate

or plot iteration statistics. Use the MATLAB cputime command to see

how long the entire iteration took for each N, c combination and the

MATLAB flops command to get an idea for the number of ¬‚oating-

point operations required by the various methods.

5.7.19. Solve the H-equation with c = 1 and N = 100. Explain your results

(see [50] and [157]). Does the trick in Exercise 5.7.11 improve the

convergence? Try one of the approaches in [51], [49], or [118] for

acceleration of the convergence. You might also look at [90], [170], [75],

[35] and [155] for more discussion of this. Try the chord method and

compare the results to part 4 of Exercise 5.7.9 (see [52]).

5.7.20. Solve the H-equation with N = 100, c = .9 and c = .9999 with ¬xed-point

iteration. Compare timings and ¬‚op counts with Newton™s method, the

chord method, and the algorithm in nsol.

5.7.21. Compute by hand the Jacobian of F , the discretized nonlinearity from

the H equation in (5.22). Prove that F can be computed at a cost of