c19-6 (779621), страница 5
Текст из файла (страница 5)
If we rewrite this asLh (u) = fh + τ(19.6.31)we see that τ can be regarded as the correction to fh so that the solution of the fine-gridequation will be the exact solution u.Now consider the relative truncation error τh , which is defined on the H-grid relativeto the h-grid:τh ≡ LH (Ruh) − RLh (uh )(19.6.32)Since Lh (uh ) = fh , this can be rewritten asLH (uH ) = fH + τh(19.6.33)In other words, we can think of τh as the correction to fH that makes the solution of thecoarse-grid equation equal to the fine-grid solution. Of course we cannot compute τh , but wedo have an approximation to it from using ueh in equation (19.6.32):τh ' τeh ≡ LH (Reuh ) − RLh (euh )(19.6.34)Replacing τh by τeh in equation (19.6.33) givesLH (uH ) = LH (Reuh ) − Rdh(19.6.35)which is just the coarse-grid equation (19.6.27)!Thus we see that there are two complementary viewpoints for the relation betweencoarse and fine grids:• Coarse grids are used to accelerate the convergence of the smooth componentsof the fine-grid residuals.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).on the coarse grid. (This is how nonzero right-hand sides appear.) Suppose the approximatesolution is ueH .
Then the coarse-grid correction is884Chapter 19.Partial Differential Equations• Fine grids are used to compute correction terms to the coarse-grid equations,yielding fine-grid accuracy on the coarse grids.One benefit of this new viewpoint is that it allows us to derive a natural stopping criterionfor a multigrid iteration. Normally the criterion would bekdh k ≤ (19.6.36)τ = Lh (u) − Lh (uh ) = h2 τ2 (x, y) + · · ·(19.6.37)Assume the solution satisfies uh = u + h2 u2 (x, y) + · · · .
Then, assuming R is of highenough order that we can neglect its effect, equation (19.6.32) givesτh ' LH (u + h2 u2 ) − Lh (u + h2 u2 )= LH (u) − Lh (u) + h2 [L0H (u2 ) − L0h (u2 )] + · · ·(19.6.38)= (H 2 − h2 )τ2 + O(h4 )For the usual case of H = 2h we therefore haveτ ' 13 τh ' 13 τeh(19.6.39)The stopping criterion is thus equation (19.6.36) with = αkeτh k,α∼13(19.6.40)We have one remaining task before implementing our nonlinear multigrid algorithm:choosing a nonlinear relaxation scheme. Once again, your first choice should probably bethe nonlinear Gauss-Seidel scheme. If the discretized equation (19.6.23) is written withsome choice of ordering asLi (u1 , .
. . , uN ) = fi ,i = 1, . . . , N(19.6.41)then the nonlinear Gauss-Seidel schemes solvesLi (u1 , . . . , ui−1 , unew, ui+1 , . . . , uN ) = fii(19.6.42)for unew. As usual new u’s replace old u’s as soon as they have been computed. Often equationi(19.6.42) is linear in unew, since the nonlinear terms are discretized by means of its neighbors.iIf this is not the case, we replace equation (19.6.42) by one step of a Newton iteration:Li (uoldi ) − fi∂Li (uoldi )/∂uiFor example, consider the simple nonlinear equationunew= uold−ii∇2 u + u2 = ρ(19.6.43)(19.6.44)In two-dimensional notation, we haveL(ui,j ) = (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j )/h2 + u2i,j − ρi,j = 0 (19.6.45)Since∂L= −4/h2 + 2ui,j(19.6.46)∂ui,jthe Newton Gauss-Seidel iteration isL(ui,j )unew(19.6.47)i,j = ui,j −−4/h2 + 2ui,jHere is a routine mgfas that solves equation (19.6.44) using the Full Multigrid Algorithmand the FAS scheme.
Restriction and prolongation are done as in mglin. We have includedthe convergence test based on equation (19.6.40). A successful multigrid solution of a problemshould aim to satisfy this condition with the maximum number of V-cycles, maxcyc, equal to1 or 2.
The routine mgfas uses the same functions copy, interp, and rstrct as mglin.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).and the question is how to choose . There is clearly no benefit in iterating beyond thepoint when the remaining error is dominated by the local truncation error τ . The computablequantity is τeh . What is the relation between τ and τeh ? For the typical case of a second-orderaccurate differencing scheme,19.6 Multigrid Methods for Boundary Value Problems#include "nrutil.h"#define NPRE 1#define NPOST 1#define ALPHA 0.33#define NGMAX 15885Number of relaxation sweeps before . .
.. . . and after the coarse-grid correction is computed.Relates the estimated truncation error to the normof the residual.nn=n;while (nn >>= 1) ng++;if (n != 1+(1L << ng)) nrerror("n-1 must be a power of 2 in mgfas.");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);itau[1]=dmatrix(1,nn,1,nn);itemp[1]=dmatrix(1,nn,1,nn);slvsm2(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);itau[j]=dmatrix(1,nn,1,nn);itemp[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<=maxcyc;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.relax2(iu[jj],irhs[jj],nf);lop(itemp[jj],iu[jj],nf);Lh (euh ).nf=nf/2+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).void mgfas(double **u, int n, int maxcyc)Full Multigrid Algorithm for FAS solution of nonlinear elliptic equation, here equation (19.6.44).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.) maxcyc is the maximum number of V-cyclesto be used at each level.{double anorm2(double **a, int n);void copy(double **aout, double **ain, int n);void interp(double **uf, double **uc, int nf);void lop(double **out, double **u, int n);void matadd(double **a, double **b, double **c, int n);void matsub(double **a, double **b, double **c, int n);void relax2(double **u, double **rhs, int n);void rstrct(double **uc, double **uf, int nc);void slvsm2(double **u, double **rhs);unsigned int j,jcycle,jj,jm1,jpost,jpre,nf,ng=0,ngrid,nn;double **irho[NGMAX+1],**irhs[NGMAX+1],**itau[NGMAX+1],**itemp[NGMAX+1],**iu[NGMAX+1];double res,trerr;886Chapter 19.Partial Differential Equations}slvsm2(iu[1],irhs[1]);nf=3;for (jj=2;jj<=j;jj++) {jm1=jj-1;rstrct(itemp[jm1],iu[jj],nf);matsub(iu[jm1],itemp[jm1],itemp[jm1],nf);nf=2*nf-1;interp(itau[jj],itemp[jm1],nf);matadd(iu[jj],itau[jj],iu[jj],nf);for (jpost=1;jpost<=NPOST;jpost++)relax2(iu[jj],irhs[jj],nf);}lop(itemp[j],iu[j],nf);matsub(itemp[j],irhs[j],itemp[j],nf);res=anorm2(itemp[j],nf);if (res < trerr) break;}Bottom of V: Solve on coarsest grid.Upward stroke of V.Reuh .ueH − Reuh .P (euH −Reuh ) stored in τeh .Form uenewh .Post-smoothing.Form residual kdh k.No more V-cycles needed ifresidual small enough.}copy(u,iu[ngrid],n);Return solution in u.for (nn=n,j=ng;j>=1;j--,nn=nn/2+1) {free_dmatrix(itemp[j],1,nn,1,nn);free_dmatrix(itau[j],1,nn,1,nn);free_dmatrix(irhs[j],1,nn,1,nn);free_dmatrix(iu[j],1,nn,1,nn);if (j != ng && j != 1) free_dmatrix(irho[j],1,nn,1,nn);}}void relax2(double **u, double **rhs, int n)Red-black Gauss-Seidel relaxation for equation (19.6.44).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
