f19-4

848 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 example, a combined advective-diffusion equation, such as ∂u ∂t=−v∂u ∂x+D∂2u ∂x2(19.3.21 ) might profitably use an explicit scheme for the advective term combined with a Crank-Nicolson or other implicit scheme for the diffusion term. The alternating-direction implicit (ADI) method, equation (19.3.16), is an example of operator splitting with a slightly different twist. Let us reinterpret(19.3.19) to have a different meaning: Let U 1now denote an updating method that includes algebraically allthe pieces of the total operator L, but which is desirably stableonly for the L1piece; likewise U2,...Um. Then a method of getting from untoun+1is un+1/m=U1(un,∆t/m) un+2/m=U2(un+1/m,∆t/m) ··· un+1=Um(un+(m−1)/m,∆t/m)(19.3.22 ) Thetimestepforeachfractionalstepin(19.3.22)isnowonly 1/mofthefulltimestep, because each partial operationacts with all the terms of the original operator. Equation(19.3.22)isusually,thoughnotalways,stableasadifferencingscheme fortheoperator L. Infact,asaruleofthumb,itisoftensufficienttohavestable Ui’s only for the operator pieces having the highest number of spatial derivatives — the otherUi’s can be unstable— to make the overall scheme stable! It is at this point that we turn our attention from initial value problems to boundaryvalueproblems. These will occupyus for the remainderof the chapter. CITED REFERENCES AND FURTHER READING: Ames, W.F. 1977, Numerical Methods for Partial Differential Equations , 2nd ed. (New York: Academic Press). 19.4 Fourier and Cyclic Reduction Methods for Boundary Value Problems As discussed in §19.0, most boundary value problems (elliptic equations, for example) reduce to solving large sparse linear systems of the form A·u=b (19.4.1 ) eitheronce,forboundaryvalueequationsthat arelinear,or iteratively,forboundary value equations that are nonlinear. 19.4FourierandCyclic ReductionMethods 849Sample 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).Two important techniques lead to “rapid” solution of equation (19.4.1) when the sparse matrix is of certain frequently occurring forms. The Fourier transform methodis directly applicable when the equations have coefficients that are constant in space. The cyclic reduction method is somewhat more general; its applicability is related to the question of whether the equations are separable (in the sense of“separation of variables”). Both methods require the boundaries to coincide with the coordinate lines. Finally, for some problems, there is a powerful combination of these two methods called FACR (Fourier Analysis and Cyclic Reduction) .W e now consider each method in turn, using equation (19.0.3), with finite-difference representation(19.0.6),asamodelexample. Generallyspeaking,themethodsinthissection are faster, when they apply, than the simpler relaxation methods discussed in§19.5; but they are not necessarily faster than the more complicated multigrid methods discussed in §19.6. FourierTransformMethod The discrete inverse Fourier transform in both xandyis ujl=1 JLJ−1/summationdisplay m=0L−1/summationdisplay n=0/hatwideumne−2πijm/Je−2πiln/L(19.4.2 ) This can be computedusingthe FFT independentlyin eachdimension,orelse all at once via the routine fournof§12.4 or the routine rlft3of§12.5. Similarly, ρjl=1 JLJ−1/summationdisplay m=0L−1/summationdisplay n=0/hatwideρmne−2πijm/Je−2πiln/L(19.4.3 ) If we substitute expressions (19.4.2) and (19.4.3) in our model problem (19.0.6), we find /hatwideumn/parenleftBig e2πim/J+e−2πim/J+e2πin/L+e−2πin/L−4/parenrightBig =/hatwideρmn∆2(19.4.4 ) or /hatwideumn=/hatwideρmn∆2 2/parenleftbigg cos2πm J+c o s2πn L−2/parenrightbigg (19.4.5 ) Thus the strategy for solving equation (19.0.6) by FFT techniques is: •Compute /hatwideρmnas the Fourier transform /hatwideρmn=J−1/summationdisplay j=0L−1/summationdisplay l=0ρjle2πimj/Je2πinl/L(19.4.6 ) •Compute /hatwideumnfrom equation (19.4.5). 850 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).•Compute ujlby the inverse Fourier transform (19.4.2). Theaboveprocedureisvalidforperiodicboundaryconditions. Inotherwords, the solution satisfies ujl=uj+J,l=uj,l+L (19.4.7 ) NextconsideraDirichletboundarycondition u=0ontherectangularboundary. Instead of the expansion (19.4.2),we now need an expansionin sine waves: ujl=2 J2 LJ−1/summationdisplay m=1L−1/summationdisplay n=1/hatwideumnsinπjm Jsinπln L(19.4.8 ) This satisfies the boundaryconditions that u=0atj=0,Jand at l=0,L.I fw e substitute this expansion and the analogous one for ρjlinto equation (19.0.6), we find that the solutionprocedureparallels that forperiodic boundaryconditions: •Compute /hatwideρmnby the sine transform /hatwideρmn=J−1/summationdisplay j=1L−1/summationdisplay l=1ρjlsinπjm Jsinπln L(19.4.9 ) (A fast sine transform algorithm was given in §12.3.) •Compute /hatwideumnfrom the expression analogous to (19.4.5), /hatwideumn=∆2/hatwideρmn 2/parenleftBig cosπm J+c o sπn L−2/parenrightBig (19.4.10 ) •Compute ujlby the inverse sine transform (19.4.8). If we have inhomogeneous boundary conditions, for example u=0on all boundariesexcept u=f(y)on the boundary x=J∆, we have to add to the above solution a solution uHof the homogeneous equation ∂2u ∂x2+∂2u ∂y2=0 ( 19.4.11 ) that satisfies the required boundary conditions. In the continuum case, this would be an expression of the form uH=/summationdisplay nAnsinhnπx L∆sinnπy L∆(19.4.12 ) where Anwould be found by requiring that u=f(y)atx=J∆. In the discrete case, we have uH jl=2 LL−1/summationdisplay n=1Ansinhπnj Lsinπnl L(19.4.13 ) 19.4FourierandCyclic ReductionMethods 851Sample 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).Iff(y=l∆)≡fl, then we get Anfrom the inverse formula An=1 sinh (πnJ/L )L−1/summationdisplay l=1flsinπnl L(19.4.14 ) The complete solution to the problem is u=ujl+uH jl (19.4.15 ) By adding appropriate terms of the form (19.4.12), we can handle inhomogeneous terms on any boundary surface. A much simpler procedure for handling inhomogeneous terms is to note that wheneverboundarytermsappearontheleft-handsideof(19.0.6),theycanbetaken over to the right-hand side since they are known. The effective source term is therefore ρjlplus a contribution from the boundary terms. To implement this idea formally, write the solution as u=u/prime+uB(19.4.16 ) where u/prime=0on the boundary, while uBvanishes everywhere excepton the boundary. There it takes on the given boundary value. In the above example, theonly nonzero values of u Bwould be uB J,l=fl (19.4.17 ) The model equation (19.0.3) becomes ∇2u/prime=−∇2uB+ρ (19.4.18 ) or, in finite-difference form, u/prime j+1,l+u/prime j−1,l+u/prime j,l+1+u/prime j,l−1−4u/prime j,l= −(uB j+1,l+uB j−1,l+uB j,l+1+uB j,l−1−4uB j,l)+∆2ρj,l(19.4.19 ) All the uBtermsin equation(19.4.19)vanishexceptwhenthe equationis evaluated atj=J−1, where u/prime J,l+u/prime J−2,l+u/prime J−1,l+1+u/prime J−1,l−1−4u/prime J−1,l=−fl+∆2ρJ−1,l(19.4.20 ) Thus the problemis now equivalentto the case of zero boundaryconditions,except that one row of the source term is modified by the replacement ∆2ρJ−1,l→∆2ρJ−1,l−fl (19.4.21 ) The case of Neumann boundary conditions ∇u=0is handled by the cosine expansion (12.3.17): ujl=2 J2 LJ/summationdisplay/prime/prime m=0L/summationdisplay/prime/prime n=0/hatwideumncosπjm Jcosπln L(19.4.22 ) 852 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).Here the double prime notation means that the terms for m=0andm=Jshould be multiplied by1 2, and similarly for n=0andn=L. Inhomogeneous terms ∇u=gcan be again included by adding a suitable solution of the homogeneous equation,or more simply by takingboundaryterms overto the right-handside. For example, the condition ∂u ∂x=g(y)atx=0 ( 19.4.23 ) becomes u1,l−u−1,l 2∆=gl (19.4.24 ) where gl≡g(y=l∆). Once again we write the solution in the form (19.4.16), where now ∇u/prime=0on the boundary. This time ∇uBtakes on the prescribed valueontheboundary,but uBvanisheseverywhereexceptjust outsidetheboundary. Thus equation (19.4.24) gives uB −1,l=−2∆gl (19.4.25 ) All the uBterms in equation (19.4.19) vanish except when j=0: u/prime 1,l+u/prime −1,l+u/prime 0,l+1+u/prime 0,l−1−4u/prime 0,l=2 ∆ gl+∆2ρ0,l (19.4.26 ) Thus u/primeis the solution of a zero-gradient problem, with the source term modified by the replacement ∆2ρ0,l→∆2ρ0,l+2 ∆ gl (19.4.27 ) Sometimes Neumann boundary conditions are handled by using a staggered grid, with the u’s defined midway between zone boundaries so that first derivatives are centered on the mesh points. You can solve such problems using similartechniques to those described above if you use the alternative form of the cosine transform, equation (12.3.23). Cyclic Reduction Evidently the FFT method works only when the original PDE has constant coefficients, and boundaries that coincide with the coordinate lines. An alternative algorithm, which can be used on somewhat more general equations, is called cyclic reduction (CR) . We illustrate cyclic reduction on the equation ∂2u ∂x2+∂2u ∂y2+b(y)∂u ∂y+c(y)u=g(x, y)( 19.4.28 ) This form arises very often in practice from the Helmholtz or Poisson equations in polar,cylindrical,orsphericalcoordinatesystems. Moregeneralseparableequations are treated in [1]. 19.4FourierandCyclic ReductionMethods 853Sample 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 finite-difference form of equation (19.4.28) can be written as a set of vector equations uj−1+T·uj+uj+1=gj∆2(19.4.29 ) Heretheindex jcomesfromdifferencinginthe x-direction,whilethe y-differencing (denoted by the index lpreviously) has been left in vector form. The matrix T has the form T=B−21 (19.4.30 ) wherethe 21comesfromthe x-differencingandthematrix Bfromthe y-differencing. The matrix B, and hence T, is tridiagonal with variable coefficients. The CR method is derived by writing down three successive equations like (19.4.29): uj−2+T·uj−1+uj=gj−1∆2 uj−1+T·uj+uj+1=gj∆2 uj+T·uj+1+uj+2=gj+1∆2(19.4.31 ) Matrix-multiplyingthe middleequationby −Tandthenaddingthethree equations, we get uj−2+T(1)·uj+uj+2=g(1) j∆2(19.4.32 ) This is an equation of the same form as (19.4.29), with T(1)=21−T2 g(1) j=∆2(gj−1−T·gj+gj+1)(19.4.33 ) After one level of CR, we have reduced the number of equations by a factor of two. Since the resultingequationsare of the same formas the originalequation,we can repeat the process. Taking the number of mesh points to be a power of 2 for simplicity,wefinally endupwitha singleequationforthe centrallineofvariables: T(f)·uJ/2=∆2g(f) J/2−u0−uJ (19.4.34 ) Here we have moved u0anduJto the right-hand side because they are known boundary values. Equation (19.4.34) can be solved for uJ/2by the standard tridiagonalalgorithm. Thetwoequationsatlevel f−1involveuJ/4andu3J/4. The equationfor uJ/4involvesu0anduJ/2, bothofwhich are known,andhencecan be solved by the usual tridiagonal routine. A similar result holds true at every stage, so we end up solving J−1tridiagonal systems. Inpractice,equations(19.4.33)shouldberewrittentoavoidnumericalinstabil- ity. For these and other practical details, refer to [2]. 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 )