834 Chapter 19.Partial Differential Equations engineering;these methods allow considerable freedom in putting computational elements where you want them,important when dealing with highly irregular geome- tries.Spectral methods [13-15]are preferred for very regular geometries and smooth functions;they converge more rapidly than finite-difference methods(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.,and Morton,K.W.1967,Difference Methods for Initial Value Problems,2nd ed. (New York:Wiley-Interscience).[2] Roache,P.J.1976,Computational Fluid Dynamics(Albuquerque:Hermosa).[3] Mitchell,A.R.,and Griffiths,D.F.1980,The Finite Difference Method in Partial Differential Equa- tions (New York:Wiley)[includes discussion of finite element methods].[4] Dorr,F.W.1970,S/AM 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] RECIPES van der Vorst,H.A.1981,Journal of Computational Physics,vol.44,pp.1-19 [review of sparse iterative methods].[7] 9 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] IENTIFIC Gottlieb,D.and Orszag.S.A.1977,Numerical Analysis of Spectral Methods:Theory and Appli- cations(Philadelphia:S.I.A.M.).[13] 6 Canuto,C.,Hussaini,M.Y.,Quarteroni,A,and Zang,T.A.1988,Spectra/Methods in Fluid Dynamics (New York:Springer-Verlag).[14] Boyd,J.P.1989,Chebyshev and Fourier Spectral Methods(New York:Springer-Verlag).[15] Numerica 10.621 19.1 Flux-Conservative Initial Value Problems 4312 Recipes A large class of initial value(time-evolution)PDEs in one space dimension can (outside be cast into the form of a flx-conservative equation, North Ou OF(u) dt=- (19.1.1) Ox where u and F are vectors,and where(in some cases)F may depend not only on u but also on spatial derivatives of u.The vector F is called the conserved flux. For example,the prototypical hyperbolic equation,the one-dimensional wave equation with constant velocity of propagation v 03=v202u 82u (19.1.2) 0x2
834 Chapter 19. Partial Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). engineering; these methods allow considerable freedom in putting computational elements where you want them, important when dealing with highly irregular geometries. Spectral methods [13-15] are preferred for very regular geometries and smooth functions; they converge more rapidly than finite-difference methods (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., and Morton, K.W. 1967, Difference Methods for Initial Value Problems, 2nd ed. (New York: Wiley-Interscience). [2] Roache, P.J. 1976, Computational Fluid Dynamics (Albuquerque: Hermosa). [3] Mitchell, A.R., and Griffiths, D.F. 1980, The Finite Difference Method in Partial Differential Equations (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, Journal of 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. and Orszag, S.A. 1977, Numerical Analysis of Spectral Methods: Theory and Applications (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 Fourier Spectral Methods (New York: Springer-Verlag). [15] 19.1 Flux-Conservative Initial Value Problems A large class of initial value (time-evolution) PDEs in one space dimension can be cast into the form of a flux-conservative equation, ∂u ∂t = −∂F(u) ∂x (19.1.1) where u and F are vectors, and where (in some cases) F may depend not only on u but also on spatial derivatives of u. The vector F is 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)
19.1 Flux-Conservative Initial Value Problems 835 can be rewritten as a set of two first-order equations Or Os Ot Ox Os or (19.1.3) Ot where Ou T三U Ox Ou 8三 t (19.1.4) In this case r and s become the two components of u,and the flux is given by the linear matrix relation ICAL F(u) 0 (19.1.5) RECIPES (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, du du Ot (19.1.6) dx es9&超 9 9 with v a constant.As it happens,we already know analytically that the general IENTIFIC( solution of this equation is a wave propagating in the positive -direction, 6 u=f(x-vt) (19.1.7) where f is an arbitrary function.However,the numerical strategies that we develop will be equally applicable to the more general equations represented by (19.1.1).In some contexts,equation (19.1.6)is called an advective equation,because the quantity u is transported by a "fluid flow"with a velocity v. Recipes Numerical 10621 How do we go about finite differencing equation(19.1.6)(or,analogously, 19.1.1)?The straightforward approach is to choose equally spaced points along both 431 the t-and z-axes.Thus denote Recipes x5=x0+j△x,j=0,1,.,J (19.1.8) tn=to+n△t, n=0,1,..,N Let u"denote u(tn,)We have several choices for representing the time derivative term.The obvious way is to set ou Ot g1-+0(△) (19.1.9) li.n △t This is called forward Euler differencing(cf.equation 16.1.1).While forward Euler is only first-order accurate in At,it has the advantage that one is able to calculate
19.1 Flux-Conservative Initial Value Problems 835 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 r and s become the two components of u, and the flux is given by the linear matrix relation F(u) = 0 −v −v 0 · 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 fluxconservative equation (19.1.1), namely the equation for a scalar u, ∂u ∂t = −v ∂u ∂x (19.1.6) with v a constant. As it happens, we already know analytically that the general solution of this equation is a wave propagating in the positive x-direction, u = f(x − vt) (19.1.7) where f is an arbitrary function. However, the numerical strategies that we develop will be equally applicable to the more general equations represented by (19.1.1). In some contexts, equation (19.1.6) is called an advective equation, because the quantity u is 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)? The straightforward approach is to choose equally spaced points along both the t- and x-axes. Thus denote xj = x0 + j∆x, j = 0, 1,...,J tn = t0 + n∆t, n = 0, 1,...,N (19.1.8) Let un j denote u(tn, xj ). We have several choices for representing the time derivative term. The obvious way is to set ∂u ∂t j,n = un+1 j − un j ∆t + O(∆t) (19.1.9) This is called forward Euler differencing (cf. equation 16.1.1). While forward Euler is only first-order accurate in ∆t, it has the advantage that one is able to calculate
836 Chapter 19.Partial Differential Equations 个 0 FTCS x orj Figure 19.1.1.Representation of the Forward Time Centered 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 points that are used to calculate spatial derivatives;the dashed lines connect points that are used to calculate time derivatives.The FTCS scheme is generally unstable for hyperbolic problems and cannot usually be used. nted for quantities at timestep n+1 in terms of only quantities known at timestep n.For the space derivative,we can use a second-order representation still using only quantities known at timestep n: Ou 0加 u吲+1-1+0（△x2)） RECIPES (19.1.10) 2△x 令 The resulting finite-difference approximation to equation(19.1.6)is called the FTCS representation (Forward Time Centered Space), ,为 u+1-吗 /u+1-u-1 At (19.1.11) 2△x which can easily be rearranged to be a formula for u+in terms of the other IENTIFIC quantities.The FTCS scheme is illustrated in Figure 19.1.1.It's a fine example of 6 an algorithm that is easy to derive,takes little storage,and executes quickly.Too bad it doesn't work!(See below.) The FTCS representation is an explicit scheme.This means thatfor each jcan be calculated explicitly from the quantities that are already known.Later we shall meet implicit schemes,which require us to solve implicit equations coupling the u for various j.(Explicit and implicit methods for ordinary differential 10.621 equations were discussed in 816.6.)The FTCS algorithm is also an example of a single-level scheme,since only values at time level n have to be stored to find E喜 43106 values at time level n+1. (outside von Neumann Stability Analysis North Software. Unfortunately,equation(19.1.11)is of very limited usefulness.It is an unstable method,which can be used only (if at all)to study waves for a short fraction of 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 coefficients of the difference equations are so slowly varying as to be considered constant in space and time.In that case,the independent solutions,or eigenmodes,of the difference equations are all of the form un=gneikjAr (19.1.12)
836 Chapter 19. Partial Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). t or n x or j FTCS Figure 19.1.1. Representation of the Forward Time Centered 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 points that are used to calculate spatial derivatives; the dashed lines connect points that are used to calculate time derivatives. The FTCS scheme is generally unstable for hyperbolic problems and cannot usually be used. quantities at timestep n + 1 in terms of only quantities known at timestep n. For the space derivative, we can use a second-order representation still using only quantities known at timestep n: ∂u ∂x j,n = un j+1 − un j−1 2∆x + O(∆x2) (19.1.10) The resulting finite-difference approximation to equation (19.1.6) is called the FTCS representation (Forward Time Centered Space), un+1 j − un j ∆t = −v un j+1 − un j−1 2∆x (19.1.11) which can easily be rearranged to be a formula for u n+1 j in terms of the other quantities. The FTCS scheme is illustrated in Figure 19.1.1. It’s a fine 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 representation is an explicit scheme. This means that un+1 j for each j can be calculated explicitly from the quantities that are already known. Later we shall meet implicit schemes, which require us to solve implicit equations coupling the un+1 j for various j. (Explicit and implicit methods for ordinary differential equations were discussed in §16.6.) The FTCS algorithm is also an example of a single-level scheme, since only values at time level n have to be stored to find values at time level n + 1. von Neumann Stability Analysis Unfortunately, equation (19.1.11) is of very limited usefulness. It is an unstable method, which can be used only (if at all) to study waves for a short fraction of 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 coefficients of the difference equations are so slowly varying as to be considered constant in space and time. In that case, the independent solutions, or eigenmodes, of the difference equations are all of the form un j = ξneikj∆x (19.1.12)
19.1 Flux-Conservative Initial Value Problems 837 ,0 Lax x or/ 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 k is 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 nted for number Therefore,the difference equations are unstable (have exponentially growing modes)if (k)>1 for some k.The number g is called the amplification factor at a given wave number k. To find (k),we simply substitute (19.1.12)back into (19.1.11).Dividing by gr,we get 2 )=1-i Ar sinkAr (19.1.13) whose modulus is 1 for allk:so the FTCS scheme is unconditionally unstable. If the velocity v were a function of t and r,then we would write v in equation (19.1.11).In the von Neumann stability analysis we would still treat v as a constant, the idea being that for v slowly varying the analysis is local.In fact,even in the SCIENTIFIC( case of strictly constant v,the von Neumann analysis does not rigorously treat the end effects atj=0 and j=N. More generally,if the equation's right-hand side were nonlinear in u,then a von Neumann analysis would linearize by writing u=uo+ou,expanding to linear order in ou.Assuming that the uo quantities already satisfy the difference equation 是 exactly,the analysis would look for an unstable eigenmode of ou. Despite its lack ofrigor,the von Neumann method generally gives valid answers Numerica 10.621 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 431 analysis.) (outside Recipes Lax Method North Software. The instability in the FTCS method can be cured by a simple change due to Lax. One replaces the term u?in the time derivative term by its average(Figure 19.1.2): 1 u巧一2(+1+u-) (19.1.14 This turns (19.1.11)into 吗1-+-小-器-- (19.1.15)
19.1 Flux-Conservative Initial Value Problems 837 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). t or n x or j Lax 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 k is 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)| > 1 for some k. The number ξ is called the amplification factor at a given wave number k. To find ξ(k), we simply substitute (19.1.12) back into (19.1.11). Dividing by ξn, we get ξ(k)=1 − i v∆t ∆x sin k∆x (19.1.13) whose modulus is > 1 for all k; so the FTCS scheme is unconditionally unstable. If the velocity v were a function of t and x, then we would write v n j in equation (19.1.11). In the von Neumann stability analysis we would still treat v as a constant, the idea being that for v slowly 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 = 0 and j = N. More generally, if the equation’s right-hand side were nonlinear in u, then a von Neumann analysis would linearize by writing u = u0 + δu, expanding to linear order in δu. Assuming that the u0 quantities already satisfy the difference equation exactly, the analysis would look for an unstable eigenmode of δu. Despite its lack of rigor, the von Neumann method generally gives valid answers 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 The instability in the FTCS method can be cured by a simple change due to Lax. One replaces the term un j in the time derivative term by its average (Figure 19.1.2): un j → 1 2 un j+1 + un j−1 (19.1.14) This turns (19.1.11) into un+1 j = 1 2 un j+1 + un j−1 − v∆t 2∆x un j+1 − un j−1 (19.1.15)
838 Chapter 19.Partial Differential Equations stable unstable or n △x x or (a) (b) 一2 83 Figure 19.1.3.Courant condition for stability of a differencing scheme.The solution of a hyperbolic 鱼 19881992 problem at a point depends on information within some domain of dependency to the past,shown here shaded.The differencing scheme (19.1.15)has its own domain of dependency determined by the choice 1.800 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 differencing domain of dependency is larger than that of the PDEs,as in (a),and unstable if the relationship is the from NUMERICAL RECIPESI reverse,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 amplification factor ξ=cosk△x-i v△t Press. sink△x (19.1.16) △x The stability condition 2<1 leads to the requirement Programs lv△t ≤1 (19.1.17) SCIENTIFIC △x This is the famous Courant-Friedrichs-Lewy stability criterion, often 61 called simply the Courant condition.Intuitively,the stability condition can be understood as follows(Figure 19.1.3):The quantity u in equation(19.1.15)is computed from information at points j-1 and j+1 at time n.In other words, j-1 and j+are the boundaries ofthe spatial region that is allowed to communicate information to Now recall that in the continuum wave equation,information Fuurggoglrion Numerical Recipes 10621 actually propagates with a maximum velocity v.If the point u is outside of the shaded region in Figure 19.1.3,then it requires information from points more 4310 distant than the differencing scheme allows.Lack of that information gives rise to g an instability.Therefore,At cannot be made too large. (outside The surprising result,that the simple replacement(19.1.14)stabilizes the FTCS scheme,is our first encounter with the fact that differencing PDEs is an art as much North Software. as a science.To see if we can demystify the art somewhat,let us compare the FTCS and Lax schemes by rewriting equation(19.1.15)so that it is in the form of equation (19.1.11)with a remainder term: -(）+(m+ (19.1.18) △t △t But this is exactly the FTCS representation of the equation 0u0u,(△x)2 72u (19.1.19) 2△t
838 Chapter 19. Partial Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). t or n ∆t x or j ∆t ∆x ∆x stable unstable (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 here shaded. 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 differencing domain of dependency is larger than that of the PDEs, as in (a), and unstable if the relationship is the reverse, 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 amplification factor ξ = cos k∆x − i v∆t ∆x sin k∆x (19.1.16) The stability condition |ξ| 2 ≤ 1 leads 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 u n+1 j in equation (19.1.15) is computed from information at points j − 1 and j + 1 at time n. In other words, xj−1 and xj+1 are the boundaries of the spatial region that is allowed to communicate information to un+1 j . Now recall that in the continuum wave equation, information actually propagates with a maximum velocity v. If the point u n+1 j is outside of the shaded region in Figure 19.1.3, then it requires information from points more distant than the differencing scheme allows. Lack of that information gives rise to an instability. Therefore, ∆t cannot be made too large. The surprising result, that the simple replacement (19.1.14) stabilizes the FTCS scheme, is our first encounter with the fact that differencing PDEs is an art as much as a science. To see if we can demystify the art somewhat, let us compare the FTCS and Lax schemes by rewriting equation (19.1.15) so that it is in the form of equation (19.1.11) with a remainder term: un+1 j − un j ∆t = −v un j+1 − un j−1 2∆x + 1 2 un j+1 − 2un j + un j−1 ∆t (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)