Heath - Scientific Computing (523150), страница 91
Текст из файла (страница 91)
×.............................. ×............... ×.......................××××××××× ××××× ×××××××× × × ×........................................× .......× ........× ........× .......× ............................. ×.............................. ×............... ×.......................×××××× ××+×× +×Amesh×+ ×+ + ×× × ×......................................................................× ..............LFigure 11.7: Finite difference mesh reordered by minimum degree, with nonzero patterns of correspondingly permuted sparse matrix A and its Cholesky factor L.to any node in the other, and hence no fill can occur in either piece as a consequence ofthe elimination of a node in the other.
In other words, dissection induces blocks of zeros inthe matrix (indicated by the squares in Fig. 11.8) that are automatically preserved duringfactorization, thereby limiting fill. The recursive nature of the algorithm can be seen in thehierarchical block structure of the matrix, which would involve many more levels in a largerproblem................................... 4 ............................................ 6 ............................................ 5 ...............................................................................................................................................................................................................7.......8.......9..............................................................................................................................................................................................1.......3.......2....................................
× 2 ×....... 2 × ×........ × × ×.... .................................... ......... ......... ......... ......... ....... .............................................. ×......×......×...........mesh...................................................................................................×× 2 × ×2 × ×× × ×××× ×××××××..................× ......................× .......................× .......× ..............A............... ×....... 2 ×........ × × ×.... .................................... ......... .........
......... ......... ....... .............................................. ×+......×......× +...........×2 ×× × ××+ ×× × ×× + + ×......................................................................× ..............LFigure 11.8: Finite difference mesh reordered by nested dissection, with nonzero patterns of correspondingly permuted sparse matrix A and its Cholesky factor L.Sparse factorization methods are accurate, reliable, and robust. They are the methodsof choice for one-dimensional problems and are usually competitive for two-dimensionalproblems, but they can be prohibitively expensive in both work and storage for very largethree-dimensional problems. We will see that iterative methods provide a viable alternativein these cases.11.4.2Fast Direct MethodsFor certain types of PDEs, special techniques can be used to solve the resulting discretizedlinear system much faster than would be expected.
For example, for certain elliptic boundary value problems having constant coefficients and simple boundaries (e.g., the Poissonequation on a rectangular domain), the fast Fourier transform, or FFT (see Chapter 12),can be used to compute the solution to the discrete system very efficiently, provided thatthe number of mesh points in each dimension is a power of two. This technique is the basis11.5.
ITERATIVE METHODS FOR LINEAR SYSTEMS341for several “fast Poisson solver” software packages. For a problem with n mesh points,such a fast Poisson solver computes the solution in O(n log2 n) operations, which is nearlyoptimal since the cost of simply writing the output is O(n).Somewhat more generally, for separable elliptic PDEs the method of cyclic reductionpermits similarly fast solutions. Cyclic reduction is a divide-and-conquer technique in whichthe even-numbered equations in the systems are solved in terms of the odd-numbered ones,and so on recursively until reaching the bottom of the recursion, where single equationscan be solved trivially.
This idea obviously works best when the order of the system is apower of two, but it can be adapted to handle systems of arbitrary order. These ideas—FFT and cyclic reduction—can be combined, for example using FFT in one dimension andcyclic reduction in the other. A more subtle combination results in the FACR (Fourieranalysis/cyclic reduction) method, which is even faster than either the FFT method orthe cyclic reduction method alone.
The computational complexity of the FACR method isO(n log log n), which is effectively optimal, since log log is essentially constant for problemsof any reasonable size.11.5Iterative Methods for Linear SystemsIterative methods for solving linear systems begin with an initial estimate for the solutionand successively improve it until the solution is as accurate as desired.
In theory, an infinitenumber of iterations might be required to converge to the exact solution, but in practicethe iteration terminates when some norm of the residual kb − Axk, or some other measureof error, is as small as desired.11.5.1Stationary Iterative MethodsPerhaps the simplest type of iterative method for solving Ax = b has the formxk+1 = Gxk + c,where the matrix G and vector c are chosen so that a fixed point of the equation x = Gx+cis a solution to Ax = b. Such a method is said to be stationary if G and c are constantover all iterations.One way to obtain a suitable matrix G is by a splitting, in which the matrix A is writtenasA = M − N,with M nonsingular. We can then take G = M −1 N and c = M −1 b, so that the iterationscheme becomesxk+1 = M −1 N xk + M −1 b,which is implemented asM xk+1 = N xk + b(i.e., we solve a linear system with matrix M at each iteration).Formally, this splitting scheme is a fixed-point iteration with iteration functiong(x) = M −1 N x + M −1 b,342CHAPTER 11.
PARTIAL DIFFERENTIAL EQUATIONSwhose Jacobian matrix isG(x) = M −1 N .Thus, the iteration scheme is convergent ifρ(G) = ρ(M −1 N ) < 1,and the smaller the spectral radius, the faster the convergence rate (see Section 5.3.1).For rapid convergence, we should choose M and N so that ρ(M −1 N ) is as small aspossible. There is a trade-off, however, as the cost per iteration is determined by the costof solving a linear system with matrix M .
As an extreme example, if M = A, then thescheme converges in a single iteration (i.e., we have a direct method), but that one iterationmay be prohibitively expensive. In practice, M is chosen to approximate A in some sense,but is usually constrained to have some simple form, such as diagonal or triangular, so thatthe linear system at each iteration is easy to solve.Example 11.6 Iterative Refinement.
We have already seen one example of a stationaryiterative method, namely, iterative refinement of a solution already computed by Gaussianelimination (see Section 2.4.3). Forward- and back-substitution using the LU factorizationin effect provide an approximation, call it B −1 , to the inverse of A (i.e., for any right-handside vector y, the solution B −1 y can be computed by forward- and back-substitution usingthe LU factors already computed). Iterative refinement then has the formxk+1 = xk + B −1 (b − Axk ),which can be rewrittenxk+1 = (I − B −1 A)xk + B −1 b.Thus, we see that iterative refinement is a stationary iterative method with G = I − B −1 Aand c = B −1 b. The scheme therefore converges if ρ(I − B −1 A) < 1, which should bethe case if B −1 is a good approximation to A−1 , such as the use of forward- and backsubstitution with the LU factors obtained by Gaussian elimination with partial pivoting.Indeed, the convergence condition may be satisfied even by a rather loose approximation tothe inverse.
For example, iterative refinement can sometimes be used to stabilize “fast butrisky” algorithms.11.5.2Jacobi MethodIn the matrix splitting A = M − N , the simplest choice for M is diagonal, specifically thediagonal of A. Let D be a diagonal matrix with the same diagonal entries as A, and let Land U be the strict lower and upper triangular portions of A, respectively, so thatM = D,N = −(L + U )gives a splitting of A. If A has no zero diagonal entries, so that D is nonsingular, we obtainthe iterative scheme known as the Jacobi method:x(k+1) = D −1 (b − (L + U )x(k) ).11.5. ITERATIVE METHODS FOR LINEAR SYSTEMS343(We use parenthesized superscripts for the iteration index when we need to reserve subscriptsto refer to individual components of a vector.) Rewriting this scheme componentwise, wesee that, beginning with an initial guess x(0) , the Jacobi method computes the next iterateby solving for each component of x in terms of the others:(k+1)xi=bi −(k)j6=i aij xjPaii,i = 1, .
. . , n.Note that the Jacobi method requires double storage for the vector x because all of theold component values are needed throughout the sweep, and therefore the new componentvalues cannot overwrite them until the sweep has been completed.To illustrate the use of the Jacobi method, if we apply it to solve the system of finitedifference equations for the Laplace equation in Example 11.4, we get(k)(k+1)ui,j=(k)(k)(k)ui−1,j + ui,j−1 + ui+1,j + ui,j+14,which means that each new approximate solution at a given grid point is simply the averageof the previous solution components at the four surrounding grid points.
In this sense,solving the elliptic problem by an iterative method adds a timelike dimension (analogous toa parabolic problem, in this case the heat equation) in which the initial solution “diffuses”until a steady state is reached at the final solution.The Jacobi method does not always converge, but it is guaranteed to converge underconditions that are often satisfied in practice (e.g., if the matrix is diagonally dominant byrows).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
