<< . .

. 14
( : 26)



. . >>

0

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

<< . .

. 14
( : 26)



. . >>