f19-5

854 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FACRMethod Thebestway to solve equations of the form (19.4.28), including the constant coefficientproblem(19.0.3),isacombinationofFourieranalysisandcyclicreduction, theFACRmethod [3-6]. Ifatthe rthstageofCRweFourieranalyzetheequationsof the form (19.4.32) along y, that is, with respect to the suppressed vector index, we will have a tridiagonal system in the x-direction for each y-Fourier mode: /hatwideuk j−2r+λ(r) k/hatwideuk j+/hatwideuk j+2r=∆2g(r)k j (19.4.35 ) Here λ(r) kis the eigenvalue of T(r)corresponding to the kth Fourier mode. For the equation (19.0.3), equation (19.4.5) shows that λ(r) kwill involve terms like cos(2 πk/L )−2raisedtoapower. Solvethetridiagonalsystemsfor /hatwideuk jatthelevels j=2r,2×2r,4×2r, ..., J−2r. Fourier synthesize to get the y-values on these x-lines. Then fill in the intermediate x-lines as in the original CR algorithm. The trick is to choose the number of levels of CR so as to minimize the total numberofarithmeticoperations. Onecanshowthatforatypicalcaseofa 128×128 mesh, the optimal level is r=2; asymptotically, r→log2(log2J). A rough estimate of running times for these algorithms for equation (19.0.3) is as follows: The FFT method (in both xandy) and the CR method are roughly comparable. FACR with r=0(that is, FFT in one dimension and solve the tridiagonal equations by the usual algorithm in the other dimension) gives about a factor of two gain in speed. The optimal FACR with r=2gives another factor of two gain in speed. CITED REFERENCES AND FURTHER READING: Swartzrauber, P.N. 1977, SIAM Review , vol. 19, pp. 490–501. [1] Buzbee,B.L, Golub,G.H., andNielson,C.W.1970, SIAM JournalonNumericalAnalysis ,vol.7, pp. 627–656; see also op. cit.vol. 11, pp. 753–763. [2] Hockney,R.W.1965, JournaloftheAssociationforComputingMachinery ,vol.12,pp.95–113.[3] Hockney, R.W.1970, in Methods ofComputational Physics , vol. 9(NewYork: Academic Press), pp. 135–211. [4] Hockney, R.W., and Eastwood, J.W. 1981, Computer Simulation Using Particles (New York: McGraw-Hill), Chapter 6. [5] Temperton, C. 1980, Journal of Computational Physics , vol. 34, pp. 314–329. [6] 19.5 Relaxation Methods for Boundary Value Problems As we mentioned in §19.0, relaxation methods involve splitting the sparse matrixthatarisesfromfinitedifferencingandtheniteratinguntilasolutionisfound. There is another way of thinking about relaxation methods that is somewhat more physical. Suppose we wish to solve the elliptic equation Lu=ρ (19.5.1 ) 19.5RelaxationMethodsforBoundaryValueProblems 855Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).where Lrepresents some elliptic operator and ρis the source term. Rewrite the equation as a diffusion equation, ∂u ∂t=Lu−ρ (19.5.2 ) An initial distribution urelaxesto an equilibrium solution as t→∞. This equilibrium has all time derivatives vanishing. Therefore it is the solution of the originalellipticproblem(19.5.1). Weseethatallthemachineryof §19.2,ondiffusive initial value equations, can be brought to bear on the solution of boundary valueproblems by relaxation methods. Let usapplythis ideato ourmodelproblem(19.0.3). Thediffusionequationis ∂u ∂t=∂2u ∂x2+∂2u ∂y2−ρ (19.5.3 ) If we use FTCS differencing (cf. equation 19.2.4), we get un+1 j,l=un j,l+∆t ∆2/parenleftbig un j+1 ,l+un j−1,l+un j,l+1+un j,l−1−4un j,l/parenrightbig −ρj,l∆t(19.5.4 ) Recallfrom(19.2.6)thatFTCSdifferencingisstableinonespatialdimensiononlyif ∆t/∆2≤1 2. In two dimensions this becomes ∆t/∆2≤1 4. Supposewe try to take thelargestpossible timestep,andset ∆t=∆2/4. Thenequation(19.5.4)becomes un+1 j,l=1 4/parenleftbig un j+1 ,l+un j−1,l+un j,l+1+un j,l−1/parenrightbig −∆2 4ρj,l (19.5.5 ) Thus the algorithm consists of using the average of uat its four nearest-neighbor points on the grid (plus the contribution from the source). This procedure is then iterated until convergence. This method is in fact a classical method with origins dating back to the last century, called Jacobi’s method (not to be confused with the Jacobi method for eigenvalues). The method is not practical because it converges too slowly.However, it is the basis for understanding the modern methods, which are always compared with it. Another classical method is the Gauss-Seidel method, which turns out to be importantinmultigridmethods( §19.6). Herewemakeuseofupdatedvaluesof uon theright-handsideof(19.5.5)as soonas theybecomeavailable. Inotherwords,theaveragingis done “in place” instead of being “copied” from an earlier timestep to a later one. Ifwe are proceedingalongthe rows,incrementing jforfixed l,weha v e u n+1 j,l=1 4/parenleftBig un j+1 ,l+un+1 j−1,l+un j,l+1+un+1 j,l−1/parenrightBig −∆2 4ρj,l (19.5.6 ) Thismethodis also slowlyconvergingandonlyoftheoreticalinterestwhenusedby itself, but some analysis of it will be instructive. Let us look at the Jacobi and Gauss-Seidel methods in terms of the matrix splitting concept. We change notationand call u“x,” to conformto standard matrix notation. To solve A·x=b (19.5.7 ) 856 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).we can consider splitting Aas A=L+D+U (19.5.8 ) whereDis the diagonal part of A,Lis the lower triangle of Awith zeros on the diagonal, and Uis the upper triangle of Awith zeros on the diagonal. In the Jacobi method we write for the rth step of iteration D·x(r)=−(L+U)·x(r−1)+b (19.5.9 ) Forourmodelproblem(19.5.5), Dis simplytheidentitymatrix. TheJacobimethod converges for matrices Athat are “diagonally dominant” in a sense that can be made mathematically precise. For matrices arising from finite differencing, this condition is usually met. What is the rate of convergence of the Jacobi method? A detailed analysis is beyond our scope, but here is some of the flavor: The matrix −D−1·(L+U)is theiterationmatrix which, apart from an additiveterm, maps one set of x’s into the next. The iteration matrix has eigenvalues, each one of which reflects the factor by which the amplitude of a particular eigenmode of undesired residual is suppressed during one iteration. Evidently those factors had better all have modulus <1for the relaxation to work at all! The rate of convergence of the method is set by the rate for the slowest-decaying eigenmode, i.e., the factor with largest modulus. The modulusofthis largestfactor,thereforelyingbetween0 and1, is calledthe spectral radiusof the relaxation operator, denoted ρs. The number of iterations rrequired to reduce the overall error by a factor 10−pis thus estimated by r≈pln 10 (−lnρs)(19.5.10 ) In general, the spectral radius ρsgoes asymptotically to the value 1as the grid sizeJis increased, so that more iterations are required. For any given equation, grid geometry, and boundary condition , the spectral radius can, in principle, be computed analytically. For example, for equation (19.5.5) on a J×Jgrid with Dirichlet boundary conditions on all four sides, the asymptotic formula for large Jturns out to be ρs/similarequal1−π2 2J2(19.5.11 ) The numberofiterations rrequiredto reducethe errorbya factorof 10−pis thus r/similarequal2pJ2ln 10 π2/similarequal1 2pJ2(19.5.12 ) Inotherwords,thenumberofiterationsisproportionaltothenumberofmeshpoints, J2. Since 100×100and larger problems are common, it is clear that the Jacobi method is only of theoretical interest. 19.5RelaxationMethodsforBoundaryValueProblems 857Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).The Gauss-Seidel method, equation (19.5.6), corresponds to the matrix de- composition (L+D)·x(r)=−U·x(r−1)+b (19.5.13 ) The fact that Lis on the left-hand side of the equation follows from the updating in place, as you can easily check if you write out (19.5.13) in components. Onecan show [1-3]that the spectral radius is just the square of the spectral radius of the Jacobi method. For our model problem, therefore, ρs/similarequal1−π2 J2(19.5.14 ) r/similarequalpJ2ln 10 π2/similarequal1 4pJ2(19.5.15 ) The factor of two improvement in the number of iterations over the Jacobi method still leaves the method impractical. Successive Overrelaxation (SOR) Wegetabetteralgorithm—onethatwasthestandardalgorithmuntilthe1970s —ifwemak ean overcorrection tothevalueof x(r)at the rth stageofGauss-Seidel iteration, thus anticipating future corrections. Solve (19.5.13) for x(r), add and subtractx(r−1)onthe right-handside,andhencewritethe Gauss-Seidelmethodas x(r)=x(r−1)−(L+D)−1·[(L+D+U)·x(r−1)−b]( 19.5.16 ) The term in square brackets is just the residual vector ξ(r−1),s o x(r)=x(r−1)−(L+D)−1·ξ(r−1)(19.5.17 ) Nowovercorrect , defining x(r)=x(r−1)−ω(L+D)−1·ξ(r−1)(19.5.18 ) Here ωis called the overrelaxation parameter , and the method is called successive overrelaxation (SOR). The following theorems can be proved [1-3]: •The method is convergent only for 0<ω< 2.I f0<ω< 1, we speak ofunderrelaxation . •Under certain mathematical restrictions generally satisfied by matrices arisingfromfinitedifferencing,onlyoverrelaxation( 1<ω< 2)cangive faster convergence than the Gauss-Seidel method. •IfρJacobiis the spectral radius of the Jacobi iteration (so that the square of it is the spectral radius of the Gauss-Seidel iteration), then the optimal choice for ωis given by ω=2 1+/radicalbig 1−ρ2 Jacobi(19.5.19 ) 858 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).•For this optimal choice, the spectral radius for SOR is ρSOR =/parenleftBigg ρJacobi 1+/radicalbig 1−ρ2 Jacobi/parenrightBigg2 (19.5.20 ) As an application of the above results, consider our model problem for which ρJacobiis givenbyequation(19.5.11). Thenequations(19.5.19)and(19.5.20)give ω/similarequal2 1+π/J(19.5.21 ) ρSOR/similarequal1−2π Jforlarge J (19.5.22 ) Equation (19.5.10) gives for the number of iterations to reduce the initial error by a factor of 10−p, r/similarequalpJln 10 2π/similarequal1 3pJ (19.5.23 ) Comparing with equation (19.5.12) or (19.5.15), we see that optimal SOR requires of order Jiterations, as opposed to of order J2. Since Jis typically 100 or larger, this makes a tremendous difference! Equation (19.5.23) leads to the mnemonic that 3-figure accuracy ( p=3) requires a number of iterations equal to the number of mesh points along a side of the grid. For 6-figure accuracy, we require about twice as many iterations. How do we choose ωfor a problem for which the answer is not known analytically? That is just the weak point of SOR! The advantages of SOR obtain only 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 known problem,replacingthecoefficientsintheequationbyaveragevalues. Note,however, thattheknownproblemmusthavethesamegridsizeandboundaryconditionsasthe actual problem. We give for reference purposes the value of ρJacobifor our model problemon a rectangular J×Lgrid,allowingforthe possibility that ∆x/negationslash=∆y: ρJacobi =cosπ J+/parenleftbigg∆x ∆y/parenrightbigg2 cosπ L 1+/parenleftbigg∆x ∆y/parenrightbigg2(19.5.24 ) Equation (19.5.24) holds for homogeneous Dirichlet or Neumann boundary condi- tions. For periodic boundary conditions, make the replacement π→2π. A second way, which is especially useful if you plan to solve many similar elliptic equations each time with slightly different coefficients, is to determine the optimum value ωempirically on the first equation and then use that value for the remaining equations. Various automated schemes for doing this and for “seeking out” 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. 19.5RelaxationMethodsforBoundaryValueProblems 859Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Consider a general second-order elliptic equation in xandy, finite differenced on a square as for our model equation. Corresponding to each row of the matrix A is an equation of the form aj,luj+1 ,l+bj,luj−1,l+cj,luj,l+1+dj,luj,l−1+ej,luj,l=fj,l (19.5.25 ) For our model equation, we had a=b=c=d=1,e=−4. The quantity fis proportional to the source term. The iterative procedure is defined by solving (19.5.25) for uj,l: u*j,l=1 ej,l(fj,l−aj,luj+1 ,l−bj,luj−1,l−cj,luj,l+1−dj,luj,l−1)(19.5.26 ) Then unew j,lis a weighted average unew j,l=ωu*j,l+( 1−ω)uold j,l (19.5.27 ) We calculate it as follows: The residual at any stage is ξj,l=aj,luj+1 ,l+bj,luj−1,l+cj,luj,l+1+dj,luj,l−1+ej,luj,l−fj,l(19.5.28 ) and the SOR algorithm (19.5.18) or (19.5.27) is unew j,l=uold j,l−ωξj,l ej,l(19.5.29 ) This formulation is very easy to program, and the norm of the residual vector ξj,l can be used as a criterion for terminating the iteration. Anotherpracticalpointconcernstheorderinwhichmeshpointsareprocessed. 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 the red and black squares of a checkerboard. Then equation (19.5.26) shows that the odd 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 another half-sweep updating the even points with the new odd values. For the version of SOR implemented below, we shall adopt odd-even ordering. The last practical point is that in practice the asymptotic rate of convergence in SOR is not attained until of order Jiterations. The error often grows by a factor of 20 before convergencesets in. A trivial modification to SOR resolves this problem. It is based on the observation that, while ωis the optimum asymptotic relaxation parameter, it is not necessarily a good initial choice. In SOR withChebyshev acceleration , one uses odd-even ordering and changes ωat each half- sweep according to the following prescription: ω (0)=1 ω(1/2)=1/(1−ρ2 Jacobi /2) ω(n+1 /2)=1/(1−ρ2 Jacobi ω(n)/4),n =1/2,1, ...,∞ ω(∞)→ωoptimal(19.5.30 ) 860 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ThebeautyofChebyshevaccelerationis thatthenormoftheerroralwaysdecreases with each iteration. (This is the norm of the actual error in uj,l. The norm of the residual ξj,lneed not decrease monotonically.) While the asymptotic rate of convergence is the same as ordinary SOR, there is never any excuse for not using Chebyshev acceleration to reduce the total number of iterations required. Here we give a routine for SOR with Chebyshev acceleration. SUBROUTINE sor(a,b,c,d,e,f,u,jmax,rjac) INTEGER jmax,MAXITSDOUBLE PRECISION rjac,a(jmax,jmax),b(jmax,jmax), * c(jmax,jmax),d(jmax,jmax),e(jmax,jmax), * f(jmax,jmax),u(jmax,jmax),EPS PARAMETER (MAXITS=1000,EPS=1.d-5) Successive overrelaxation solution of equation (19.5.25) with Chebyshev acceleration. a, b,c,d,e,a n d fare input as the coefficients of the equation, each dimensioned to the grid size JMAX×JMAX.uis input as the initial guess to the solution, usually zero, and returns with the final value. rjacis input as the spectral radius of the Jacobi iteration, or an estimate of it. INTEGER ipass,j,jsw,l,lsw,nDOUBLE PRECISION anorm,anormf, * omega,resid Doubleprecisionisagoodideafor JMAXbiggerthanabout25. anormf=0.d0 Compute initial norm of residual and terminate iteration when normhasbeenreducedbyafactor EPS. do 12j=2,jmax-1 do11l=2,jmax-1 anormf=anormf+abs(f(j,l)) Assumesinitial uiszero. enddo 11 enddo 12 omega=1.d0 do16n=1,MAXITS anorm=0.d0jsw=1 do 15ipass=1,2 Odd-evenordering. lsw=jswdo 14j=2,jmax-1 do13l=lsw+1,jmax-1,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=anorm+abs(resid) u(j,l)=u(j,l)-omega*resid/e(j,l) enddo 13 lsw=3-lsw enddo 14 jsw=3-jsw if(n.eq.1.and.ipass.eq.1) then omega=1.d0/(1.d0-.5d0*rjac**2) else omega=1.d0/(1.d0-.25d0*rjac**2*omega) endif enddo 15 if(anorm.lt.EPS*anormf)return enddo 16 pause ’MAXITS exceeded in sor’ END The main advantage of SOR is that it is very easy to program. Its main disadvantage is that it is still very inefficient on large problems. 19.5RelaxationMethodsforBoundaryValueProblems 861Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ADI(Alternating-Direction Implicit)Method TheADImethodof §19.3fordiffusionequationscanbeturnedintoarelaxation method for elliptic equations [1-4].I n§19.3, we discussed ADI as a method for solving the time-dependent heat-flow equation ∂u ∂t=∇2u−ρ (19.5.31 ) By letting t→∞onealso gets aniterativemethodforsolvingtheelliptic equation ∇2u=ρ (19.5.32 ) In either case, the operator splitting is of the form L=Lx+Ly (19.5.33 ) where Lxrepresents the differencing in xandLythat in y. For example,in our modelproblem(19.0.6)with ∆x=∆y=∆,w eh a v e Lxu=2uj,l−uj+1 ,l−uj−1,l Lyu=2uj,l−uj,l+1−uj,l−1(19.5.34 ) More 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. Usually one tries to base the splitting on the physical nature of the problem. We know for our model problem that an initial transient diffuses away, and we set up the xand ysplitting to mimic diffusion in each dimension. Havingchosenasplitting,wedifferencethetime-dependentequation(19.5.31) implicitly in two half-steps: un+1 /2−un ∆t/2=−Lxun+1 /2+Lyun ∆2−ρ un+1−un+1 /2 ∆t/2=−Lxun+1 /2+Lyun+1 ∆2−ρ(19.5.35 ) (cf. equation 19.3.16). Here we have suppressed the spatial indices ( j, l). In matrix notation, equations (19.5.35) are (Lx+r1)·un+1 /2=(r1−Ly)·un−∆2ρ (19.5.36 ) (Ly+r1)·un+1=(r1−Lx)·un+1 /2−∆2ρ (19.5.37 ) where r≡2∆2 ∆t(19.5.38 ) The matrices on the left-hand sides of equations (19.5.36) and (19.5.37) are tridiagonal (and usually positive definite), so the equations can be solved by the 862 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).standardtridiagonalalgorithm. Given un,onesolves(19.5.36)for un+1 /2,substitutes on the right-hand side of (19.5.37), and then solves for un+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 rfor 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 computertime. 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 LargeLinearSystems (NewYork: Academic Press). [1] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork: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 Practicalmultigridmethodswerefirstintroducedinthe1970sbyBrandt. These methods can solve elliptic PDEs discretized on Ngrid points in O(N)operations. The “rapid” direct elliptic solvers discussed in §19.4 solve special kinds of elliptic equations in O(NlogN)operations. The numerical coefficients in these estimates are such that multigrid methods are comparable to the rapid methods in executionspeed. Unlike the rapid methods, however,the multigridmethods 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 solvingtheseproblems. Youhaveto adjustthevariouscomponentsofthealgorithm within this framework to solve your specific problem. We can only give a brief