Heath - Scientific Computing (523150), страница 92
Текст из файла (страница 92)
Unfortunately, the convergence rate of the Jacobi method is usually unacceptablyslow.11.5.3Gauss-Seidel MethodOne reason for the slow convergence of the Jacobi method is that it does not make use ofthe latest information available: new component values are used only after the entire sweephas been completed. The Gauss-Seidel method remedies this drawback by using each newcomponent of the solution as soon as it has been computed:(k+1)xi=bi −(k+1)j<i aij xjPaii−(k)j>i aij xjP,i = 1, . . . , n.In the same notation as in Section 11.5.2, the Gauss-Seidel method corresponds to thesplittingM = D + L, N = −Uand can be written in matrix terms asx(k+1) = D −1 (b − Lx(k+1) − U x(k) )= (D + L)−1 (b − U x(k) ).In addition to faster convergence, another benefit of the Gauss-Seidel method is that duplicate storage is not needed for the vector x, since the newly computed component values can344CHAPTER 11.
PARTIAL DIFFERENTIAL EQUATIONSoverwrite the old ones immediately (a programmer would have invented this method in thefirst place because of its more natural and convenient implementation). On the other hand,the updating of the unknowns must now be done successively, in contrast to the Jacobimethod, in which the unknowns can be updated in any order or even simultaneously. Thelatter feature may make Jacobi preferable on a parallel computer.To illustrate the use of the Gauss-Seidel method, if we apply it to solve the system offinite difference equations for the Laplace equation in Example 11.4, we get(k+1)(k+1)ui,j=(k+1)(k)(k)ui−1,j + ui,j−1 + ui+1,j + ui,j+14,assuming that we sweep from left to right and bottom to top in the grid.
Thus, we againaverage the solution values at the four surrounding grid points but always use new component values as soon as they become available rather than waiting until the current iterationhas been completed.The Gauss-Seidel method does not always converge, but it is guaranteed to convergeunder conditions that are often satisfied in practice and that are somewhat weaker thanthose for the Jacobi method (e.g., if the matrix is symmetric and positive definite).
Althoughthe Gauss-Seidel method converges more rapidly than the Jacobi method, it is often stilltoo slow to be practical.11.5.4Successive Over-RelaxationThe convergence rate of the Gauss-Seidel method can be accelerated by a technique calledsuccessive over-relaxation (SOR), which in effect uses the step to the next Gauss-Seideliterate as a search direction, but with a fixed search parameter denoted by ω. Starting with(k+1)x(k) , we first compute the next iterate that would be given by Gauss-Seidel, xGS , theninstead take the next iterate to be(k+1)x(k+1) = x(k) + ω(xGS− x(k) ).Equivalently, we can think of this scheme as taking a weighted average of the current iterateand the next Gauss-Seidel iterate:(k+1)x(k+1) = (1 − ω)x(k) + ωxGS .In either case, ω is a fixed relaxation parameter chosen to accelerate convergence.
Avalue ω > 1 gives over -relaxation, whereas ω < 1 gives under -relaxation (ω = 1 simply givesthe Gauss-Seidel method). We always have 0 < ω < 2 (otherwise the method diverges),but choosing a specific value of ω to attain the best possible convergence rate is a difficultproblem in general and is the subject of an elaborate theory for special classes of matrices.In the same notation as in Section 11.5.2, the SOR method corresponds to the splittingM = D + ωL,N = (1 − ω)D − ωU ,and can be written in matrix terms asx(k+1) = x(k) + ω[D −1 (b − Lx(k+1) − U x(k) ) − x(k) ]= (D + ωL)−1 [(1 − ω)D − ωU ]x(k) + ω(D + ωL)−1 b.11.5.
ITERATIVE METHODS FOR LINEAR SYSTEMS345Like the Gauss-Seidel method, the SOR method makes repeated forward sweeps throughthe unknowns, updating them successively. A variant of SOR, known as SSOR (symmetricSOR), alternates forward and backward sweeps through the unknowns. SSOR is not necessarily faster than SOR (indeed SSOR is often slower), but it has the theoretical advantagethat its iteration matrix, G = M −1 N , which is too complicated to express here, is similarto a symmetric matrix when A is symmetric (which is not true of the iteration matrix forSOR).
For example, this makes SSOR useful as a preconditioner (see Section 11.5.5).11.5.5Conjugate Gradient MethodWe now turn from stationary iterative methods to methods based on optimization. If A isan n × n symmetric positive definite matrix, then the quadratic functionφ(x) = 12 xT Ax − xT battains a minimum precisely when Ax = b. Thus, we can apply any of the optimizationmethods discussed in Section 6.3 to obtain a solution to the corresponding linear system.Recall from Section 6.3 that most multidimensional optimization methods progress from oneiteration to the next by performing a one-dimensional search along some search directionsk , so thatxk+1 = xk + αsk ,where α is a search parameter chosen to minimize the objective function φ(xk + αsk ) alongsk .We note some special features of such a quadratic optimization problem. First, thenegative gradient is simply the residual vector:−∇φ(x) = b − Ax = r.Second, for any search direction sk , we need not perform a line search, because the optimalchoice for α can be determined analytically.
Specifically, the minimum over α occurs whenthe new residual is orthogonal to the search direction:0=dddTφ(xk+1 ) = ∇φ(xk+1 )Txk+1 = (Axk+1 − b)T ( (xk + αsk )) = −rk+1sk .dαdαdαSince the new residual can be expressed in terms of the old residual and the search direction,rk+1 = b − Axk+1 = b − A(xk + αsk ) = (b − Axk ) − αAsk = rk − αAsk ,we can thus solve forα=rkT sk.sTk AskIf we take advantage of these properties in the algorithm of Section 6.3.6, we obtainthe conjugate gradient (CG) method for solving symmetric positive definite linear systems.Starting with an initial guess x0 and taking s0 = r0 = b − Ax0 , the following steps arerepeated for k = 0, 1, .
. . until convergence:1. αk = rkT rk /sTk Ask .3462.3.4.5.CHAPTER 11. PARTIAL DIFFERENTIAL EQUATIONSxk+1 = xk + αk sk .rk+1 = rk − αk Ask .T rTβk+1 = rk+1k+1 /rk rk .sk+1 = rk+1 + βk+1 sk .Each iteration of the algorithm requires only a single matrix-vector multiplication, Ask ,plus a small number of inner products. The storage requirements are also very modest,since the vectors x, r, and s can be overwritten.Although the foregoing algorithm is not terribly difficult to derive, we content ourselveshere with the following intuitive motivation.
The features noted earlier for the quadraticoptimization problem would make it extremely easy to apply the steepest descent method,using the negative gradient—in this case the residual—as search direction at each iteration.Unfortunately, we have already observed that its convergence rate is often very poor owingto repeated searches in the same directions (zigzagging).
We could avoid this repetition byorthogonalizing each new search direction against all of the previous ones (see Section 3.4.6),leaving only components in “new” directions, but this would appear to be prohibitivelyexpensive computationally and would also require excessive storage to save all of the searchdirections. However, if instead of using the standard inner product we make the searchdirections mutually A-orthogonal (vectors y and z are A-orthogonal if y T Az = 0), orconjugate, then it can be shown that the successive A-orthogonal search directions satisfy athree-term recurrence (this is the role played by β in the algorithm).
This short recurrencemakes the computation very cheap, and, most important, it means that we do not need tosave all of the previous gradients, only the most recent two, which makes a huge differencein storage requirements.In addition to the other special properties already mentioned, it turns out that in thequadratic case the residual at each step is minimal (with respect to the norm induced by A)over the space spanned by the search directions generated so far. Since the search directionsare A-orthogonal, and hence linearly independent, this property implies that after at mostn steps the solution is exact, because the n search directions must span the whole space.Thus, in theory, the conjugate gradient method is direct, but in practice rounding errorcauses a loss of orthogonality, which spoils this finite termination property.
As a result,the conjugate gradient method is usually used in an iterative manner and halted whenthe residual, or some other measure of error, is sufficiently small. In practice, the methodoften converges in far fewer than n iterations. We will consider its convergence rate inSection 11.5.6.Although it is a significant improvement over steepest descent, the conjugate gradientalgorithm can still converge very slowly if the matrix A is ill-conditioned. Convergence canoften be substantially accelerated by preconditioning, which can be thought of as implicitlymultiplying A by M −1 , where M is a matrix for which systems of the form M z = y areeasily solved, and whose inverse approximates that of A, so that M −1 A is relatively wellconditioned.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
