f19-2

838 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).Woodward,P., andColella,P. 1984, Journalof ComputationalPhysics , vol. 54,pp. 115–173.[8] Rizzi, A., and Engquist, B. 1987, Journal of Computational Physics , vol. 72, pp. 1–69. [9] 19.2 Diffusive Initial Value Problems Recall the model parabolic equation, the diffusion equation in one space dimension, ∂u ∂t=∂ ∂x/parenleftbigg D∂u ∂x/parenrightbigg (19.2.1 ) where Dis the diffusion coefficient. Actually, this equation is a flux-conservative equation of the form considered in the previous section, with F=−D∂u ∂x(19.2.2 ) the flux in the x-direction. We will assume D≥0, otherwise equation (19.2.1)has physicallyunstablesolutions: Asmalldisturbanceevolvestobecomemoreandmoreconcentratedinsteadofdispersing. (Don’tmakethemistakeoftryingtofindastable differencingschemeforaproblemwhoseunderlyingPDEsarethemselvesunstable!) Even though (19.2.1)is of the form already considered, it is useful to consider it as a model in its own right. The particular form of flux (19.2.2), and its direct generalizations, occur quite frequentlyin practice. Moreover,we have already seenthat numerical viscosity and artificial viscosity can introduce diffusive pieces like the right-hand side of (19.2.1) in many other situations. Consider first the case when Dis a constant. Then the equation ∂u ∂t=D∂2u ∂x2(19.2.3 ) can be differenced in the obvious way: un+1 j−un j ∆t=D/bracketleftbiggun j+1−2un j+un j−1 (∆x)2/bracketrightbigg (19.2.4 ) This is the FTCS scheme again, except that it is a second derivative that has been differencedontheright-handside. But this makesaworldofdifference! TheFTCS schemewasunstableforthehyperbolicequation;however,aquickcalculationshows that the amplification factor for equation (19.2.4) is ξ=1−4D∆t (∆x)2sin2/parenleftbiggk∆x 2/parenrightbigg (19.2.5 ) The requirement |ξ|≤1leads to the stability criterion 2D∆t (∆x)2≤1( 19.2.6 ) 19.2DiffusiveInitialValueProblems 839Sample 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 physical interpretation of the restriction (19.2.6) is that the maximum allowed timestep is, up to a numerical factor, the diffusion time across a cell ofwidth ∆x. Moregenerally,the diffusiontime τacross a spatial scale ofsize λis oforder τ∼λ 2 D(19.2.7 ) Usually we are interested in modeling accurately the evolution of features with spatial scales λ/greatermuch∆x. If we are limited to timesteps satisfying (19.2.6), we will needtoevolvethroughoforder λ2/(∆x)2steps beforethingsstartto happenonthe scale of interest. This number of steps is usually prohibitive. We must therefore find a stable way of taking timesteps comparable to, or perhaps — for accuracy — somewhat smaller than, the time scale of (19.2.7). This goal poses an immediate “philosophical” question. Obviously the large timesteps that we propose to take are going to be woefully inaccurate for the small scales that we have decided not to be interested in. We want those scales to dosomething stable, “innocuous,” and perhaps not too physically unreasonable. We want to build this innocuous behavior into our differencing scheme. What should it be? There are two different answers, each of which has its pros and cons. The first answer is to seek a differencingscheme that drives small-scale features to theirequilibrium forms, e.g., satisfying equation (19.2.3) with the left-hand side set to zero. Thisanswergenerallymakesthebestphysicalsense;but,aswewillsee,itleads toadifferencingscheme(“fullyimplicit”)thatisonly first-order accurateintimefor the scales that we are interested in. The second answer is to let small-scale features maintain their initial amplitudes, so that the evolution of the larger-scale features of interest takes place superposed with a kind of “frozen in” (though fluctuating) backgroundof small-scale stuff. This answer gives a differencingscheme (“Crank- Nicolson”) that is second-order accurate in time. Toward the end of an evolution calculation,however,onemightwanttoswitch overtosomesteps oftheotherkind, to drive the small-scale stuff into equilibrium. Let us now see where these distinct differencing schemes come from: Consider the following differencing of (19.2.3), u n+1 j−un j ∆t=D/bracketleftBigg un+1 j+1−2un+1 j+un+1 j−1 (∆x)2/bracketrightBigg (19.2.8 ) This is exactly like the FTCS scheme (19.2.4),except that the spatial derivatives on the right-handside are evaluatedat timestep n+1. Schemes with this characterare calledfully implicit orbackward time , by contrast with FTCS (which is called fully explicit). To solve equation (19.2.8) one has to solve a set of simultaneous linear equationsateachtimestepforthe un+1 j. Fortunately,thisisasimpleproblembecause the system is tridiagonal: Just groupthe terms in equation(19.2.8)appropriately: −αun+1 j−1+( 1+2 α)un+1 j−αun+1 j+1=un j,j =1,2...J−1(19.2.9 ) where α≡D∆t (∆x)2(19.2.10 ) 840 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).Supplemented by Dirichlet or Neumann boundary conditions at j=0andj=J, equation(19.2.9)is clearly a tridiagonalsystem, which can easily be solved at eachtimestep by the method of §2.4. What is the behavior of (19.2.8) for very large timesteps? The answer is seen mostclearlyin(19.2.9),inthelimit α→∞(∆t→∞). Dividingby α, wesee that thedifferenceequationsarejustthefinite-differenceformoftheequilibriumequation ∂ 2u ∂x2=0 ( 19.2.11 ) What about stability? The amplification factor for equation (19.2.8)is ξ=1 1+4 αsin2/parenleftbiggk∆x 2/parenrightbigg (19.2.12 ) Clearly |ξ|<1foranystepsize ∆t. Theschemeisunconditionallystable. Thedetails of the small-scale evolution from the initial conditions are obviously inaccurate for large ∆t. But, as advertised, the correct equilibrium solution is obtained. This is the characteristic feature of implicit methods. Here,ontheotherhand,ishowonegetstothesecondofourabovephilosophical answers,combiningthestabilityofanimplicitmethodwiththeaccuracyofamethodthat is second-orderin bothspace and time. Simply formthe averageof the explicit and implicit FTCS schemes: u n+1 j−un j ∆t=D 2/bracketleftBigg (un+1 j+1−2un+1 j+un+1 j−1)+(un j+1−2un j+un j−1) (∆x)2/bracketrightBigg (19.2.13 ) Hereboththeleft-andright-handsidesarecenteredattimestep n+1 2,sothemethod is second-orderaccurate in time as claimed. The amplification factor is ξ=1−2αsin2/parenleftbiggk∆x 2/parenrightbigg 1+2 αsin2/parenleftbiggk∆x 2/parenrightbigg (19.2.14 ) so the method is stable for any size ∆t. This scheme is called the Crank-Nicolson scheme,andisourrecommendedmethodforanysimplediffusionproblem(perhaps supplementedby a few fully implicit steps at the end). (See Figure 19.2.1.) Now turn to some generalizations of the simple diffusion equation (19.2.3). Supposefirstthatthediffusioncoefficient Disnotconstant,say D=D(x). Wecan adopteither of two strategies. First, we can make an analyticchangeof variable y=/integraldisplaydx D(x)(19.2.15 ) 19.2DiffusiveInitialValueProblems 841Sample 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).t or n x or jFTCS (a) Fully Implicit (b) Crank-Nicolson (c) Figure 19.2.1. Three differencing schemes for diffusive problems (shown as in Figure 19.1.2). (a) ForwardTimeCenterSpaceis first-orderaccurate,butstableonlyforsuf ficientlysmalltimesteps. (b)Fully Implicit is stable for arbitrarily large timesteps, but is still only first-order accurate. (c) Crank-Nicolson is second-order accurate, and is usually stable for large timesteps. Then ∂u ∂t=∂ ∂xD(x)∂u ∂x(19.2.16 ) becomes ∂u ∂t=1 D(y)∂2u ∂y2(19.2.17 ) andweevaluate Dattheappropriate yj. Heuristically,thestabilitycriterion(19.2.6) in an explicit scheme becomes ∆t≤min j/bracketleftBigg (∆y)2 2D−1 j/bracketrightBigg (19.2.18 ) Note that constant spacing ∆yinydoes not imply constant spacing in x. An alternative method that does not require analytically tractable forms for Dis simply to difference equation (19.2.16) as it stands, centering everything appropriately. Thus the FTCS method becomes un+1 j−un j ∆t=Dj+1 /2(un j+1−un j)−Dj−1/2(un j−un j−1) (∆x)2(19.2.19 ) where Dj+1 /2≡D(xj+1 /2)( 19.2.20 ) 842 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).and the heuristic stability criterion is ∆t≤min j/bracketleftbigg(∆x)2 2Dj+1 /2/bracketrightbigg (19.2.21 ) The Crank-Nicolson method can be generalized similarly. The second complication one can consider is a nonlinear diffusion problem, for example where D=D(u). Explicit schemes can be generalized in the obvious way. For example, in equation (19.2.19) write Dj+1 /2=1 2/bracketleftbig D(un j+1)+D(un j)/bracketrightbig (19.2.22 ) Implicitschemesare notas easy. Thereplacement(19.2.22)with n→n+1leaves us with a nasty set of coupled nonlinear equations to solve at each timestep. Often there is an easier way: If the form of D(u)allows us to integrate dz=D(u)du (19.2.23 ) analytically for z(u), then the right-handside of (19.2.1)becomes ∂2z/∂x2, which we difference implicitly as zn+1 j+1−2zn+1 j+zn+1 j−1 (∆x)2(19.2.24 ) Now linearizeeachterm onthe right-handside of equation(19.2.24),forexample zn+1 j≡z(un+1 j)=z(un j)+(un+1 j−un j)∂z ∂u/vextendsingle/vextendsingle/vextendsingle/vextendsingle j,n =z(un j)+(un+1 j−un j)D(un j)(19.2.25 ) This reduces the problem to tridiagonal form again and in practice usually retains the stability advantages of fully implicit differencing. Schr¨odingerEquation Sometimes the physical problem being solved imposes constraints on the differencingscheme that we havenot yet takeninto account. For example,consider thetime-dependentSchr ¨odingerequationofquantummechanics. Thisis basicallya parabolicequationforthe evolutionofa complexquantity ψ. Forthescatteringofa wavepacket by a one-dimensionalpotential V(x), the equation has the form i∂ψ ∂t=−∂2ψ ∂x2+V(x)ψ (19.2.26 ) (Here we have chosen units so that Planck ’s constant ¯h=1and the particle mass m=1/2.) Oneis giventheinitialwavepacket, ψ(x, t=0 ),togetherwithboundary 19.2DiffusiveInitialValueProblems 843Sample 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).conditions that ψ→0atx→± ∞. Suppose we content ourselves with first- order accuracy in time, but want to use an implicit scheme, for stability. A slightgeneralization of (19.2.8) leads to i/bracketleftBigg ψ n+1 j−ψn j ∆t/bracketrightBigg =−/bracketleftBigg ψn+1 j+1−2ψn+1 j+ψn+1 j−1 (∆x)2/bracketrightBigg +Vjψn+1 j (19.2.27 ) for which ξ=1 1+i/bracketleftbigg4∆t (∆x)2sin2/parenleftbiggk∆x 2/parenrightbigg +Vj∆t/bracketrightbigg (19.2.28 ) This is unconditionallystable, butunfortunatelyis not unitary. Theunderlying physicalproblemrequiresthatthetotalprobabilityof findingtheparticlesomewhere remains unity. This is represented formally by the modulus-square norm of ψ remaining unity: /integraldisplay∞ −∞|ψ|2dx=1 ( 19.2.29 ) Theinitialwavefunction ψ(x,0)isnormalizedtosatisfy(19.2.29). TheSchr ¨odinger equation(19.2.26)thenguaranteesthat this conditionis satis fied at all later times. Let us write equation (19.2.26) in the form i∂ψ ∂t=Hψ (19.2.30 ) where the operator His H=−∂2 ∂x2+V(x)( 19.2.31 ) The formal solution of equation (19.2.30) is ψ(x, t)=e−iHtψ(x,0) ( 19.2.32 ) where the exponentialof the operatoris de fined by its power series expansion. The unstable explicit FTCS scheme approximates (19.2.32) as ψn+1 j=( 1−iH∆t)ψn j (19.2.33 ) where His represented by a centered finite-difference approximation in x. The stable implicit scheme (19.2.27) is, by contrast, ψn+1 j=( 1+ iH∆t)−1ψn j (19.2.34 ) These are both first-order accurate in time, as can be seen by expanding equation (19.2.32). However, neither operator in (19.2.33)or (19.2.34)is unitary. 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 itsfinite-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