f19-1

19.1Flux-ConservativeInitialValueProblems 825Sample 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).engineering; these methods allow considerable freedom in putting computational elementswhereyouwantthem,importantwhendealingwithhighlyirregulargeome-tries. Spectral methods [13-15]are preferredfor veryregulargeometriesand smooth functions;theyconvergemorerapidlythanfinite-differencemethods(cf. §19.4),but they do not work well for problems with discontinuities. CITED REFERENCES AND FURTHER READING: Ames, W.F. 1977, Numerical Methods for Partial Differential Equations , 2nd ed. (New York: Academic Press). [1] Richtmyer, R.D., andMorton, K.W.1967, DifferenceMethods for InitialValue Problems ,2nded. (New York: Wiley-Interscience). [2] Roache, P.J. 1976, Computational Fluid Dynamics (Albuquerque: Hermosa). [3] Mitchell,A.R., andGriffiths, D.F.1980, TheFiniteDifferenceMethod inPartialDifferentialEqua- tions(New York: Wiley) [includes discussion of finite element methods]. [4] Dorr, F.W. 1970, SIAM Review , vol. 12, pp. 248–263. [5] Meijerink, J.A., and van der Vorst, H.A. 1977, Mathematics of Computation , vol. 31, pp. 148– 162. [6] van der Vorst, H.A. 1981, Journalof Computational Physics , vol. 44, pp. 1–19[review of sparse iterative methods]. [7] Kershaw, D.S. 1970, Journal of Computational Physics , vol. 26, pp. 43–65. [8] Stone, H.J. 1968, SIAM Journal on Numerical Analysis , vol. 5, pp. 530–558. [9] Jesshope, C.R. 1979, Computer Physics Communications , vol. 17, pp. 383–391. [10] Strang, G., and Fix, G. 1973, An Analysis of the Finite Element Method (Englewood Cliffs, NJ: Prentice-Hall). [11] Burnett, D.S. 1987, Finite Element Analysis: From Concepts to Applications (Reading, MA: Addison-Wesley). [12] Gottlieb, D. andOrszag, S.A. 1977, NumericalAnalysis of Spectral Methods: Theoryand Appli- cations(Philadelphia: S.I.A.M.). [13] Canuto, C., Hussaini, M.Y., Quarteroni, A., and Zang, T.A. 1988, Spectral Methods in Fluid Dynamics (New York: Springer-Verlag). [14] Boyd, J.P. 1989, Chebyshev and FourierSpectral Methods (New York: Springer-Verlag). [15] 19.1 Flux-Conservative InitialValue Problems Alargeclassofinitialvalue(time-evolution)PDEsinonespacedimensioncan be cast into the form of a flux-conservative equation , ∂u ∂t=−∂F(u) ∂x(19.1.1 ) whereuandFare vectors, and where (in some cases) Fmay depend not only on u but also on spatial derivativesof u. The vector Fis called the conserved flux . For example, the prototypical hyperbolic equation, the one-dimensional wave equation with constant velocity of propagation v ∂2u ∂t2=v2∂2u ∂x2(19.1.2 ) 826 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).can be rewritten as a set of two first-order equations ∂r ∂t=v∂s ∂x ∂s ∂t=v∂r ∂x(19.1.3 ) where r≡v∂u ∂x s≡∂u ∂t(19.1.4 ) In this case randsbecome the two components of u, and the flux is given by the linear matrix relation F(u)=/parenleftbigg 0−v −v0/parenrightbigg ·u (19.1.5 ) (The physicist-reader may recognize equations (19.1.3) as analogous to Maxwell’s equations for one-dimensional propagation of electromagnetic waves.) We will consider, in this section, a prototypical example of the general flux- conservative equation (19.1.1), namely the equation for a scalar u, ∂u ∂t=−v∂u ∂x(19.1.6 ) with va constant. As it happens, we already know analytically that the general solution of this equation is a wave propagatingin the positive x-direction, u=f(x−vt)( 19.1.7 ) where fis an arbitraryfunction. However,the numericalstrategies that we develop will be equally applicableto the moregeneral equationsrepresentedby (19.1.1). In somecontexts,equation(19.1.6)iscalledan advective equation,becausethequantity uis transported by a “fluid flow” with a velocity v. How do we go about finite differencing equation (19.1.6) (or, analogously, 19.1.1)? Thestraightforwardapproachistochooseequallyspacedpointsalongboth thet- and x-axes. Thus denote xj=x0+j∆x, j =0,1,...,J tn=t0+n∆t, n =0,1,...,N(19.1.8 ) Letun jdenote u(tn,x j). We have several choices for representing the time derivative term. The obvious way is to set ∂u ∂t/vextendsingle/vextendsingle/vextendsingle/vextendsingle j,n=un+1 j−un j ∆t+O(∆t)( 19.1.9 ) Thisis called forwardEuler differencing(cf.equation16.1.1). WhileforwardEuler is only first-order accurate in ∆t, it has the advantage that one is able to calculate 19.1Flux-ConservativeInitialValueProblems 827Sample 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 Figure 19.1.1. Representation ofthe Forward TimeCentered Space (FTCS)differencing scheme. In this and subsequent figures, the open circle is the new point at which the solution is desired; filled circles are known points whose function values are used in calculating the new point; the solid lines connect pointsthat areusedtocalculate spatial derivatives; thedashed linesconnect pointsthatareusedtocalculate timederivatives. TheFTCSschemeisgenerally unstable forhyperbolic problemsandcannot usuallybeused. quantitiesat timestep n+1intermsofonlyquantitiesknownat timestep n. Forthe spacederivative,wecanuseasecond-orderrepresentationstill usingonlyquantities known at timestep n: ∂u ∂x/vextendsingle/vextendsingle/vextendsingle/vextendsingle j,n=un j+1−un j−1 2∆x+O(∆x2)( 19.1.10 ) Theresulting finite-differenceapproximationtoequation(19.1.6)iscalledtheFTCS representation (Forward Time Centered Space), un+1 j−un j ∆t=−v/parenleftbiggun j+1−un j−1 2∆x/parenrightbigg (19.1.11 ) which can easily be rearranged to be a formula for un+1 jin terms of the other quantities. The FTCS scheme is illustrated in Figure 19.1.1. It ’safine example of an algorithm that is easy to derive, takes little storage, and executes quickly. Too bad it doesn ’t work! (See below.) The FTCS representationis an explicitscheme. This means that un+1 jforeach jcan be calculated explicitly from the quantities that are already known. Later we shall meet implicitschemes, which require us to solve implicit equations coupling theun+1 jfor various j. (Explicit and implicit methods for ordinary differential equations were discussed in §16.6.) The FTCS algorithm is also an example of asingle-level scheme, since only values at time level nhave to be stored to find values at time level n+1. von NeumannStability Analysis Unfortunately,equation(19.1.11)isofverylimitedusefulness. Itisan unstable method,which can be used only (if at all) to study waves for a short fractionof one oscillation period. To find alternative methods with more general applicability, we must introduce the von Neumann stability analysis . The von Neumann analysis is local: We imagine that the coef ficients of the difference equations are so slowly varying as to be considered constant in spaceand time. In that case, the independent solutions, or eigenmodes , of the difference equations are all of the form u n j=ξneikj∆x(19.1.12 ) 828 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).t or n x or jLax Figure 19.1.2. Representation of the Lax differencing scheme, as in the previous figure. The stability criterion for this scheme is the Courant condition. where kis a real spatial wave number (which can have any value) and ξ=ξ(k)is a complex number that depends on k. The key fact is that the time dependence of a single eigenmode is nothing more than successive integer powers of the complex number ξ. Therefore, the difference equations are unstable (have exponentially growing modes) if |ξ(k)|>1forsome k. The number ξis called the amplification factorat a given wave number k. Tofindξ(k), we simply substitute (19.1.12) back into (19.1.11). Dividing byξn, we get ξ(k)=1−iv∆t ∆xsink∆x (19.1.13 ) whose modulusis >1forallk; so the FTCS scheme is unconditionallyunstable. Ifthevelocity vwereafunctionof tandx,thenwewouldwrite vn jinequation (19.1.11). InthevonNeumannstabilityanalysiswewouldstill treat vasaconstant, the idea being that for vslowly varying the analysis is local. In fact, even in the case of strictly constant v, the von Neumann analysis does not rigorously treat the end effects at j=0andj=N. More generally, if the equation ’s right-hand side were nonlinear in u, then a vonNeumannanalysiswouldlinearizebywriting u=u0+δu,expandingtolinear order in δu. Assuming that the u0quantities already satisfy the difference equation exactly, the analysis would look for an unstable eigenmode of δu. Despiteitslackofrigor,thevonNeumannmethodgenerallygivesvalidanswers and is much easier to apply than more careful methods. We accordingly adopt it exclusively. (See, for example, [1]for a discussion of other methods of stability analysis.) Lax Method TheinstabilityintheFTCSmethodcanbecuredbyasimplechangeduetoLax. Onereplacesthe term un jinthe timederivativetermbyits average(Figure19.1.2): un j→1 2/parenleftbig un j+1+un j−1/parenrightbig (19.1.14 ) This turns (19.1.11) into un+1 j=1 2/parenleftbig un j+1+un j−1/parenrightbig −v∆t 2∆x/parenleftbig un j+1−un j−1/parenrightbig (19.1.15 ) 19.1Flux-ConservativeInitialValueProblems 829Sample 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∆t x or j∆t ∆x ∆xunstable stable (a) (b) Figure 19.1.3. Courant condition for stability of a differencing scheme. The solution of a hyperbolic problem at a point depends on information within some domain of dependency to the past, shown hereshaded. The differencing scheme (19.1.15) has its own domain of dependency determined by the choice of points on one time slice (shown as connected solid dots) whose values are used in determining a new point (shown connected by dashed lines). A differencing scheme is Courant stable if the differencingdomain of dependency is larger than that of the PDEs, as in (a), and unstable if the relationship is thereverse, as in (b). For more complicated differencing schemes, the domain of dependency might not be determined simply by the outermost points. Substituting equation (19.1.12), we find for the ampli fication factor ξ=c o s k∆x−iv∆t ∆xsink∆x (19.1.16 ) The stability condition |ξ|2≤1leads to the requirement |v|∆t ∆x≤1( 19.1.17 ) This is the famous Courant-Friedrichs-Lewy stability criterion, often called simply the Courant condition . Intuitively, the stability condition can be understood as follows (Figure 19.1.3): The quantity un+1 jin equation (19.1.15) is computed from information at points j−1andj+1at time n. In other words, xj−1andxj+1aretheboundariesofthespatialregionthatisallowedtocommunicate information to un+1 j. Now recall that in the continuumwave equation, information actually propagates with a maximum velocity v. If the point un+1 jis outside of the shaded region in Figure 19.1.3, then it requires information from points moredistant than the differencing scheme allows. Lack of that information gives rise to an instability. Therefore, ∆tcannot be made too large. Thesurprisingresult,thatthesimplereplacement(19.1.14)stabilizestheFTCS scheme, is our first encounterwith the fact that differencingPDEs is an art as much as a science. Tosee if we candemystifytheart somewhat,let us comparetheFTCS andLaxschemesbyrewritingequation(19.1.15)sothatitis intheformofequation (19.1.11) with a remainder term: u n+1 j−un j ∆t=−v/parenleftbiggun j+1−un j−1 2∆x/parenrightbigg +1 2/parenleftbiggun j+1−2un j+un j−1 ∆t/parenrightbigg (19.1.18 ) But this is exactly the FTCS representation of the equation ∂u ∂t=−v∂u ∂x+(∆x)2 2∆t∇2u (19.1.19 ) 830 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).where∇2=∂2/∂x2inonedimension. Wehave,ineffect,addedadiffusiontermto theequation,or,ifyourecalltheformoftheNavier-Stokesequationforviscous fluid flow, a dissipative term. The Lax schemeis thus said to have numericaldissipation , ornumericalviscosity . Wecanseethisalsointheampli ficationfactor. Unless |v|∆t is exactlyequalto ∆x,|ξ|<1andtheamplitudeofthewave decreasesspuriously. Isn’ta spuriousdecreaseas badas aspuriousincrease? No. Thescales that we hopetostudyaccuratelyarethosethatencompassmanygridpoints,sothattheyhave k∆x/lessmuch1. (The spatial wave number kis defined by equation 19.1.12.) For these scales,theampli ficationfactorcanbeseentobeveryclosetoone,inboththestable and unstable schemes. The stable and unstable schemes are thereforeabout equallyaccurate. For the unstable scheme, however, short scales with k∆x∼1,which we are not interested in , will blow up and swamp the interesting part of the solution. Much better to have a stable scheme in which these short wavelengths die awayinnocuously. Boththestableandtheunstableschemesare inaccurate fortheseshort wavelengths,buttheinaccuracyisofatolerablecharacterwhentheschemeisstable. When the independent variable uis a vector, then the von Neumann analysis is slightly more complicated. For example, we can consider equation (19.1.3), rewritten as ∂ ∂t/bracketleftbigg r s/bracketrightbigg =∂ ∂x/bracketleftbigg vs vr/bracketrightbigg (19.1.20 ) The Lax method for this equation is rn+1 j=1 2(rn j+1+rn j−1)+v∆t 2∆x(sn j+1−sn j−1) sn+1 j=1 2(sn j+1+sn j−1)+v∆t 2∆x(rn j+1−rn j−1)(19.1.21 ) The von Neumann stability analysis now proceeds by assuming that the eigenmode is of the following (vector) form, /bracketleftbigg rn j sn j/bracketrightbigg =ξneikj∆x/bracketleftbigg r0 s0/bracketrightbigg (19.1.22 ) Here the vector on the right-hand side is a constant (both in space and in time) eigenvector, and ξis a complex number, as before. Substituting (19.1.22) into (19.1.21),anddividingbythe power ξn, givesthe homogeneousvectorequation  (cosk∆x)−ξiv∆t ∆xsink∆x iv∆t ∆xsink∆x (cosk∆x)−ξ · r0 s0 = 0 0  (19.1.23 ) This admits a solution only if the determinant of the matrix on the left vanishes, a condition easily shown to yield the two roots ξ ξ=c o s k∆x±iv∆t ∆xsink∆x (19.1.24 ) The stability condition is that both roots satisfy |ξ|≤1. This again turns out to be simply the Courant condition (19.1.17). 19.1Flux-ConservativeInitialValueProblems 831Sample 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).OtherVarieties of Error Thus far we have beenconcernedwith amplitudeerror , because of its intimate connectionwith the stability or instability of a differencingscheme. Other varieties of errorare relevantwhen we shift our concernto accuracy,ratherthan stability. Finite-difference schemes for hyperbolic equations can exhibit dispersion, or phase errors . For example, equation (19.1.16) can be rewritten as ξ=e−ik∆x+i/parenleftbigg 1−v∆t ∆x/parenrightbigg sink∆x (19.1.25 ) An arbitrary initial wave packet is a superposition of modes with different k’s. At each timestep the modes get multiplied by different phase factors (19.1.25), dependingontheirvalueof k.I f∆t=∆x/v,thentheexactsolutionforeachmode ofawavepacket f(x−vt)isobtainedifeachmodegetsmultipliedby exp(−ik∆x). For this value of ∆t, equation (19.1.25) shows that the finite-difference solution gives the exact analytic result. However, if v∆t/∆xis not exactly 1, the phase relationsofthemodescanbecomehopelesslygarbledandthewavepacketdisperses. Note from (19.1.25) that the dispersion becomes large as soon as the wavelength becomes comparable to the grid spacing ∆x. A third typeof erroris one associated with nonlinearhyperbolicequationsand isthereforesometimescalled nonlinearinstability . Forexample,apieceoftheEuler or Navier-Stokes equations for fluidflow looks like ∂v ∂t=−v∂v ∂x+... (19.1.26 ) The nonlinear term in vcan cause a transfer of energy in Fourier space from long wavelengths to short wavelengths. This results in a wave pro file steepening until a vertical pro file or“shock”develops. Since the von Neumann analysis suggests that the stability can dependon k∆x, a scheme that was stable forshallow profiles can become unstable for steep pro files. This kind of dif ficulty arises in a differencing scheme where the cascade in Fourier space is halted at the shortest wavelength representable on the grid, that is, at k∼1/∆x. If energy simply accumulatesinthesemodes,iteventuallyswampstheenergyinthelongwavelength modes of interest. Nonlinear instability and shock formation is thus somewhat controlled by numerical viscosity such as that discussed in connection with equation (19.1.18) above. Insome fluidproblems,however,shockformationisnotmerelyanannoyance, but an actual physical behavior of the fluid whose detailed study is a goal. Then, numericalviscosityalonemaynotbe adequateorsuf ficientlycontrollable. Thisis a complicated subject which we discuss further in the subsection on fluid dynamics, below. For wave equations, propagation errors (amplitude or phase) are usually most worrisome. For advectiveequations, on the other hand, transport errors are usually of greater concern. In the Lax scheme, equation (19.1.15), a disturbance in the advected quantity uat mesh point jpropagates to mesh points j+1andj−1at the next timestep. In reality, however, if the velocity vis positive then only mesh point j+1should be affected. 832 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).t or n x or jv upwind v Figure 19.1.4. Representation of upwind differencing schemes. The upper scheme is stable when the advection constant vis negative, as shown; the lower scheme is stable when the advection constant vis positive, also as shown. The Courant condition must, of course, also be satis fied. The simplest way to model the transport properties “better”is to use upwind differencing (see Figure 19.1.4): un+1 j−un j ∆t=−vn j  un j−un j−1 ∆x,vn j>0 un j+1−un j ∆x,vn j<0(19.1.27 ) Note that this scheme is only first-order, not second-order, accurate in the calculation of the spatial derivatives. How can it be “better”? The answer is one that annoys the mathematicians: The goal of numerical simulations is not always“accuracy”in a strictly mathematical sense, but sometimes “fidelity”to the underlying physics in a sense that is looser and more pragmatic. In such contexts, some kinds of error are much more tolerable than others. Upwind differencing generallyadds fidelitytoproblemswheretheadvectedvariablesareliabletoundergo sudden changes of state, e.g., as they pass through shocks or other discontinuities. You will have to be guided by the speci fic nature of your own problem. Forthedifferencingscheme(19.1.27),theampli ficationfactor(forconstant v)is ξ=1−/vextendsingle/vextendsingle/vextendsingle/vextendsinglev∆t ∆x/vextendsingle/vextendsingle/vextendsingle/vextendsingle(1−cosk∆x)−iv∆t ∆xsink∆x (19.1.28 ) |ξ|2=1−2/vextendsingle/vextendsingle/vextendsingle/vextendsinglev∆t ∆x/vextendsingle/vextendsingle/vextendsingle/vextendsingle/parenleftbigg 1−/vextendsingle/vextendsingle/vextendsingle/vextendsinglev∆t ∆x/vextendsingle/vextendsingle/vextendsingle/vextendsingle/parenrightbigg (1−cosk∆x)( 19.1.29 ) So thestability criterion |ξ|2≤1is (again)simplytheCourantcondition(19.1.17). There are various ways of improving the accuracy of first-order upwind differencing. In the continuum equation, material originally a distance v∆taway 19.1Flux-ConservativeInitialValueProblems 833Sample 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).staggered leapfrogt or n x or j Figure 19.1.5. Representation of the staggered leapfrog differencing scheme. Note that information from two previous time slices is used in obtaining the desired point. This scheme is second-order accurate in both space and time. arrives at a given point after a time interval ∆t. In the first-order method, the material always arrives from ∆xaway. If v∆t/lessmuch∆x(to insure accuracy),this can cause a largeerror. One wayof reducingthis erroris tointerpolate ubetween j−1 andjbefore transporting it. This gives effectively a second-ordermethod. Various schemesforsecond-orderupwinddifferencingare discussedandcomparedin [2-3]. Second-Order Accuracy inTime When using a method that is first-order accurate in time but second-order accurate in space, one generally has to take v∆tsignificantly smaller than ∆xto achieve desired accuracy, say, by at least a factor of 5. Thus the Courant condition is not actually the limiting factorwith such schemes in practice. However,there are schemesthataresecond-orderaccurateinbothspaceandtime,andthesecanoftenbe pushedrighttotheirstabilitylimit,withcorrespondinglysmallercomputationtimes. For example, the staggered leapfrog method for the conservation equation (19.1.1) is de fined as follows (Figure 19.1.5): Using the values of unat time tn, compute the fluxes Fn j. Then compute new values un+1using the time-centered values of the fluxes: un+1 j−un−1 j=−∆t ∆x(Fn j+1−Fn j−1)( 19.1.30 ) The name comes from the fact that the time levels in the time derivative term “leapfrog”over the time levels in the space derivative term. The method requires thatun−1andunbe stored to compute un+1. For our simple modelequation(19.1.6),staggered leapfrogtakes the form un+1 j−un−1 j=−v∆t ∆x(un j+1−un j−1)( 19.1.31 ) ThevonNeumannstabilityanalysisnowgivesaquadraticequationfor ξ,ratherthan a linear one, because of the occurrence of three consecutive powers of ξwhen the 834 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).form (19.1.12) for an eigenmode is substituted into equation (19.1.31), ξ2−1=−2iξv∆t ∆xsink∆x (19.1.32 ) whose solution is ξ=−iv∆t ∆xsink∆x±/radicalBigg 1−/parenleftbiggv∆t ∆xsink∆x/parenrightbigg2 (19.1.33 ) Thus the Courant condition is again required for stability. In fact, in equation (19.1.33), |ξ|2=1for any v∆t≤∆x. This is the great advantageof the staggered leapfrog method: There is no amplitude dissipation. Staggered leapfrog differencingof equations like (19.1.20)is most transparent if the variables are centered on appropriate half-mesh points: rn j+1 /2≡v∂u ∂x/vextendsingle/vextendsingle/vextendsingle/vextendsinglen j+1 /2=vun j+1−un j ∆x sn+1 /2 j≡∂u ∂t/vextendsingle/vextendsingle/vextendsingle/vextendsinglen+1 /2 j=un+1 j−un j ∆t(19.1.34 ) This is purely a notational convenience: we can think of the mesh on which rand sare defined as being twice as fine as the mesh on which the original variable uis defined. The leapfrog differencing of equation (19.1.20) is rn+1 j+1 /2−rn j+1 /2 ∆t=sn+1 /2 j+1−sn+1 /2 j ∆x sn+1 /2 j−sn−1/2 j ∆t=vrn j+1 /2−rn j−1/2 ∆x(19.1.35 ) If you substitute equation (19.1.22) in equation (19.1.35), you will find that once again the Courant condition is required for stability, and that there is no amplitude dissipation when it is satis fied. If we substitute equation(19.1.34)in equation(19.1.35),we find that equation (19.1.35) is equivalent to un+1 j−2un j+un−1 j (∆t)2=v2un j+1−2un j+un j−1 (∆x)2(19.1.36 ) Thisis justthe “usual”second-orderdifferencingofthewaveequation(19.1.2). We see that it is a two-level scheme, requiring both unandun−1to obtain un+1.I n equation (19.1.35) this shows up as both sn−1/2andrnbeing needed to advance the solution. For equations more complicated than our simple model equation, especially nonlinearequations, the leapfrogmethod usually becomesunstable when the gradi- entsgetlarge. Theinstabilityis relatedtothefactthatoddandevenmeshpointsare completely decoupled, like the black and white squares of a chess board, as shown 19.1Flux-ConservativeInitialValueProblems 835Sample 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).Figure 19.1.6. Origin of mesh-drift instabilities in a staggered leapfrog scheme. If the mesh points are imagined to lie in the squares of a chess board, then white squares couple to themselves, black tothemselves, but there is no coupling between white and black. The fix is to introduce a small diffusive mesh-coupling piece. in Figure 19.1.6. This mesh driftinginstability is curedbycouplingthe two meshes throughanumericalviscosityterm,e.g.,addingtotherightsideof(19.1.31)asmall coefficient (/lessmuch1) times un j+1−2un j+un j−1. For more on stabilizing difference schemes by adding numerical dissipation, see, e.g., [4]. TheTwo-Step Lax-Wendroff scheme is a second-order in time method that avoids large numerical dissipation and mesh drifting. One de fines intermediate values uj+1 /2at the half timesteps tn+1 /2and the half mesh points xj+1 /2. These are calculated by the Lax scheme: un+1 /2 j+1 /2=1 2(un j+1+un j)−∆t 2∆x(Fn j+1−Fn j)( 19.1.37 ) Using these variables, one calculates the fluxes Fn+1 /2 j+1 /2. Then the updated values un+1 jare calculated by the properly centered expression un+1 j=un j−∆t ∆x/parenleftBig Fn+1 /2 j+1 /2−Fn+1 /2 j−1/2/parenrightBig (19.1.38 ) The provisional values un+1 /2 j+1 /2are now discarded. (See Figure 19.1.7.) Letusinvestigatethestabilityofthismethodforourmodeladvectiveequation, where F=vu. Substitute (19.1.37) in (19.1.38) to get un+1 j=un j−α/bracketleftbigg1 2(un j+1+un j)−1 2α(un j+1−un j) −1 2(un j+un j−1)+1 2α(un j−un j−1)/bracketrightbigg (19.1.39 ) 836 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).t or n x or jhalfstep pointstwo-step Lax Wendroff Figure 19.1.7. Representation of the two-step Lax-Wendroff differencing scheme. Two halfstep points (⊗) are calculated by the Lax method. These, plus one of the original points, produce the new point via staggered leapfrog. Halfstep points are usedonly temporarily anddonot require storage allocation onthegrid. This scheme is second-order accurate in both space and time. where α≡v∆t ∆x(19.1.40 ) Then ξ=1−iαsink∆x−α2(1−cosk∆x)( 19.1.41 ) so |ξ|2=1−α2(1−α2)(1−cosk∆x)2(19.1.42 ) The stability criterion |ξ|2≤1is therefore α2≤1,orv∆t≤∆xas usual. Incidentally, you should not think that the Courant condition is the only stability requirement that ever turns up in PDEs. It keeps doing so in our model examples just because those examples are so simple in form. The method of analysis is,however, general. Except when α=1,|ξ| 2<1in (19.1.42), so some amplitude damping does occur. The effect is relativelysmall, however,for wavelengths large comparedwiththe mesh size ∆x. If we expand (19.1.42) for small k∆x,w efind |ξ| 2=1−α2(1−α2)(k∆x)4 4+... (19.1.43 ) Thedeparturefromunityoccursonlyat fourthorderin k. Thisshouldbecontrasted with equation (19.1.16) for the Lax method, which shows that |ξ|2=1−(1−α2)(k∆x)2+... (19.1.44 ) for small k∆x. 19.1Flux-ConservativeInitialValueProblems 837Sample 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).In summary,our recommendationfor initial value problemsthat can be cast in flux-conservativeform,andespeciallyproblemsrelatedtothewaveequation,istouse thestaggeredleapfrogmethodwhenpossible. Wehavepersonallyhadbettersuccess with it than with the Two-Step Lax-Wendroff method. For problems sensitive to transporterrors,upwinddifferencingoroneofits re finementsshouldbeconsidered. FluidDynamics withShocks As we alludedto earlier,the treatment of fluid dynamicsproblemswith shocks has become a very complicated and very sophisticated subject. All we can attempt to do here is to guide you to some starting points in the literature. There are basically three important general methods for handling shocks. The oldest and simplest method, invented by von Neumann and Richtmyer, is to add artificial viscosity to the equations, modeling the way Nature uses real viscosity to smooth discontinuities. A good starting point for trying out this method is the differencingschemein §12.11of [1]. Thisschemeisexcellentfornearlyallproblems in one spatial dimension. Thesecondmethodcombinesahigh-orderdifferencingschemethatisaccurate for smooth flows with a low order scheme that is very dissipative and can smooth the shocks. Typically, various upwind differencing schemes are combined usingweights chosento zero the low orderschemeunless steep gradientsare present, and also chosen to enforce various “monotonicity ”constraints that prevent nonphysical oscillations from appearing in the numerical solution. References [2-3,5]are a good place to start with these methods. Thethird,andpotentiallymostpowerfulmethod,is Godunov ’sapproach. Here one gives up the simple linearization inherent in finite differencingbased on Taylor series and includes the nonlinearity of the equations explicitly. There is an analytic solutionfortheevolutionoftwouniformstatesofa fluidseparatedbyadiscontinuity, the Riemann shock problem. Godunov ’s idea was to approximate the fluid by a large numberof cells of uniform states, and piece them togetherusing the Riemann solution. There have been many generalizations of Godunov ’s approach, of which the most powerful is probably the PPM method [6]. Readable reviews of all these methods, discussing the dif ficulties arising when one-dimensionalmethods are generalizedto multidimensions,are givenin [7-9]. CITED REFERENCES AND FURTHER READING: Ames, W.F. 1977, Numerical Methods for Partial Differential Equations , 2nd ed. (New York: Academic Press), Chapter 4. Richtmyer, R.D., andMorton, K.W.1967, DifferenceMethods for InitialValue Problems ,2nded. (New York: Wiley-Interscience). [1] Centrella, J., and Wilson, J.R. 1984, Astrophysical Journal Supplement , vol. 54, pp. 229–249, Appendix B. [2] Hawley, J.F., Smarr, L.L., and Wilson, J.R. 1984, Astrophysical Journal Supplement , vol. 55, pp. 211–246, §2c. [3] Kreiss, H.-O. 1978, NumericalMethods for SolvingTime-DependentProblemsfor Partial Differ- ential Equations (Montreal: University of Montreal Press), pp. 66ff. [4] Harten, A., Lax, P.D., and Van Leer, B. 1983, SIAM Review , vol. 25, pp. 36–61. [5] Woodward,P., andColella,P. 1984, Journalof ComputationalPhysics , vol. 54,pp. 174–201.[6] Roache, P.J. 1976, Computational Fluid Dynamics (Albuquerque: Hermosa). [7] 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 coef ficient. Actually, this equation is a flux-conservative equation of the form considered in the previous section, with F=−D∂u ∂x(19.2.2 ) theflux in the x-direction. We will assume D≥0, otherwise equation (19.2.1)has physicallyunstablesolutions: Asmalldisturbanceevolvestobecomemoreandmoreconcentratedinsteadofdispersing. (Don ’tmakethemistakeoftryingto findastable 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 arti ficial 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 ampli fication 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 )