Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 39
Текст из файла (страница 39)
. . , tN = 0, 0.1, 0.3, 0.6, 1. (In the early stages, iterate only long enough todetect quadratic convergence.)5.2-6 Solve the problem in Example 5.4 for n = 10 and G(x) = x5.*5.3FIXED-POINT ITERATION AND RELAXATION METHODS223*5.3 FIXED-POINT ITERATION AND RELAXATIONMETHODSNewton’s method and some of its variants discussed in Sec. 5.2 areexamples of fixed-point iteration. Here, one rewrites the equationf(ξ)(ξ) = 0into an equivalent one of the formξ = g(( ξ )and then, starting from some initial guess x (0), generates the sequencewhich, so one hopes, converges to the fixed point ξ of g.For example, Newton’s method is such a fixed-point iteration, with theiteration function g given byg(x) = x - f´(x)- 1f(x)More generally, the quasi-Newton methods use an iteration function of theform(5.12)g(x) = x - Cf(x)with C = C(x) some matrix. Relaxation, as discussed later in this section,provides a different idea for constructing iteration functions for solvingf(ξ) = 0 by fixed-point iteration.The analysis of fixed-point iteration for systems differs little from thatgiven for the case of one equation in Chap.
3, the only difference beingthat we now measure the size of the error ξ - x ( m ) in the mth iterate bynorms rather than absolute values.Theorem 5.1 Suppose the iteration function g maps some closed set Sinto itself, i.e., g(x) belongs to S if x does, and suppose further that g iscontractive on S, i.e.,||g(x) - g(y)|| < K||x - y||for all x and y in S and some K < 1. Then(i) g has a fixed point in S.(ii) If ξ is any fixed point of g in S, then fixed-point iteration startingwith any x ( 0 ) in S converges to ξ, i.e.,forsuch a sequence x ( m +1) = g(x( m )), m = 0, 1, 2, . .
. . More explicitly,(5.13)hence(5.14)224*SYSTEMS OF EQUATIONS AND UNCONSTRAINED OPTIMIZATIONThe assumptions ensure that we can start with any x(0) in S andcontinue the iteration x( m +1) = g(x( m )), m = 0, 1, 2, . . . , indefinitely, witheach x( m ) in S. Further, by an argument which goes beyond the level of thisbook (namely using the completeness of n-dimensional space), (i) follows.Finally, to get the estimate (5.13) and thereby (5.14), observe that(5.15)since g is contractive, hence, by the triangle inequality (4.33iii),orNow combine this inequality with (5.15) to get (5.13).Example 5.6 Newton’s method is fixed-point iteration with the iteration functiong(x) = x - f´(x)- 1 f(x).Thusf´(x)[g(x) - x] = -f(x)while, by (5.7) we find0 = f(ξ)(ξ) = f(x + ξ - x)assuming that f has continuous first and second partial derivatives.
Hence, substitutinghere - f´(x)[g(x) - x] for f(x), we get0 = f´(x)[ - (g(x) - x) + (ξξ - x) +orf´(x)[g(x) - ξ] =This says thatfor some constant c.If now f´(ξ)(ξ) is invertible, then, since f´(x) is continuous by assumption, we can finda positive δ and an M so that f´(x) -1 exists for all x within δ of ξ and has a matrix normno bigger than M.
But then, choosing ε to be the smaller of δ and (M c)-1, we have, forall x in the closed setthat f´(x)-1 exists (hence g(x) is defined) andThus g maps the closed set S into itself. Further, if || ξ - x|| < ε, then K = Mc||ξ - x||< 1, hence ξ is an attracting fixed point of g, and iteration starting with any x (0) withinless than ε of ξ will converge to ξ.As a further illustration, we now consider the solution of the linearsystem(5.16)Aξ=b*5.3FIXED-POINT ITERATION AND RELAXATION METHODS225by fixed-point iteration. Such iteration schemes can all be based on thenotion of approximate inverse.
By this we mean any matrix C for which(5.17)||I - CA|| < 1in some matrix norm.Lemma 5.1 If C is an approximate inverse for the matrix A, i.e., if||I - CA|| < 1 in some matrix norm, then both C and A are invertible.Indeed, if C or A were not invertible, then neither would the matrixCA be (see Exercise 4.1-8). By Theorem 4.4, we could then find x0 sothat CAx = 0.
But thenwhich is nonsense.In particular, (5.16) has exactly one solution if A has an approximateinverse.Corresponding to an approximate inverse C for A, we consider theiteration functiong(x) = Cb + (I - CA)x= x + C(b - Ax)Note that this iteration function is of quasi-Newton type, i.e., of the formg(x) = x - Cf(x), if we take f(x) = Ax - b.
Also,g(x) - g(y) = Cb + (I - CA) x - [Cb + (I - CA)y]= (I - CA)(x - y)Consequently,||g(x) - g(y)|| < || I - CA|| ||xshowing g to be contractive, with-y||(5.18)K = ||I - CA|| < 1Therefore, fixed point iterationx ( m +1) = x ( m ) + C(b - Ax( m ) )m = 0, 1, 2, . ..starting from any x , will converge to the unique solution ξ of (5.16), withthe error at each step reduced by at least a factor of K = ||I - CA||.(0)Example 5.7 Suppose the matrix A is strictly row diagonally dominant, i.e.,Let D - diag(a11, a22, .
. . ,ann) be the diagonal of A. Then226*SYSTEMS OF EQUATIONS AND UNCONSTRAINED OPTIMIZATIONTable 5.1showing that D is then an approximate inverse for A. The corresponding iterationschemex (m+1)= x( m)+ D - l (b - Ax( m ) )m = 0, 1, 2, . . .(5.19)is Jacobi iteration. Note that x( m +1) can be obtained from x( m ) by solving, for each i, thei th equation for the ith unknown, giving all the other unknowns their current values. Informulas,For the particular linear system1 0x1 + x2 + x3 = 12x 1 + 10x2 + x3 = 12x1 + x2 + 10x3 = 12(0)Jacobi iteration starting with x = 0 produces the vectors x(l), x(2), . .
. ,x(6) listedabove in Table 5.1. The sequence seems to converge nicely to the solution [1 1 l] T ofthe system. For this example,so that we would expect a reduction in error by at least a factor of 0.2 per step, which isborne out by the numbers in Table 5.1.It is, of course, easy in principle to find an approximate inverse C forA. For example, C = A - 1 would do, and the corresponding iterationwould converge in one step. But, the point of using iteration for solvingAξ = b in the first place is that one might obtain an approximate solutionof acceptable accuracy much faster by iteration than by solving Aξ = bdirectly. For this, it is important to choose C so that we can calculate thevector Cr for any particular r with much less work than it would take tocalculate the vector A- 1 r.
Typically, one chooses C as the inverse of adiagonal matrix (as in Jacobi iteration), or the inverse of a triangularmatrix (as in Gauss-Seidel iteration discussed below), or as the inverse ofthe product of two triangular matrices (as in the iterative improvementalgorithm 4.5), or even as the inverse of a tridiagonal matrix, etc.*5.3FIXED-POINT ITERATION AND RELAXATION METHODS227Algorithm 5.5: Fixed-point iteration for linear systems Given the linearsystem Aξ = b of order n.Pick a matrix C of order n such that(i) For given r, the vector Cr is “easily” calculated(ii) In some matrix norm, ||I - CA|| < 1Pick an n-vector x(0), for example, x(0) = 0For m =0, 1, 2, .
. . , until satisfied, do:In the absence of round-off error, the resulting sequencex (0), x(1), x(2), . . . converges to the solution of the given linear system.As in Chap. 3, we employ here the phrase “until satisfied” to stress theincompleteness of the description given. To complete the algorithm, onehas to specify precise termination criteria. Typical criteria are:Terminate if (a):for some prescribed εor if (b):or if (c):for some given MThe last criterion should always be present in any program implementingthe algorithm.
We repeat the warning first voiced in Sec. 1.6: The fact that||x ( m ) - x( m - 1 ) || < εdoes not imply that||x( m ) - ξ|| < εBut we do know from (5.13) and (5.18) that(5.20)with K = ||I - CA||.To give an example, we found for the Jacobi iteration in Example 5.7 thatTherefore, (5.20) gives the estimateIn fact, ||ξ - x(6)|| = 0.000064, so that the error is overestimated by only 50 percent.Unfortunately, it is usually difficult to obtain good estimates for ||I - CA ||, or else theestimate for ||I - CA|| is so close to 1 as to make the denominator 1 - K in (5.20)excessively small and the resulting bound on ||ξ - x( m )|| useless.It should be pointed out that C may be an approximate inverse for Aeven though| |I - CA|| > 1228*SYSTEMS OF EQUATIONS AND UNCONSTRAINED OPTIMIZATIONfor some particular matrix norm.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
