f19-3

844 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).The correct way to difference Schr ¨odinger’s equation [1,2]is to use Cayley’s formforthefinite-differencerepresentationof e−iHt,whichissecond-orderaccurate andunitary: e−iHt/similarequal1−1 2iH∆t 1+1 2iH∆t(19.2.35 ) In other words, /parenleftbig 1+1 2iH∆t/parenrightbig ψn+1 j=/parenleftbig 1−1 2iH∆t/parenrightbig ψn j (19.2.36 ) On replacing Hby its finite-difference approximation in x, we have a complex tridiagonalsystemtosolve. Themethodisstable,unitary,andsecond-orderaccurate in space and time. In fact, it is simply the Crank-Nicolsonmethod onceagain! CITED REFERENCES AND FURTHER READING: Ames, W.F. 1977, Numerical Methods for Partial Differential Equations , 2nd ed. (New York: Academic Press), Chapter 2. Goldberg, A., Schey, H.M., and Schwartz, J.L. 1967, American Journal of Physics , vol. 35, pp. 177–186. [1] Galbraith, I., Ching, Y.S., and Abraham, E. 1984, American Journal of Physics , vol. 52, pp. 60– 68. [2] 19.3 InitialValue Problems in Multidimensions The methods described in §19.1 and §19.2 for problems in 1+1dimension (onespaceandonetimedimension)caneasily begeneralizedto N+1dimensions. However, the computing power necessary to solve the resulting equations is enor- mous. If you have solved a one-dimensional problem with 100 spatial grid points, solving the two-dimensional version with 100 ×100 mesh points requires at least 100 times as much computing. You generally have to be content with very modest spatial resolution in multidimensional problems. Indulge us in offering a bit of advice about the development and testing of multidimensional PDE codes: You should always first run your programs on very smallgrids, e.g., 8×8, even though the resulting accuracy is so poor as to be useless. When yourprogramis all debuggedand demonstrablystable, thenyoucan increase the grid size to a reasonable one and start looking at the results. We have actually heard someone protest, “my program would be unstable for a crude grid,but I am sure the instability will go away on a larger grid.” That is nonsense of a most pernicious sort, evidencing total confusion between accuracy and stability. In fact, new instabilities sometimes do show up on largergrids; but old instabilities never (in our experience) just go away. Forcedtolivewithmodestgridsizes,somepeoplerecommendgoingtohigher- ordermethodsinanattempttoimproveaccuracy. Thisisverydangerous. Unlessthe solution you are lookingfor is knownto be smooth, and the high-ordermethodyou 19.3InitialValueProblemsinMultidimensions 845Sample 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).are using is known to be extremely stable, we do not recommend anything higher thansecond-orderintime(forsetsoffirst-orderequations). Forspatialdifferencing,we recommend the order of the underlying PDEs, perhaps allowing second-order spatial differencing for first-order-in-space PDEs. When you increase the order of a differencingmethod to greater than the order of the original PDEs, you introducespurioussolutionstothedifferenceequations. Thisdoesnotcreateaproblemifthey allhappentodecayexponentially;otherwiseyouaregoingtoseeallhellbreakloose! Lax Method for a Flux-Conservative Equation As an example, we show how to generalize the Lax method (19.1.15) to two dimensions for the conservation equation ∂u ∂t=−∇ ·F=−/parenleftbigg∂F x ∂x+∂F y ∂y/parenrightbigg (19.3.1 ) Use a spatial grid with xj=x0+j∆ yl=y0+l∆(19.3.2 ) We have chosen ∆x=∆y≡∆for simplicity. Then the Lax scheme is un+1 j,l=1 4(un j+1 ,l+un j−1,l+un j,l+1+un j,l−1) −∆t 2∆(Fn j+1 ,l−Fn j−1,l+Fn j,l+1−Fn j,l−1)(19.3.3 ) Note that as an abbreviated notation Fj+1andFj−1refer to Fx, while Fl+1and Fl−1refer to Fy. Let us carry out a stability analysis for the model advective equation (analog of 19.1.6) with Fx=vxu, F y=vyu (19.3.4 ) Thisrequiresaneigenmodewithtwodimensionsinspace,thoughstillonlyasimple dependence on powers of ξin time, un j,l=ξneik xj∆eik yl∆(19.3.5 ) Substituting in equation (19.3.3), we find ξ=1 2(coskx∆+c o s ky∆)−iαxsinkx∆−iαysinky∆( 19.3.6 ) where αx=vx∆t ∆,α y=vy∆t ∆(19.3.7 ) 846 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).The expression for |ξ|2can be manipulated into the form |ξ|2=1−(sin2kx∆+s i n2ky∆)/bracketleftbigg1 2−(α2 x+α2 y)/bracketrightbigg −1 4(coskx∆−cosky∆)2−(αysinkx∆−αxsinky∆)2(19.3.8 ) Thelast two termsare negative,and so thestability requirement |ξ|2≤1becomes 1 2−(α2 x+α2 y)≥0( 19.3.9 ) or ∆t≤∆√ 2(v2x+v2y)1/2(19.3.10 ) This is an example of the general result for the N-dimensional Courant condition: If |v|is the maximum propagationvelocity in the problem, then ∆t≤∆√ N|v|(19.3.11 ) is the Courant condition. DiffusionEquation in Multidimensions Let us consider the two-dimensional diffusion equation, ∂u ∂t=D/parenleftbigg∂2u ∂x2+∂2u ∂y2/parenrightbigg (19.3.12 ) An explicit method, such as FTCS, can be generalized from the one-dimensional case in the obviousway. However,we haveseen that diffusiveproblemsare usuallybest treated implicitly. Supposewe try to implement the Crank-Nicolsonscheme in two dimensions. This would give us u n+1 j,l=un j,l+1 2α/parenleftBig δ2 xun+1 j,l+δ2 xunj,l+δ2 yun+1 j,l+δ2 yunj,l/parenrightBig (19.3.13 ) Here α≡D∆t ∆2∆≡∆x=∆y (19.3.14 ) δ2 xunj,l≡un j+1 ,l−2un j,l+un j−1,l (19.3.15 ) and similarly for δ2 yun j,l. This is certainly a viable scheme; the problem arises in solving the coupled linear equations. Whereas in one space dimension the system was tridiagonal, that is no longer true, though the matrix is still very sparse. Onepossibility is to use a suitable sparse matrix technique (see §2.7 and §19.0). Another possibility, which we generally prefer, is a slightly different way of generalizing the Crank-Nicolson algorithm. It is still second-orderaccurate in time and space, and unconditionally stable, but the equations are easier to solve than 19.3InitialValueProblemsinMultidimensions 847Sample 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).(19.3.13). Called the alternating-directionimplicitmethod(ADI) ,this embodiesthe powerful concept of operator splitting ortime splitting , about which we will say more below. Here, the idea is to divide each timestep into two steps of size ∆t/2. In each substep, a different dimension is treated implicitly: un+1 /2 j,l=un j,l+1 2α/parenleftBig δ2 xun+1 /2 j,l+δ2 yunj,l/parenrightBig un+1 j,l=un+1 /2 j,l+1 2α/parenleftBig δ2 xun+1 /2 j,l+δ2 yun+1 j,l/parenrightBig (19.3.16 ) The advantage of this method is that each substep requires only the solution of a simple tridiagonal system. OperatorSplitting Methods Generally The basic idea of operator splitting, which is also called time splitting orthe method of fractional steps , is this: Suppose you have an initial value equation of the form ∂u ∂t=Lu (19.3.17 ) where Lis some operator. While Lis not necessarily linear, suppose that it can at least be written as a linear sum of mpieces, which act additively on u, Lu=L1u+L2u+···+Lmu (19.3.18 ) Finally,supposethatfor eachofthepieces,youalreadyknowadifferencingscheme for updating the variable ufrom timestep nto timestep n+1, valid if that piece of the operator were the onlyone on the right-hand side. We will write these updatings symbolically as un+1=U1(un,∆t) un+1=U2(un,∆t) ··· un+1=Um(un,∆t)(19.3.19 ) Now, one form of operator splitting would be to get from nton+1by the following sequence of updatings: un+(1 /m)=U1(un,∆t) un+(2 /m)=U2(un+(1 /m),∆t) ··· un+1=Um(un+(m−1)/m,∆t)(19.3.20 ) 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.