882 Chapter 19.Partial Differential Equations The defect dh vanishes identically at all black mesh points after a red-black Gauss-Seidel step.Thus d=Rdh for half-weighting reduces to simply copying half the defect from the fine grid to the corresponding coarse-grid point.The calls to resid followed by rstrct in the first part of the V-cycle can be replaced by a routine that loops only over the coarse grid, filling it with half the defect. Similarly,the quantity unew=h+Pg need not be computed at red mesh points,since they will immediately be redefined in the subsequent Gauss-Seidel sweep.This means that addint need only loop over black points. You can speed up relax in several ways.First,you can have a special form when the initial guess is zero,and omit the routine fill0.Next, B you can store h2fh on the various grids and save a multiplication.Finally, it is possible to save an addition in the Gauss-Seidel formula by rewriting Cam it with intermediate variables. On typical problems,mglin with ncycle =1 will return a solution with the iteration error bigger than the truncation error for the given size of h. RECIPES I To knock the error down to the size of the truncation error,you have to set ncycle =2 or,more cheaply,npre =2.A more efficient way turns 9 out to be to use a higher-order P in(19.6.20)than the linear interpolation used in the V-cycle. Implementing all the above features typically gives up to a factor of two improvement in execution time and is certainly worthwhile in a production code. 芒a 9 Nonlinear Multigrid:The FAS Algorithm Now turn to solving a nonlinear elliptic equation,which we write symbolically as 61 C(u)=0 (19.6.21) Any explicit source term has been moved to the left-hand side.Suppose equation (19.6.21) is suitably discretized: Ch(uh)=0 (19.6.22) We will see below that in the multigrid algorithm we will have to consider equations where a 方03 10621 nonzero right-hand side is generated during the course of the solution: Numerica Ch(uh)=fh (19.6.23) 431 One way of solving nonlinear problems with multigrid is to use Newton's method,which Recipes produces linear equations for the correction term at each iteration.We can then use linear multigrid to solve these equations.A great strength of the multigrid idea,however,is that it can be applied directly to nonlinear problems.All we need is a suitable nonlinear relaxation 拿 method to smooth the errors,plus a procedure for approximating corrections on coarser grids. This direct approach is Brandt's Full Approximation Storage Algorithm(FAS).No nonlinear equations need be solved,except perhaps on the coarsest grid. To develop the nonlinear algorithm,suppose we have a relaxation procedure that can smooth the residual vector as we did in the linear case.Then we can seek a smooth correction Uh to solve (19.6.23): Ch(uh +Uh)=fh (19.6.24) To find vh,note that Ch(un +Uh)-Ch(un)=fh-Ch(un) (19.6.25) =-dh882 Chapter 19. Partial Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). • The defect dh vanishes identically at all black mesh points after a red-black Gauss-Seidel step. Thus dH = Rdh for half-weighting reduces to simply copying half the defect from the fine grid to the corresponding coarse-grid point. The calls to resid followed by rstrct in the first part of the V-cycle can be replaced by a routine that loops only over the coarse grid, filling it with half the defect. • Similarly, the quantity u
new h = u
h + Pv
H need not be computed at red mesh points, since they will immediately be redefined in the subsequent Gauss-Seidel sweep. This means that addint need only loop over black points. • You can speed up relax in several ways. First, you can have a special form when the initial guess is zero, and omit the routine fill0. Next, you can store h2fh on the various grids and save a multiplication. Finally, it is possible to save an addition in the Gauss-Seidel formula by rewriting it with intermediate variables. • On typical problems, mglin with ncycle = 1 will return a solution with the iteration error bigger than the truncation error for the given size of h. To knock the error down to the size of the truncation error, you have to set ncycle = 2 or, more cheaply, npre = 2. A more efficient way turns out to be to use a higher-order P in (19.6.20) than the linear interpolation used in the V-cycle. Implementing all the above features typically gives up to a factor of two improvement in execution time and is certainly worthwhile in a production code. Nonlinear Multigrid: The FAS Algorithm Now turn to solving a nonlinear elliptic equation, which we write symbolically as L(u)=0 (19.6.21) Any explicit source term has been moved to the left-hand side. Suppose equation (19.6.21) is suitably discretized: Lh(uh)=0 (19.6.22) We will see below that in the multigrid algorithm we will have to consider equations where a nonzero right-hand side is generated during the course of the solution: Lh(uh) = fh (19.6.23) One way of solving nonlinear problems with multigrid is to use Newton’s method, which produces linear equations for the correction term at each iteration. We can then use linear multigrid to solve these equations. A great strength of the multigrid idea, however, is that it can be applied directly to nonlinear problems. All we need is a suitable nonlinear relaxation method to smooth the errors, plus a procedure for approximating corrections on coarser grids. This direct approach is Brandt’s Full Approximation Storage Algorithm (FAS). No nonlinear equations need be solved, except perhaps on the coarsest grid. To develop the nonlinear algorithm, suppose we have a relaxation procedure that can smooth the residual vector as we did in the linear case. Then we can seek a smooth correction vh to solve (19.6.23): Lh(u
h + vh) = fh (19.6.24) To find vh, note that Lh(u
h + vh) − Lh(u
h) = fh − Lh(u
h) = −dh (19.6.25)