c19-6 (779621), страница 4
Текст из файла (страница 4)
The fine-grid solution isreturned 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. Thecoarse-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 = 12 .
The right-hand side is inputin 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);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).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];}19.6 Multigrid Methods for Boundary Value Problems881u[2][2] = -h*h*rhs[2][2]/4.0;}h=1.0/(n-1);h2=h*h;for (ipass=1;ipass<=2;ipass++,jsw=3-jsw) {Red and black sweeps.isw=jsw;for (j=2;j<n;j++,isw=3-isw)for (i=isw+1;i<n;i+=2)Gauss-Seidel formula.u[i][j]=0.25*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]-h2*rhs[i][j]);}}void resid(double **res, double **u, double **rhs, int n)Returns minus the residual for the model problem.
Input quantities are u[1..n][1..n] andrhs[1..n][1..n], while res[1..n][1..n] is returned.{int i,j;double h,h2i;h=1.0/(n-1);h2i=1.0/(h*h);for (j=2;j<n;j++)Interior points.for (i=2;i<n;i++)res[i][j] = -h2i*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]4.0*u[i][j])+rhs[i][j];for (i=1;i<=n;i++)Boundary points.res[i][1]=res[i][n]=res[1][i]=res[n][i]=0.0;}void copy(double **aout, double **ain, int n)Copies ain[1..n][1..n] to aout[1..n][1..n].{int i,j;for (i=1;i<=n;i++)for (j=1;j<=n;j++)aout[j][i]=ain[j][i];}void fill0(double **u, int n)Fills u[1..n][1..n] with zeros.{int i,j;for (j=1;j<=n;j++)for (i=1;i<=n;i++)u[i][j]=0.0;}The routine mglin is written for clarity, not maximum efficiency, so that it iseasy to modify. Several simple changes will speed up the execution time: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 relax(double **u, double **rhs, int n)Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solutionu[1..n][1..n], using the right-hand side function rhs[1..n][1..n].{int i,ipass,isw,j,jsw=1;double h,h2;882Chapter 19.Partial Differential EquationsNonlinear Multigrid: The FAS AlgorithmNow turn to solving a nonlinear elliptic equation, which we write symbolically asL(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 anonzero 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, whichproduces linear equations for the correction term at each iteration. We can then use linearmultigrid to solve these equations.
A great strength of the multigrid idea, however, is that itcan be applied directly to nonlinear problems. All we need is a suitable nonlinear relaxationmethod 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 nonlinearequations need be solved, except perhaps on the coarsest grid.To develop the nonlinear algorithm, suppose we have a relaxation procedure that cansmooth the residual vector as we did in the linear case.
Then we can seek a smooth correctionvh to solve (19.6.23):Lh (euh + vh ) = fh(19.6.24)To find vh, note thatLh (euh + vh ) − Lh (euh ) = fh − Lh (euh )= −dh(19.6.25)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).• The defect dh vanishes identically at all black mesh points after a red-blackGauss-Seidel step. Thus dH = Rdh for half-weighting reduces to simplycopying half the defect from the fine grid to the corresponding coarse-gridpoint.
The calls to resid followed by rstrct in the first part of theV-cycle can be replaced by a routine that loops only over the coarse grid,filling it with half the defect.• Similarly, the quantity uenew=ueh + P veH need not be computed at redhmesh points, since they will immediately be redefined in the subsequentGauss-Seidel sweep. This means that addint need only loop over blackpoints.• You can speed up relax in several ways. First, you can have a specialform when the initial guess is zero, and omit the routine fill0. Next, youcan store h2 fh on the various grids and save a multiplication. Finally, itis possible to save an addition in the Gauss-Seidel formula by rewritingit with intermediate variables.• On typical problems, mglin with ncycle = 1 will return a solution withthe 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 toset ncycle = 2 or, more cheaply, npre = 2.
A more efficient way turnsout to be to use a higher-order P in (19.6.20) than the linear interpolationused in the V-cycle.Implementing all the above features typically gives up to a factor of twoimprovement in execution time and is certainly worthwhile in a production code.88319.6 Multigrid Methods for Boundary Value ProblemsThe right-hand side is smooth after a few nonlinear relaxation sweeps. Thus we can transferthe left-hand side to a coarse grid:LH (uH ) − LH (Reuh ) = −Rdh(19.6.26)LH (uH ) = LH (Reuh ) − Rdh(19.6.27)that is, we solvevH = ueeH − Reuh(19.6.28)uenew=ueh + P(euH − Reuh )h(19.6.29)andNote that PR 6= 1 in general, so uenew6= P ueH .
This is a key point: In equation (19.6.29) thehinterpolation error comes only from the correction, not from the full solution ueH .Equation (19.6.27) shows that one is solving for the full approximation uH , not just theerror as in the linear algorithm. This is the origin of the name FAS.The FAS multigrid algorithm thus looks very similar to the linear multigrid algorithm.The only differences are that both the defect dh and the relaxed approximation uh have to berestricted to the coarse grid, where now it is equation (19.6.27) that is solved by recursiveinvocation of the algorithm. However, instead of implementing the algorithm this way, wewill first describe the so-called dual viewpoint, which leads to a powerful alternative wayof looking at the multigrid idea.The dual viewpoint considers the local truncation error, defined asτ ≡ Lh (u) − fh(19.6.30)where u is the exact solution of the oiginal continuum equation.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
