c19-5 (779620), страница 2
Текст из файла (страница 2)
Onecan show [1-3] that the spectral radius is just the square of the spectral radius of theJacobi method. For our model problem, therefore,19.5 Relaxation Methods for Boundary Value Problems867• For this optimal choice, the spectral radius for SOR isρSOR =ρJacobip1 + 1 − ρ2Jacobi!2(19.5.20)21 + π/J2π'1−Jω'ρSOR(19.5.21)for large J(19.5.22)Equation (19.5.10) gives for the number of iterations to reduce the initial error bya factor of 10−p ,r'1pJ ln 10' pJ2π3(19.5.23)Comparing with equation (19.5.12) or (19.5.15), we see that optimal SOR requiresof order J iterations, as opposed to of order J 2 . Since J is typically 100 or larger,this makes a tremendous difference! Equation (19.5.23) leads to the mnemonicthat 3-figure accuracy (p = 3) requires a number of iterations equal to the numberof mesh points along a side of the grid. For 6-figure accuracy, we require abouttwice as many iterations.How do we choose ω for a problem for which the answer is not knownanalytically? That is just the weak point of SOR! The advantages of SOR obtainonly in a fairly narrow window around the correct value of ω.
It is better to take ωslightly too large, rather than slightly too small, but best to get it right.One way to choose ω is to map your problem approximately onto a knownproblem, replacing the coefficients in the equation by average values. Note, however,that the known problem must have the same grid size and boundary conditions as theactual problem. We give for reference purposes the value of ρJacobi for our modelproblem on a rectangular J × L grid, allowing for the possibility that ∆x 6= ∆y:cosρJacobi =2π∆xcos∆yL2∆x1+∆yπ+J(19.5.24)Equation (19.5.24) holds for homogeneous Dirichlet or Neumann boundary conditions.
For periodic boundary conditions, make the replacement π → 2π.A second way, which is especially useful if you plan to solve many similarelliptic equations each time with slightly different coefficients, is to determine theoptimum value ω empirically on the first equation and then use that value for theremaining equations.
Various automated schemes for doing this and for “seekingout” the best values of ω are described in the literature.While the matrix notation introduced earlier is useful for theoretical analyses,for practical implementation of the SOR algorithm we need explicit formulas.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).As an application of the above results, consider our model problem for whichρJacobi is given by equation (19.5.11). Then equations (19.5.19) and (19.5.20) give868Chapter 19.Partial Differential EquationsConsider a general second-order elliptic equation in x and y, finite differenced ona square as for our model equation.
Corresponding to each row of the matrix Ais an equation of the formaj,l uj+1,l + bj,l uj−1,l + cj,l uj,l+1 + dj,l uj,l−1 + ej,l uj,l = fj,l(19.5.25)u*j,l =1(fj,l − aj,l uj+1,l − bj,l uj−1,l − cj,l uj,l+1 − dj,l uj,l−1 )ej,l(19.5.26)Then unewj,l is a weighted averageoldunewj,l = ωu*j,l + (1 − ω)uj,l(19.5.27)We calculate it as follows: The residual at any stage isξj,l = aj,l uj+1,l + bj,l uj−1,l + cj,l uj,l+1 + dj,l uj,l−1 + ej,l uj,l − fj,l (19.5.28)and the SOR algorithm (19.5.18) or (19.5.27) isoldunewj,l = uj,l − ωξj,lej,l(19.5.29)This formulation is very easy to program, and the norm of the residual vector ξj,lcan be used as a criterion for terminating the iteration.Another practical point concerns the order in which mesh points are processed.The obvious strategy is simply to proceed in order down the rows (or columns).Alternatively, suppose we divide the mesh into “odd” and “even” meshes, like thered and black squares of a checkerboard.
Then equation (19.5.26) shows that theodd points depend only on the even mesh values and vice versa. Accordingly,we can carry out one half-sweep updating the odd points, say, and then anotherhalf-sweep updating the even points with the new odd values. For the version ofSOR implemented below, we shall adopt odd-even ordering.The last practical point is that in practice the asymptotic rate of convergencein SOR is not attained until of order J iterations. The error often grows by afactor of 20 before convergence sets in.
A trivial modification to SOR resolves thisproblem. It is based on the observation that, while ω is the optimum asymptoticrelaxation parameter, it is not necessarily a good initial choice. In SOR withChebyshev acceleration, one uses odd-even ordering and changes ω at each halfsweep according to the following prescription:ω(0) = 1ω(1/2) = 1/(1 − ρ2Jacobi /2)ω(n+1/2) = 1/(1 − ρ2Jacobi ω(n) /4),ω(∞) → ωoptimaln = 1/2, 1, ..., ∞(19.5.30)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 our model equation, we had a = b = c = d = 1, e = −4. The quantityf is proportional to the source term. The iterative procedure is defined by solving(19.5.25) for uj,l :19.5 Relaxation Methods for Boundary Value Problems869The beauty of Chebyshev acceleration is that the norm of the error always decreaseswith each iteration. (This is the norm of the actual error in uj,l. The norm ofthe residual ξj,l need not decrease monotonically.) While the asymptotic rate ofconvergence is the same as ordinary SOR, there is never any excuse for not usingChebyshev acceleration to reduce the total number of iterations required.#include <math.h>#define MAXITS 1000#define EPS 1.0e-5void sor(double **a, double **b, double **c, double **d, double **e,double **f, double **u, int jmax, double rjac)Successive overrelaxation solution of equation (19.5.25) with Chebyshev acceleration.
a, b, c,d, e, and f are input as the coefficients of the equation, each dimensioned to the grid size[1..jmax][1..jmax]. u is input as the initial guess to the solution, usually zero, and returnswith the final value. rjac is input as the spectral radius of the Jacobi iteration, or an estimateof it.{void nrerror(char error_text[]);int ipass,j,jsw,l,lsw,n;double anorm,anormf=0.0,omega=1.0,resid;Double precision is a good idea for jmax bigger than about 25.for (j=2;j<jmax;j++)Compute initial norm of residual and terminate iteration when norm has been reduced bya factor EPS.for (l=2;l<jmax;l++)anormf += fabs(f[j][l]);Assumes initial u is zero.for (n=1;n<=MAXITS;n++) {anorm=0.0;jsw=1;for (ipass=1;ipass<=2;ipass++) {Odd-even ordering.lsw=jsw;for (j=2;j<jmax;j++) {for (l=lsw+1;l<jmax;l+=2) {resid=a[j][l]*u[j+1][l]+b[j][l]*u[j-1][l]+c[j][l]*u[j][l+1]+d[j][l]*u[j][l-1]+e[j][l]*u[j][l]-f[j][l];anorm += fabs(resid);u[j][l] -= omega*resid/e[j][l];}lsw=3-lsw;}jsw=3-jsw;omega=(n == 1 && ipass == 1 ? 1.0/(1.0-0.5*rjac*rjac) :1.0/(1.0-0.25*rjac*rjac*omega));}if (anorm < EPS*anormf) return;}nrerror("MAXITS exceeded");}The main advantage of SOR is that it is very easy to program.
Its maindisadvantage is that it is still very inefficient on large problems.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).Here we give a routine for SOR with Chebyshev acceleration.870Chapter 19.Partial Differential EquationsADI (Alternating-Direction Implicit) MethodThe ADI method of §19.3 for diffusion equations can be turned into a relaxationmethod for elliptic equations [1-4] .
In §19.3, we discussed ADI as a method forsolving the time-dependent heat-flow equation(19.5.31)By letting t → ∞ one also gets an iterative method for solving the elliptic equation∇2 u = ρ(19.5.32)In either case, the operator splitting is of the formL = Lx + Ly(19.5.33)where Lx represents the differencing in x and Ly that in y.For example, in our model problem (19.0.6) with ∆x = ∆y = ∆, we haveLxu = 2uj,l − uj+1,l − uj−1,l(19.5.34)Ly u = 2uj,l − uj,l+1 − uj,l−1More complicated operators may be similarly split, but there is some art involved.A bad choice of splitting can lead to an algorithm that fails to converge.
Usuallyone tries to base the splitting on the physical nature of the problem. We know forour model problem that an initial transient diffuses away, and we set up the x andy splitting to mimic diffusion in each dimension.Having chosen a splitting, we difference the time-dependent equation (19.5.31)implicitly in two half-steps:Lx un+1/2 + Ly unun+1/2 − un=−−ρ∆t/2∆2Lx un+1/2 + Ly un+1un+1 − un+1/2=−−ρ∆t/2∆2(19.5.35)(cf.
equation 19.3.16). Here we have suppressed the spatial indices (j, l). In matrixnotation, equations (19.5.35) are(Lx + r1) · un+1/2 = (r1 − Ly ) · un − ∆2 ρ(Ly + r1) · un+1= (r1 − Lx) · un+1/2(19.5.36)−∆ ρ2(19.5.37)wherer≡2∆2∆t(19.5.38)The matrices on the left-hand sides of equations (19.5.36) and (19.5.37) aretridiagonal (and usually positive definite), so the equations can be solved by theSample 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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
