c19-6 (779621), страница 2
Текст из файла (страница 2)
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).• Compute the next approximation by (19.6.11).Let’s contrast the advantages and disadvantages of relaxation and the coarse-gridcorrection scheme. Consider the error vh expanded into a discrete Fourier series. Callthe components in the lower half of the frequency spectrum the smooth componentsand the high-frequency components the nonsmooth components. We have seen thatrelaxation becomes very slowly convergent in the limit h → 0, i.e., when there are alarge number of mesh points. The reason turns out to be that the smooth componentsare only slightly reduced in amplitude on each iteration.
However, many relaxationmethods reduce the amplitude of the nonsmooth components by large factors oneach iteration: They are good smoothing operators.For the two-grid iteration, on the other hand, components of the error withwavelengths <∼ 2H are not even representable on the coarse grid and so cannot bereduced to zero on this grid. But it is exactly these high-frequency components thatcan be reduced by relaxation on the fine grid! This leads us to combine the ideasof relaxation and coarse-grid correction:Two-Grid Iteration• Pre-smoothing: Compute ūh by applying ν1 ≥ 0 steps of a relaxationmethod to ueh .• Coarse-grid correction: As above, using ūh to give ūnewh .• Post-smoothing: Compute uenewby applying ν2 ≥ 0 steps of the relaxationhmethod to ūnewh .It is only a short step from the above two-grid method to a multigrid method.Instead of solving the coarse-grid defect equation (19.6.8) exactly, we can getan approximate solution of it by introducing an even coarser grid and using thetwo-grid iteration method.
If the convergence factor of the two-grid method issmall enough, we will need only a few steps of this iteration to get a good enoughapproximate solution. We denote the number of such iterations by γ. Obviouslywe can apply this idea recursively down to some coarsest grid. There the solutionis found easily, for example by direct matrix inversion or by iterating the relaxationscheme to convergence.One iteration of a multigrid method, from finest grid to coarser grids and backto finest grid again, is called a cycle.
The exact structure of a cycle depends onthe value of γ, the number of two-grid iterations at each intermediate stage. Thecase γ = 1 is called a V-cycle, while γ = 2 is called a W-cycle (see Figure 19.6.1).These are the most important cases in practice.Note that once more than two grids are involved, the pre-smoothing steps afterthe first one on the finest grid need an initial approximation for the error v. Thisshould be taken to be zero.87519.6 Multigrid Methods for Boundary Value ProblemsSS2-gridSSS3-gridSSSESSSSSSEγ=1SSEES4-gridSSSSSESSSESESEγ=2Figure 19.6.1. Structure of multigrid cycles.
S denotes smoothing, while E denotes exact solutionon the coarsest grid. Each descending line \ denotes restriction (R) and each ascending line / denotesprolongation (P ). The finest grid is at the top level of each diagram. For the V-cycles (γ = 1) the Estep is replaced by one 2-grid iteration each time the number of grid levels is increased by one. For theW-cycles (γ = 2), each E step gets replaced by two 2-grid iterations.where new values of u are used on the right-hand side as they become available. Theexact form of the Gauss-Seidel method depends on the ordering chosen for the meshpoints. For typical second-order elliptic equations like our model problem equation(19.0.3), as differenced in equation (19.0.8), it is usually best to use red-blackordering, making one pass through the mesh updating the “even” points (like the redsquares of a checkerboard) and another pass updating the “odd” points (the blacksquares). When quantities are more strongly coupled along one dimension thananother, one should relax a whole line along that dimension simultaneously.
Linerelaxation for nearest-neighbor coupling involves solving a tridiagonal system, andso is still efficient. Relaxing odd and even lines on successive passes is called zebrarelaxation and is usually preferred over simple line relaxation.Note that SOR should not be used as a smoothing operator. The overrelaxationdestroys the high-frequency smoothing that is so crucial for the multigrid method.A succint notation for the prolongation and restriction operators is to give theirsymbol. The symbol of P is found by considering vH to be 1 at some mesh point(x, y), zero elsewhere, and then asking for the values of PvH . The most popularprolongation operator is simple bilinear interpolation.
It gives nonzero values atthe 9 points (x, y), (x + h, y), . . . , (x − h, y − h), where the values are 1, 12 , . . . , 14 .Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).SE876Chapter 19.Partial Differential EquationsIts symbol is therefore14121412112141214(19.6.13)x,yThen the adjoint of P, denoted P † , is defined byhuH |P † vh iH = hPuH |vh ih(19.6.15)Now take P to be bilinear interpolation, and choose uH = 1 at (x, y), zero elsewhere.Set P † = R in (19.6.15) and H = 2h. You will find that(Rvh )(x,y) = 14 vh (x, y) + 18 vh (x + h, y) +116 vh (x+ h, y + h) + · · · (19.6.16)so that the symbol of R is1161811618141811618116(19.6.17)Note the simple rule: The symbol of R is 14 the transpose of the matrix defining thesymbol of P, equation (19.6.13).
This rule is general whenever R = P † and H = 2h.The particular choice of R in (19.6.17) is called full weighting. Another popularchoice for R is half weighting, “halfway” between full weighting and straightinjection. Its symbol is0 18 01 1 1(19.6.18)8 2 80180A similar notation can be used to describe the difference operator Lh .
Forexample, the standard differencing of the model problem, equation (19.0.6), isrepresented by the five-point difference star01 01 (19.6.19)Lh = 2 1 −4 1 h01 0Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The symbol of R is defined by considering vh to be defined everywhere on thefine grid, and then asking what is Rvh at (x, y) as a linear combination of thesevalues.
The simplest possible choice for R is straight injection, which means simplyfilling each coarse-grid point with the value from the corresponding fine-grid point.Its symbol is “[1].” However, difficulties can arise in practice with this choice. Itturns out that a safe choice for R is to make it the adjoint operator to P. To define theadjoint, define the scalar product of two grid functions uh and vh for mesh size h asXuh (x, y)vh (x, y)(19.6.14)huh |vh ih ≡ h219.6 Multigrid Methods for Boundary Value Problems877Full Multigrid AlgorithmSo far we have described multigrid as an iterative scheme, where one startswith some initial guess on the finest grid and carries out enough cycles (V-cycles,W-cycles,.















