19.6 Multigrid Methods for Boundary Value Problems 871 standard tridiagonal algorithm.Given u",one solves(19.5.36)for u+1/2,substitutes on the right-hand side of(19.5.37),and then solves for u+.The key question is how to choose the iteration parameter r,the analog of a choice of timestep for an initial value problem. As usual,the goal is to minimize the spectral radius of the iteration matrix. Although it is beyond our scope to go into details here,it turns out that,for the optimal choice of r,the ADI method has the same rate of convergence as SOR. The individual iteration steps in the ADI method are much more complicated than in SOR,so the ADI method would appear to be inferior.This is in fact true if we choose the same parameter r for every iteration step.However,it is possible to 81 choose a different r for each step.If this is done optimally,then ADI is generally more efficient than SOR.We refer you to the literature [1-4]for details. Our reason for not fully implementing ADI here is that,in most applications, it has been superseded by the multigrid methods described in the next section.Our advice is to use SOR for trivial problems (e.g.,20 x 20),or for solving a larger problem once only,where ease of programming outweighs expense of computer time.Occasionally,the sparse matrix methods of 82.7 are useful for solving a set of difference equations directly.For production solution of large elliptic problems, however,multigrid is now almost always the method of choice. 9 CITED REFERENCES AND FURTHER READING: Hockney,R.W.,and Eastwood,J.W.1981,Computer Simulation Using Particles (New York: McGraw-Hill),Chapter 6. 85总6 Young,D.M.1971,Iterative Solution of Large Linear Systems (New York:Academic Press).[1] Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). 558.3-8.6.2] Varga,R.S.1962,Matrix /terative Analysis (Englewood Cliffs,NJ:Prentice-Hall).[3] Spanier,J.1967,in Mathematical Methods for Digital Computers,Volume 2(New York:Wiley). Chapter 11.[4] 19.6 Multigrid Methods for Boundary Value Problems 、之花 61 Numerical 10.621 43108 Recipes Practical multigrid methods were first introduced in the 1970s by Brandt.These (outside methods can solve elliptic PDEs discretized on N grid points in O(N)operations North The "rapid"direct elliptic solvers discussed in $19.4 solve special kinds of elliptic equations in O(N log N)operations.The numerical coefficients in these estimates are such that multigrid methods are comparable to the rapid methods in execution speed.Unlike the rapid methods,however,the multigrid methods can solve general elliptic equations with nonconstant coefficients with hardly any loss in efficiency. Even nonlinear equations can be solved with comparable speed. Unfortunately there is not a single multigrid algorithm that solves all elliptic problems.Rather there is a multigrid technique that provides the framework for solving these problems.You have to adjust the various components of the algorithm within this framework to solve your specific problem.We can only give a brief
19.6 Multigrid Methods for Boundary Value Problems 871 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 machinereadable 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). standard tridiagonal algorithm. Given un, one solves (19.5.36) for un+1/2, substitutes on the right-hand side of (19.5.37), and then solves for u n+1. The key question is how to choose the iteration parameter r, the analog of a choice of timestep for an initial value problem. As usual, the goal is to minimize the spectral radius of the iteration matrix. Although it is beyond our scope to go into details here, it turns out that, for the optimal choice of r, the ADI method has the same rate of convergence as SOR. The individual iteration steps in the ADI method are much more complicated than in SOR, so the ADI method would appear to be inferior. This is in fact true if we choose the same parameter r for every iteration step. However, it is possible to choose a different r for each step. If this is done optimally, then ADI is generally more efficient than SOR. We refer you to the literature [1-4] for details. Our reason for not fully implementing ADI here is that, in most applications, it has been superseded by the multigrid methods described in the next section. Our advice is to use SOR for trivial problems (e.g., 20 × 20), or for solving a larger problem once only, where ease of programming outweighs expense of computer time. Occasionally, the sparse matrix methods of §2.7 are useful for solving a set of difference equations directly. For production solution of large elliptic problems, however, multigrid is now almost always the method of choice. CITED REFERENCES AND FURTHER READING: Hockney, R.W., and Eastwood, J.W. 1981, Computer Simulation Using Particles (New York: McGraw-Hill), Chapter 6. Young, D.M. 1971, Iterative Solution of Large Linear Systems (New York: Academic Press). [1] Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §§8.3–8.6. [2] Varga, R.S. 1962, Matrix Iterative Analysis (Englewood Cliffs, NJ: Prentice-Hall). [3] Spanier, J. 1967, in Mathematical Methods for Digital Computers, Volume 2 (New York: Wiley), Chapter 11. [4] 19.6 Multigrid Methods for Boundary Value Problems Practical multigrid methods were first introduced in the 1970s by Brandt. These methods can solve elliptic PDEs discretized on N grid points in O(N) operations. The “rapid” direct elliptic solvers discussed in §19.4 solve special kinds of elliptic equations in O(N log N) operations. The numerical coefficients in these estimates are such that multigrid methods are comparable to the rapid methods in execution speed. Unlike the rapid methods, however, the multigrid methods can solve general elliptic equations with nonconstant coefficients with hardly any loss in efficiency. Even nonlinear equations can be solved with comparable speed. Unfortunately there is not a single multigrid algorithm that solves all elliptic problems. Rather there is a multigrid technique that provides the framework for solving these problems. You have to adjust the various components of the algorithm within this framework to solve your specific problem. We can only give a brief
872 Chapter 19.Partial Differential Equations introduction to the subject here.In particular,we will give two sample multigrid routines,one linear and one nonlinear.By following these prototypes and by perusing the references [1-41,you should be able to develop routines to solve your own problems. There are two related,but distinct,approaches to the use ofmultigrid techniques. The first,termed"the multigrid method,"is a means for speeding up the convergence of a traditional relaxation method,as defined by you on a grid of pre-specified fineness.In this case,you need define your problem(e.g.,evaluate its source terms) only on this grid.Other,coarser,grids defined by the method can be viewed as temporary computational adjuncts. The second approach,termed(perhaps confusingly)"the full multigrid(FMG) method,"requires you to be able to define your problem on grids of various sizes (generally by discretizing the same underlying PDE into different-sized sets of finite- difference equations).In this approach,the method obtains successive solutions on 不 finer and finer grids.You can stop the solution either at a pre-specified fineness,or you can monitor the truncation error due to the discretization,quitting only when it is tolerably small. RECIPESI In this section we will first discuss the"multigrid method,"then use the concepts 器3d 2 developed to introduce the FMG method.The latter algorithm is the one that we implement in the accompanying programs. From One-Grid,through Two-Grid,to Multigrid The key idea of the multigrid method can be understood by considering the simplest case of a two-grid method.Suppose we are trying to solve the linear elliptic problem 6 Cu=f (19.6.1) where C is some linear elliptic operator and f is the source term.Discretize equation (19.6.1)on a uniform grid with mesh size h.Write the resulting set of linear algebraic equations as Chuh fh (19.6.2) 多，aoao Numerica 10.621 43106 Let ih denote some approximate solution to equation(19.6.2).We will use the symbol uh to denote the exact solution to the difference equations(19.6.2).Then the error in uh or the correction is Uh th uh (19.6.3 The residual or defect is dh=Chuh-fh (19.6.4) (Beware:some authors define residual as minus the defect,and there is not universal agreement about which of these two quantities 19.6.4 defines.)Since Ch is linear, the error satisfies ChUh =-dh (19.6.5)
872 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 machinereadable 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). introduction to the subject here. In particular, we will give two sample multigrid routines, one linear and one nonlinear. By following these prototypes and by perusing the references [1-4], you should be able to develop routines to solve your own problems. There are two related, but distinct, approaches to the use of multigrid techniques. The first, termed “the multigrid method,” is a means for speeding up the convergence of a traditional relaxation method, as defined by you on a grid of pre-specified fineness. In this case, you need define your problem (e.g., evaluate its source terms) only on this grid. Other, coarser, grids defined by the method can be viewed as temporary computational adjuncts. The second approach, termed (perhaps confusingly) “the full multigrid (FMG) method,” requires you to be able to define your problem on grids of various sizes (generally by discretizing the same underlying PDE into different-sized sets of finitedifference equations). In this approach, the method obtains successive solutions on finer and finer grids. You can stop the solution either at a pre-specified fineness, or you can monitor the truncation error due to the discretization, quitting only when it is tolerably small. In this section we will first discuss the “multigrid method,” then use the concepts developed to introduce the FMG method. The latter algorithm is the one that we implement in the accompanying programs. From One-Grid, through Two-Grid, to Multigrid The key idea of the multigrid method can be understood by considering the simplest case of a two-grid method. Suppose we are trying to solve the linear elliptic problem Lu = f (19.6.1) where L is some linear elliptic operator and f is the source term. Discretize equation (19.6.1) on a uniform grid with mesh size h. Write the resulting set of linear algebraic equations as Lhuh = fh (19.6.2) Let uh denote some approximate solution to equation (19.6.2). We will use the symbol uh to denote the exact solution to the difference equations (19.6.2). Then the error in uh or the correction is vh = uh − uh (19.6.3) The residual or defect is dh = Lhuh − fh (19.6.4) (Beware: some authors define residual as minus the defect, and there is not universal agreement about which of these two quantities 19.6.4 defines.) Since Lh is linear, the error satisfies Lhvh = −dh (19.6.5)
19.6 Multigrid Methods for Boundary Value Problems 873 At this point we need to make an approximation to Ch in order to find uh.The classical iteration methods,such as Jacobi or Gauss-Seidel,do this by finding,at each stage,an approximate solution of the equation Chin =-dh (19.6.6) where Ch is a"simpler"operator than Ch.For example,Ch is the diagonal part of Ch for Jacobi iteration,or the lower triangle for Gauss-Seidel iteration.The next approximation is generated by 好w=h+th (19.6.7) 菲 Now consider,as an alternative,a completely different type of approximation ICAL for Ch,one in which we“coarsify”rather than“simplify.”That is,.we form some 3 appropriate approximation Ci of Ch on a coarser grid with mesh size H(we will always take H=2h,but other choices are possible).The residual equation (19.6.5) RECIPES is now approximated by 令 CHUH =-dH (19.6.8) Since C has smaller dimension,this equation will be easier to solve than equation 9」 (19.6.5).To define the defect d on the coarse grid,we need a restriction operator R that restricts dh to the coarse grid: dH Rdh (19.6.9) The restriction operator is also called the fine-to-coarse operator or the injection operator.Once we have a solution n to equation (19.6.8),we need a prolongation operator P that prolongates or interpolates the correction to the fine grid: 三 Uh =PUH (19.6.10) Numerica 10621 431 The prolongation operator is also called the coarse-to-fine operator or the inter- Recipes polation operator.Both R and P are chosen to be linear operators.Finally the approximation un can be updated: 裙ew=h十th (19.6.11) One step of this coarse-grid correction scheme is thus: Coarse-Grid Correction Compute the defect on the fine grid from (19.6.4). .Restrict the defect by (19.6.9). .Solve (19.6.8)exactly on the coarse grid for the correction .Interpolate the correction to the fine grid by (19.6.10)
19.6 Multigrid Methods for Boundary Value Problems 873 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 machinereadable 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). At this point we need to make an approximation to Lh in order to find vh. The classical iteration methods, such as Jacobi or Gauss-Seidel, do this by finding, at each stage, an approximate solution of the equation Lhvh = −dh (19.6.6) where Lh is a “simpler” operator than Lh. For example, Lh is the diagonal part of Lh for Jacobi iteration, or the lower triangle for Gauss-Seidel iteration. The next approximation is generated by unew h = uh + vh (19.6.7) Now consider, as an alternative, a completely different type of approximation for Lh, one in which we “coarsify” rather than “simplify.” That is, we form some appropriate approximation LH of Lh on a coarser grid with mesh size H (we will always take H = 2h, but other choices are possible). The residual equation (19.6.5) is now approximated by LH vH = −dH (19.6.8) Since LH has smaller dimension, this equation will be easier to solve than equation (19.6.5). To define the defect dH on the coarse grid, we need a restriction operator R that restricts dh to the coarse grid: dH = Rdh (19.6.9) The restriction operator is also called the fine-to-coarse operator or the injection operator. Once we have a solution vH to equation (19.6.8), we need a prolongation operator P that prolongates or interpolates the correction to the fine grid: vh = PvH (19.6.10) The prolongation operator is also called the coarse-to-fine operator or the interpolation operator. Both R and P are chosen to be linear operators. Finally the approximation uh can be updated: unew h = uh + vh (19.6.11) One step of this coarse-grid correction scheme is thus: Coarse-Grid Correction • Compute the defect on the fine grid from (19.6.4). • Restrict the defect by (19.6.9). • Solve (19.6.8) exactly on the coarse grid for the correction. • Interpolate the correction to the fine grid by (19.6.10).
874 Chapter 19.Partial Differential Equations Compute the next approximation by (19.6.11). Let's contrast the advantages and disadvantages ofrelaxation and the coarse-grid correction scheme.Consider the error vh expanded into a discrete Fourier series.Call the components in the lower half of the frequency spectrum the smooth components and the high-frequency components the nonsmooth components.We have seen that relaxation becomes very slowly convergent in the limit h-0,i.e.,when there are a large number of mesh points.The reason turns out to be that the smooth components are only slightly reduced in amplitude on each iteration.However,many relaxation methods reduce the amplitude of the nonsmooth components by large factors on 三 each iteration:They are good smoothing operators. For the two-grid iteration,on the other hand,components of the error with g wavelengths S 2H are not even representable on the coarse grid and so cannot be reduced to zero on this grid.But it is exactly these high-frequency components that can be reduced by relaxation on the fine grid!This leads us to combine the ideas of relaxation and coarse-grid correction: Two-Grid Iteration 影 Pre-smoothing:Compute ih by applying v>0 steps of a relaxation method to uh. 9 .Coarse-grid correction:As above,using in to give inew Post-smoothing:Compute unew by applying v2>0 steps of the relaxation method to inew a以N 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 get an approximate solution of it by introducing an even coarser grid and using the two-grid OF SCIENTIFIC iteration method.If the convergence factor of the two-grid method is small enough, we will need only a few steps of this iteration to get a good enough approximate 6 solution.We denote the number of such iterations by y.Obviously we can apply this idea recursively down to some coarsest grid.There the solution is found easily,for example by direct matrix inversion or by iterating the relaxation scheme to convergence. One iteration of a multigrid method,from finest grid to coarser grids and back to finest grid again,is called a cycle.The exact structure of a cycle depends on the 10621 value ofy.the number of two-grid iterations at each intermediate stage.The case y=1 is called a V-cycle,while y=2 is called a W-cycle(see Figure 19.6.1).These Numerical Recipes 43106 are the most important cases in practice. Note that once more than two grids are involved,the pre-smoothing steps after (outside 腿 the first one on the finest grid need an initial approximation for the error v.This should be taken to be zero. North Smoothing,Restriction,and Prolongation Operators The most popular smoothing method,and the one you should try first,is Gauss-Seidel,since it usually leads to a good convergence rate.If we order the mesh points from I to N.then the Gauss-Seidel scheme is i=1,,N (19.6.12)
874 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 machinereadable 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). • Compute the next approximation by (19.6.11). Let’s contrast the advantages and disadvantages of relaxation and the coarse-grid correction scheme. Consider the error vh expanded into a discrete Fourier series. Call the components in the lower half of the frequency spectrum the smooth components and the high-frequency components the nonsmooth components. We have seen that relaxation becomes very slowly convergent in the limit h → 0, i.e., when there are a large number of mesh points. The reason turns out to be that the smooth components are only slightly reduced in amplitude on each iteration. However, many relaxation methods reduce the amplitude of the nonsmooth components by large factors on each iteration: They are good smoothing operators. For the two-grid iteration, on the other hand, components of the error with wavelengths <∼ 2H are not even representable on the coarse grid and so cannot be reduced to zero on this grid. But it is exactly these high-frequency components that can be reduced by relaxation on the fine grid! This leads us to combine the ideas of relaxation and coarse-grid correction: Two-Grid Iteration • Pre-smoothing: Compute u¯h by applying ν1 ≥ 0 steps of a relaxation method to uh. • Coarse-grid correction: As above, using u¯h to give u¯new h . • Post-smoothing: Compute unew h by applying ν2 ≥ 0 steps of the relaxation method to u¯new h . 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 get an approximate solution of it by introducing an even coarser grid and using the two-grid iteration method. If the convergence factor of the two-grid method is small enough, we will need only a few steps of this iteration to get a good enough approximate solution. We denote the number of such iterations by γ. Obviously we can apply this idea recursively down to some coarsest grid. There the solution is found easily, for example by direct matrix inversion or by iterating the relaxation scheme to convergence. One iteration of a multigrid method, from finest grid to coarser grids and back to finest grid again, is called a cycle. The exact structure of a cycle depends on the value of γ, the number of two-grid iterations at each intermediate stage. The case γ = 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 after the first one on the finest grid need an initial approximation for the error v. This should be taken to be zero. Smoothing, Restriction, and Prolongation Operators The most popular smoothing method, and the one you should try first, is Gauss-Seidel, since it usually leads to a good convergence rate. If we order the mesh points from 1 to N, then the Gauss-Seidel scheme is ui = − N j=1 j=i Lijuj − fi 1 Lii i = 1,...,N (19.6.12)
19.6 Multigrid Methods for Boundary Value Problems 875 ⑤ 2-grid E ⑤ ⑤ ⑨ 3-grid ⊙ © ⑨ E E E 83 旦 granted for from NUMERICAL RECIPES IN C: 19881992 ⑨ 1.200 4-grid ⑤ ⊙ server THE E E E E E (North America computer, ART Y=1 Y=2 Figure 19.6.1. Structure of multigrid cycles.S denotes smoothing,while E denotes exact solution 9 Progra on the coarsest grid.Each descending line denotes restriction (R)and each ascending line denotes prolongation (P).The finest grid is at the top level of each diagram.For the V-cycles (y 1)the E step is replaced by one 2-grid iteration each time the number of grid levels is increased by one.For the W-cycles(=2),each E step gets replaced by two 2-grid iterations. 受 to dir OF SCIENTIFIC COMPUTING (ISBN 1988-1992 where new values of u are used on the right-hand side as they become available.The exact form of the Gauss-Seidel method depends on the ordering chosen for the mesh 10-621 points.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-black Fuurggoglrion Numerical Recipes 43108 ordering,making one pass through the mesh updating the"even"points(like the red squares of a checkerboard)and another pass updating the "odd"points(the black squares).When quantities are more strongly coupled along one dimension than (outside another,one should relax a whole line along that dimension simultaneously.Line North Software. relaxation for nearest-neighbor coupling involves solving a tridiagonal system,and so is still efficient.Relaxing odd and even lines on successive passes is called zebra relaxation and is usually preferred over simple line relaxation. visit website Note that SOR should not be used as a smoothing operator.The overrelaxation machine destroys the high-frequency smoothing that is so crucial for the multigrid method. A succint notation for the prolongation and restriction operators is to give their symbol.The symbol of p is found by considering vr to be 1 at some mesh point (y),zero elsewhere,and then asking for the values of Pun.The most popular prolongation operator is simple bilinear interpolation.It gives nonzero values at the 9 points (,y),(+h,y),...,(-h,y-h),where the values are 1
19.6 Multigrid Methods for Boundary Value Problems 875 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 machinereadable 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). E γ = 1 γ = 2 2-grid 3-grid 4-grid S S S S S S E S S S S E S S E S S E S S S E S S E S S S S E S S S E S S Figure 19.6.1. Structure of multigrid cycles. S denotes smoothing, while E denotes exact solution on the coarsest grid. Each descending line \ denotes restriction (R) and each ascending line / denotes prolongation (P). The finest grid is at the top level of each diagram. For the V-cycles (γ = 1) the E step is replaced by one 2-grid iteration each time the number of grid levels is increased by one. For the W-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. The exact form of the Gauss-Seidel method depends on the ordering chosen for the mesh points. 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-black ordering, making one pass through the mesh updating the “even” points (like the red squares of a checkerboard) and another pass updating the “odd” points (the black squares). When quantities are more strongly coupled along one dimension than another, one should relax a whole line along that dimension simultaneously. Line relaxation for nearest-neighbor coupling involves solving a tridiagonal system, and so is still efficient. Relaxing odd and even lines on successive passes is called zebra relaxation and is usually preferred over simple line relaxation. Note that SOR should not be used as a smoothing operator. The overrelaxation destroys the high-frequency smoothing that is so crucial for the multigrid method. A succint notation for the prolongation and restriction operators is to give their symbol. The symbol of P is found by considering v H to be 1 at some mesh point (x, y), zero elsewhere, and then asking for the values of Pv H. The most popular prolongation operator is simple bilinear interpolation. It gives nonzero values at the 9 points (x, y), (x + h, y),...,(x − h, y − h), where the values are 1, 1 2 ,..., 1 4
876 Chapter 19.Partial Differential Equations Its symbol is therefore 12 12 (19.6.13) 12 The symbol of R is defined by considering vn to be defined everywhere on the fine grid,and then asking what is Ruh at (y)as a linear combination of these values.The simplest possible choice for R is straight injection,which means simply filling each coarse-grid point with the value from the corresponding fine-grid point. 81 Its symbol is "[1]."However,difficulties can arise in practice with this choice.It turns out that a safe choice for R is to make it the adjoint operator to P.To define the adjoint,define the scalar product of two grid functions uh and vh for mesh size h as (uhnh三h2∑h(c,y)h(,) (19.6.14) 3 工，y RECIPES Then the adjoint of P,denoted Pt,is defined by ≥ (uHIPivh)H =(PuHlUh)h (19.6.15) Now take P to be bilinear interpolation,and choose u=1 at(,y),zero elsewhere. Set Pt R in (19.6.15)and H 2h.You will find that e9兰公 (Rh)红，)=h(z,)+h(z+h,)+h(z+h,y+h)+…(19.6.16) IENTIFIC so that the symbol of R is 1 18 (19.6.17) 」 Note the simple rule:The symbol of R is the transpose of the matrix defining the symbol ofP,equation(19.6.13).This rule is general wheneverR=Pt and H=2h. Recipes Numerical 10521 The particular choice of R in(19.6.17)is called full weighting.Another popular 431 choice for R is half weighting,"halfway"between full weighting and straight Recipes injection.Its symbol is (outside 0 01 North 18 18 (19.6.18) 0 0 A similar notation can be used to describe the difference operator Ch For example,the standard differencing of the model problem,equation (19.0.6),is represented by the five-point difference star 01 0 Ch= h2 -4 (19.6.19) 0 10
876 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 machinereadable 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). Its symbol is therefore 1 4 1 2 1 4 1 2 1 1 2 1 4 1 2 1 4 (19.6.13) The symbol of R is defined by considering vh to be defined everywhere on the fine grid, and then asking what is Rvh at (x, y) as a linear combination of these values. The simplest possible choice for R is straight injection, which means simply filling 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. It turns out that a safe choice for R is to make it the adjoint operator to P. To define the adjoint, define the scalar product of two grid functions u h and vh for mesh size h as uh|vh h ≡ h2 x,y uh(x, y)vh(x, y) (19.6.14) Then the adjoint of P, denoted P†, is defined by uH|P† vh H = PuH|vh h (19.6.15) Now take P to be bilinear interpolation, and choose u H = 1 at(x, y), zero elsewhere. Set P† = R in (19.6.15) and H = 2h. You will find that (Rvh)(x,y) = 1 4 vh(x, y) + 1 8 vh(x + h, y) + 1 16 vh(x + h, y + h) + ··· (19.6.16) so that the symbol of R is 1 16 1 8 1 16 1 8 1 4 1 8 1 16 1 8 1 16 (19.6.17) Note the simple rule: The symbol of R is 1 4 the transpose of the matrix defining the symbol ofP, 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 popular choice for R is half weighting, “halfway” between full weighting and straight injection. Its symbol is 0 1 8 0 1 8 1 2 1 8 0 1 8 0 (19.6.18) A similar notation can be used to describe the difference operator Lh. For example, the standard differencing of the model problem, equation (19.0.6), is represented by the five-point difference star Lh = 1 h2 0 10 1 −4 1 0 10 (19.6.19)
19.6 Multigrid Methods for Boundary Value Problems 877 If you are confronted with a new problem and you are not sure what P and R choices are likely to work well,here is a safe rule:Suppose mp is the order of the interpolation P(ie.,it interpolates polynomials of degree mp-1 exactly).Suppose mr is the order of R,and that R is the adjoint of some P(not necessarily the P you intend to use).Then if m is the order of the differential operator Ch,you should satisfy the inequality mp+m>m.For example,bilinear interpolation and its adjoint,full weighting,for Poisson's equation satisfy mp+mr=4>m =2. Of course the p and R operators should enforce the boundary conditions for your problem.The easiest way to do this is to rewrite the difference equation to have homogeneous boundary conditions by modifying the source term if necessary 81 (cf.$19.4).Enforcing homogeneous boundary conditions simply requires the P operator to produce zeros at the appropriate boundary points.The corresponding R is then found by R=Pt Full Multigrid Algorithm So far we have described multigrid as an iterative scheme,where one starts with some initial guess on the finest grid and carries out enough cycles(V-cycles, 9 W-cycles,...)to achieve convergence.This is the simplest way to use multigrid: Simply apply enough cycles until some appropriate convergence criterion is met However,efficiency can be improved by using the Full Multigrid Algorithm(FMG), also known as nested iteration. Instead of starting with an arbitrary approximation on the finest grid (e.g., uh=0),the first approximation is obtained by interpolating from a coarse-grid 、屋多5 solution: uh =PuH (19.6.20) 61 The coarse-grid solution itself is found by a similar FMG process from even coarser grids.At the coarsest level,you start with the exact solution.Rather than proceed as in Figure 19.6.1,then,FMG gets to its solution by a series of increasingly tall "N's," each taller one probing a finer grid(see Figure 19.6.2). Numerica 10621 Note that P in(19.6.20)need not be the same P used in the multigrid cycles. It should be at least of the same order as the discretization Ch,but sometimes a 431 higher-order operator leads to greater efficiency. Recipes It turns out that you usually need one or at most two multigrid cycles at each level before proceeding down to the next finer grid.While there is theoretical (outside 腿 guidance on the required number of cycles(e.g.,[21),you can easily determine it North empirically.Fix the finest level and study the solution values as you increase the number of cycles per level.The asymptotic value of the solution is the exact solution of the difference equations.The difference between this exact solution and the solution for a small number of cycles is the iteration error.Now fix the number of cycles to be large,and vary the number oflevels,i.e.,the smallest value of h used.In this way you can estimate the truncation error for a given h.In your final production code,there is no point in using more cycles than you need to get the iteration error down to the size of the truncation error. The simple multigrid iteration (cycle)needs the right-hand side f only at the finest level.FMG needs f at all levels.If the boundary conditions are homogeneous
19.6 Multigrid Methods for Boundary Value Problems 877 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 machinereadable 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). If you are confronted with a new problem and you are not sure what P and R choices are likely to work well, here is a safe rule: Suppose mp is the order of the interpolation P (i.e., it interpolates polynomials of degree mp − 1 exactly). Suppose mr is the order of R, and that R is the adjoint of some P (not necessarily the P you intend to use). Then if m is the order of the differential operator Lh, you should satisfy the inequality mp + mr > m. For example, bilinear interpolation and its adjoint, full weighting, for Poisson’s equation satisfy mp + mr = 4 > m = 2. Of course the P and R operators should enforce the boundary conditions for your problem. The easiest way to do this is to rewrite the difference equation to have homogeneous boundary conditions by modifying the source term if necessary (cf. §19.4). Enforcing homogeneous boundary conditions simply requires the P operator to produce zeros at the appropriate boundary points. The corresponding R is then found by R = P†. Full Multigrid Algorithm So far we have described multigrid as an iterative scheme, where one starts with some initial guess on the finest grid and carries out enough cycles (V-cycles, W-cycles,...) to achieve convergence. This is the simplest way to use multigrid: Simply apply enough cycles until some appropriate convergence criterion is met. However, efficiency can be improved by using the Full Multigrid Algorithm (FMG), also known as nested iteration. Instead of starting with an arbitrary approximation on the finest grid (e.g., uh = 0), the first approximation is obtained by interpolating from a coarse-grid solution: uh = PuH (19.6.20) The coarse-grid solution itself is found by a similar FMG process from even coarser grids. At the coarsest level, you start with the exact solution. Rather than proceed as in Figure 19.6.1, then, FMG gets to its solution by a series of increasingly tall “N’s,” each taller one probing a finer grid (see Figure 19.6.2). Note that P in (19.6.20) need not be the same P used in the multigrid cycles. It should be at least of the same order as the discretization Lh, but sometimes a higher-order operator leads to greater efficiency. It turns out that you usually need one or at most two multigrid cycles at each level before proceeding down to the next finer grid. While there is theoretical guidance on the required number of cycles (e.g., [2]), you can easily determine it empirically. Fix the finest level and study the solution values as you increase the number of cycles per level. The asymptotic value of the solution is the exact solution of the difference equations. The difference between this exact solution and the solution for a small number of cycles is the iteration error. Now fix the number of cycles to be large, and vary the number of levels, i.e., the smallest value of h used. In this way you can estimate the truncation error for a given h. In your final production code, there is no point in using more cycles than you need to get the iteration error down to the size of the truncation error. The simple multigrid iteration (cycle) needs the right-hand side f only at the finest level. FMG needs f at all levels. If the boundary conditions are homogeneous
878 Chapter 19. Partial Differential Equations ③ ⑤ ⑨ © © 4-grid ⑨⑤ ncycle =1 自间 E E http://www.nr read able files ⑨ © ⑤ 83 granted for ③ ⑨ ⑤ © §③ ⑤ 1100 9 ⊙ © 4-grid 872 ncycle =2 自自 自 自 E E (Nort server Figure 19.6.2.Structure of cycles for the full multigrid(FMG)method.This method starts on the coarsest grid,interpolates,and then refines (by "V's"),the solution onto grids of increasing fineness. Americ University Press. make one paper computer, THE ART you can use fr=Rfh.This prescription is not always safe for inhomogeneous boundary conditions.In that case it is better to discretize f on each coarse grid. Progra Note that the FMG algorithm produces the solution on all levels.It can therefore be combined with techniques like Richardson extrapolation. We now give a routine mglin that implements the Full Multigrid Algorithm 6 %州 for a linear equation,the model problem(19.0.6).It uses red-black Gauss-Seidel as the smoothing operator,bilinear interpolation for P,and half-weighting for R.To change the routine to handle another linear problem,all you need do is modify the OF SCIENTIFIC COMPUTING (ISBN functions relax,resid,and slvsml appropriately.A feature of the routine is the dynamical allocation of storage for variables defined on the various grids. 、 1988-18920 10-621 #include "nrutil.h" #define NPRE 1 Number of relaxation sweeps before... Numerical Recipes #define NPOST 1 ..and after the coarse-grid correction is com- #define NGMAX 15 puted. (outside 膜 void mglin(double **u,int n,int ncycle) oftware. Full Multigrid Algorithm for solution of linear elliptic equation,here the model problem(19.0.6) On input u[1..n][1..n]contains the right-hand side p,while on output it returns the solution. The dimension n must be of the form 2+1 for some integer j.(j is actually the number of grid levels used in the solution,called ng below.)ncycle is the number of V-cycles to be B used at each level. void addint(double **uf,double **uc,double **res,int nf); void copy(double **aout,double **ain,int n); void fill0(double **u,int n); void interp(double **uf,double **uc,int nf); void relax(double **u,double +*rhs,int n); void resid(double **res,double *u,double **rhs,int n); void rstrct(double **uc,double **uf,int nc); void slvsml(double **u,double **rhs);
878 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 machinereadable 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). 4-grid ncycle = 1 4-grid ncycle = 2 S S S S S S S S S S S E E S S S S S E E E S S S S S E E S S S E S S S S S E E S S S S E S S S S S E S S S E E S S S S S S S S E Figure 19.6.2. Structure of cycles for the full multigrid (FMG) method. This method starts on the coarsest grid, interpolates, and then refines (by “V’s”), the solution onto grids of increasing fineness. you can use fH = Rfh. This prescription is not always safe for inhomogeneous boundary conditions. In that case it is better to discretize f on each coarse grid. Note that the FMG algorithm produces the solution on all levels. It can therefore be combined with techniques like Richardson extrapolation. We now give a routine mglin that implements the Full Multigrid Algorithm for a linear equation, the model problem (19.0.6). It uses red-black Gauss-Seidel as the smoothing operator, bilinear interpolation for P, and half-weighting for R. To change the routine to handle another linear problem, all you need do is modify the functions relax, resid, and slvsml appropriately. A feature of the routine is the dynamical allocation of storage for variables defined on the various grids. #include "nrutil.h" #define NPRE 1 Number of relaxation sweeps before ... #define NPOST 1 ... and after the coarse-grid correction is com- #define NGMAX 15 puted. void mglin(double **u, int n, int ncycle) Full Multigrid Algorithm for solution of linear elliptic equation, here the model problem (19.0.6). On input u[1..n][1..n] contains the right-hand side ρ, while on output it returns the solution. The dimension n must be of the form 2j + 1 for some integer j. (j is actually the number of grid levels used in the solution, called ng below.) ncycle is the number of V-cycles to be used at each level. { void addint(double **uf, double **uc, double **res, int nf); void copy(double **aout, double **ain, int n); void fill0(double **u, int n); void interp(double **uf, double **uc, int nf); void relax(double **u, double **rhs, int n); void resid(double **res, double **u, double **rhs, int n); void rstrct(double **uc, double **uf, int nc); void slvsml(double **u, double **rhs);
19.6 Multigrid Methods for Boundary Value Problems 879 unsigned int j,jcycle,jj,jpost,jpre,nf,ng=0,ngrid,nn; double **ires [NGMAX+1],**irho [NGMAX+1],**irhs [NGMAX+1],**iu[NGMAX+1]; nn=n; while (nn >>1)ng++; if (n !1+(1L =2;jj--){ Downward stoke of the V. for (ipre=1;ipre=2;j--,nn=nn/2+1){ free_dmatrix(ires[i],1,nn,1,nn); free_dmatrix(irhs[j],1,nn,1,nn); free_dmatrix(iu[j],1,nn,1,nn); if (j !ng)free_dmatrix(irho[j],1,nn,1,nn); 2 free_dmatrix(irhs[1],1,3,1,3); free_dmatrix(iu[1],1,3,1,3);
19.6 Multigrid Methods for Boundary Value Problems 879 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 machinereadable 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). unsigned int j,jcycle,jj,jpost,jpre,nf,ng=0,ngrid,nn; double **ires[NGMAX+1],**irho[NGMAX+1],**irhs[NGMAX+1],**iu[NGMAX+1]; nn=n; while (nn >>= 1) ng++; if (n != 1+(1L NGMAX) nrerror("increase NGMAX in mglin."); nn=n/2+1; ngrid=ng-1; irho[ngrid]=dmatrix(1,nn,1,nn); Allocate storage for r.h.s. on grid ng − 1, rstrct(irho[ngrid],u,nn); and fill it by restricting from the fine grid. while (nn > 3) { Similarly allocate storage and fill r.h.s. on all nn=nn/2+1; coarse grids. irho[--ngrid]=dmatrix(1,nn,1,nn); rstrct(irho[ngrid],irho[ngrid+1],nn); } nn=3; iu[1]=dmatrix(1,nn,1,nn); irhs[1]=dmatrix(1,nn,1,nn); slvsml(iu[1],irho[1]); Initial solution on coarsest grid. free_dmatrix(irho[1],1,nn,1,nn); ngrid=ng; for (j=2;j=2;jj--) { Downward stoke of the V. for (jpre=1;jpre=2;j--,nn=nn/2+1) { free_dmatrix(ires[j],1,nn,1,nn); free_dmatrix(irhs[j],1,nn,1,nn); free_dmatrix(iu[j],1,nn,1,nn); if (j != ng) free_dmatrix(irho[j],1,nn,1,nn); } free_dmatrix(irhs[1],1,3,1,3); free_dmatrix(iu[1],1,3,1,3); }
880 Chapter 19.Partial Differential Equations void rstrct(double **uc,double **uf,int nc) Half-weighting restriction.nc is the coarse-grid dimension.The fine-grid solution is input in uf [1..2*nc-1][1..2*nc-1],the coarse-grid solution is returned in uc [1..nc][1..nc]. int ic,iif,jc,jf,ncc=2*nc-1; for (jf=3,jc=2;jc<nc;jc++,jf+=2) Interior points for(1if=3,1c=2;ic<nc;1c++,11f+=2)[ uc[1c][jc]=0.5*uf[iif][f]+0.125*(uf[i1f+1][jf]+uf[1if-1][jf] +uf [iif][if+1]+uf[iif][if-1]); for(jc=1,1c=1;1c<=nc;1c++,jc+=2)[ Boundary points. uc[ic][1]=uf[jc][1]; uc [ic][nc]=uf [jc][ncc]: granted for 19881992 for(jc=1,1c=1;ic<=nc;1c++,jc+=2){ uc[1][ic]=uf[1][jc]; 1600 9 uc [ne][ic]=uf [ncc][jc]; 2 NUMERICAL RECIPES I void interp(double **uf,double **uc,int nf) (Nort server Coarse-to-fine prolongation by bilinear interpolation.nf is the fine-grid dimension.The coarse- to make grid solution is input as uc [1..nc][1..nc],where nc nf/2+1.The fine-grid solution is returned in uf [1..nf][1..nf]. int ic,iif,jc,jf,nci nc=nf/2+1; for (jc=1,jf=1;jc<=nc;jc++,jf+=2) Do elements that are copies. 9 for (ic=1;ic<=nc;ic++)uf[2*ic-1][jf]=uc[ic][jc]; for (jf=1;jf<=nf;jf+=2) Do odd-numbered columns,interpolat- for (iif=2;iif<nf;iif+=2) ing vertically. IENTIFIC uf[iif][if]=0.5*(uf [iif+1][if]+uf[iif-1][if]); for (jf=2;jf<nf;jf+=2) Do even-numbered columns,interpolat- for (iif=1;iif <nf;iif++) ing horizontally. uf[iif][if]=0.5*(uf [iif][if+1]+uf [iif][if-1]); 1992 COMPUTING(ISBN void addint(double **uf,double **uc,double **res,int nf) Does coarse-to-fine interpolation and adds result to uf.nf is the fine-grid dimension.The 021 coarse-grid solution is input as uc [1..nc][1..nc],where nc nf/2+1.The fine-grid solu- tion is returned in uf [1..nf][1..nf].res [1..nf][1..nf]is used for temporary storage. Numerical Recipes 43108 void interp(double **uf,double **uc,int nf); int i,j; (outside interp(res,uc,nf); for (j=1;j<=nf;j++) North Software. for (i=1;i<=nf;i++) uf [i][j]+res[i][j]; void slvsml(double *u,double **rhs) Solution of the model problem on the coarsest grid,where h=.The right-hand side is input in rhs [1..3][1..3]and the solution is returned in u[1..3][1..3]. void fill0(double **u,int n); double h=0.5; f1110(u,3);
880 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 machinereadable 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). void rstrct(double **uc, double **uf, int nc) Half-weighting restriction. nc is the coarse-grid dimension. The fine-grid solution is input in uf[1..2*nc-1][1..2*nc-1], the coarse-grid solution is returned in uc[1..nc][1..nc]. { int ic,iif,jc,jf,ncc=2*nc-1; for (jf=3,jc=2;jc<nc;jc++,jf+=2) { Interior points. for (iif=3,ic=2;ic<nc;ic++,iif+=2) { uc[ic][jc]=0.5*uf[iif][jf]+0.125*(uf[iif+1][jf]+uf[iif-1][jf] +uf[iif][jf+1]+uf[iif][jf-1]); } } for (jc=1,ic=1;ic<=nc;ic++,jc+=2) { Boundary points. uc[ic][1]=uf[jc][1]; uc[ic][nc]=uf[jc][ncc]; } for (jc=1,ic=1;ic<=nc;ic++,jc+=2) { uc[1][ic]=uf[1][jc]; uc[nc][ic]=uf[ncc][jc]; } } void interp(double **uf, double **uc, int nf) Coarse-to-fine prolongation by bilinear interpolation. nf is the fine-grid dimension. The coarsegrid solution is input as uc[1..nc][1..nc], where nc = nf/2+1. The fine-grid solution is returned in uf[1..nf][1..nf]. { int ic,iif,jc,jf,nc; nc=nf/2+1; for (jc=1,jf=1;jc<=nc;jc++,jf+=2) Do elements that are copies. for (ic=1;ic<=nc;ic++) uf[2*ic-1][jf]=uc[ic][jc]; for (jf=1;jf<=nf;jf+=2) Do odd-numbered columns, interpolatfor (iif=2;iif<nf;iif+=2) ing vertically. uf[iif][jf]=0.5*(uf[iif+1][jf]+uf[iif-1][jf]); for (jf=2;jf<nf;jf+=2) Do even-numbered columns, interpolatfor (iif=1;iif <= nf;iif++) ing horizontally. uf[iif][jf]=0.5*(uf[iif][jf+1]+uf[iif][jf-1]); } void addint(double **uf, double **uc, double **res, int nf) Does coarse-to-fine interpolation and adds result to uf. nf is the fine-grid dimension. The coarse-grid solution is input as uc[1..nc][1..nc], where nc = nf/2+1. The fine-grid solution is returned in uf[1..nf][1..nf]. res[1..nf][1..nf] is used for temporary storage. { void interp(double **uf, double **uc, int nf); int i,j; interp(res,uc,nf); for (j=1;j<=nf;j++) for (i=1;i<=nf;i++) uf[i][j] += res[i][j]; } void slvsml(double **u, double **rhs) Solution of the model problem on the coarsest grid, where h = 1 2 . The right-hand side is input in rhs[1..3][1..3] and the solution is returned in u[1..3][1..3]. { void fill0(double **u, int n); double h=0.5; fill0(u,3);