c19-6 (779621), страница 3
Текст из файла (страница 3)
. . ) 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-gridsolution:uh = PuH(19.6.20)The coarse-grid solution itself is found by a similar FMG process from even coarsergrids.
At the coarsest level, you start with the exact solution. Rather than proceed asin 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 ahigher-order operator leads to greater efficiency.It turns out that you usually need one or at most two multigrid cycles at eachlevel before proceeding down to the next finer grid.
While there is theoreticalguidance on the required number of cycles (e.g., [2]), you can easily determine itempirically. Fix the finest level and study the solution values as you increase thenumber of cycles per level. The asymptotic value of the solution is the exact solutionof the difference equations. The difference between this exact solution and thesolution for a small number of cycles is the iteration error. Now fix the number ofcycles to be large, and vary the number of levels, i.e., the smallest value of h used.
Inthis way you can estimate the truncation error for a given h. In your final productioncode, there is no point in using more cycles than you need to get the iteration errordown to the size of the truncation error.The simple multigrid iteration (cycle) needs the right-hand side f only at thefinest level. FMG needs f at all levels. If the boundary conditions are homogeneous,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).If you are confronted with a new problem and you are not sure what P and Rchoices are likely to work well, here is a safe rule: Suppose mp is the order of theinterpolation P (i.e., it interpolates polynomials of degree mp − 1 exactly).
Supposemr is the order of R, and that R is the adjoint of some P (not necessarily the P youintend to use). Then if m is the order of the differential operator Lh , you shouldsatisfy the inequality mp + mr > m. For example, bilinear interpolation and itsadjoint, full weighting, for Poisson’s equation satisfy mp + mr = 4 > m = 2.Of course the P and R operators should enforce the boundary conditions foryour problem.
The easiest way to do this is to rewrite the difference equation tohave homogeneous boundary conditions by modifying the source term if necessary(cf. §19.4). Enforcing homogeneous boundary conditions simply requires the Poperator to produce zeros at the appropriate boundary points. The correspondingR is then found by R = P † .878Chapter 19.Partial Differential EquationsSSSSSSESSSESESSSESESESSSESSSE4-gridncycle = 1SSSSSESSSSS4-gridncycle = 2EFigure 19.6.2. Structure of cycles for the full multigrid (FMG) method. This method starts on thecoarsest 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 inhomogeneousboundary 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 thereforebe combined with techniques like Richardson extrapolation.We now give a routine mglin that implements the Full Multigrid Algorithmfor a linear equation, the model problem (19.0.6).
It uses red-black Gauss-Seidel asthe smoothing operator, bilinear interpolation for P, and half-weighting for R. Tochange the routine to handle another linear problem, all you need do is modify thefunctions relax, resid, and slvsml appropriately. A feature of the routine is thedynamical allocation of storage for variables defined on the various grids.#include "nrutil.h"#define NPRE 1#define NPOST 1#define NGMAX 15Number of relaxation sweeps before .
. .. . . and after the coarse-grid correction is computed.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 ofgrid levels used in the solution, called ng below.) ncycle is the number of V-cycles to beused 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);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).ESS19.6 Multigrid Methods for Boundary Value Problems879unsigned int j,jcycle,jj,jpost,jpre,nf,ng=0,ngrid,nn;double **ires[NGMAX+1],**irho[NGMAX+1],**irhs[NGMAX+1],**iu[NGMAX+1];}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).nn=n;while (nn >>= 1) ng++;if (n != 1+(1L << ng)) nrerror("n-1 must be a power of 2 in mglin.");if (ng > 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 allnn=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<=ngrid;j++) {Nested iteration loop.nn=2*nn-1;iu[j]=dmatrix(1,nn,1,nn);irhs[j]=dmatrix(1,nn,1,nn);ires[j]=dmatrix(1,nn,1,nn);interp(iu[j],iu[j-1],nn);Interpolate from coarse grid to next finer grid.copy(irhs[j],(j != ngrid ? irho[j] : u),nn);Set up r.h.s.for (jcycle=1;jcycle<=ncycle;jcycle++) {V-cycle loop.nf=nn;for (jj=j;jj>=2;jj--) {Downward stoke of the V.for (jpre=1;jpre<=NPRE;jpre++)Pre-smoothing.relax(iu[jj],irhs[jj],nf);resid(ires[jj],iu[jj],irhs[jj],nf);nf=nf/2+1;rstrct(irhs[jj-1],ires[jj],nf);Restriction of the residual is the next r.h.s.fill0(iu[jj-1],nf);Zero for initial guess in next}relaxation.slvsml(iu[1],irhs[1]);Bottom of V: solve on coarsnf=3;est grid.for (jj=2;jj<=j;jj++) {Upward stroke of V.nf=2*nf-1;addint(iu[jj],iu[jj-1],ires[jj],nf);Use res for temporary storage inside addint.for (jpost=1;jpost<=NPOST;jpost++)Post-smoothing.relax(iu[jj],irhs[jj],nf);}}}copy(u,iu[ngrid],n);Return solution in u.for (nn=n,j=ng;j>=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);880Chapter 19.Partial Differential Equationsvoid rstrct(double **uc, double **uf, int nc)Half-weighting restriction.
nc is the coarse-grid dimension. The fine-grid solution is input inuf[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;}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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
