c19-1 (779616), страница 2
Текст из файла (страница 2)
We accordinglyadopt it exclusively. (See, for example, [1] for a discussion of other methods ofstability analysis.)Lax MethodThe instability in the FTCS method can be cured by a simple change due to Lax.One replaces the term unj in the time derivative term by its average (Figure 19.1.2):unj →1 nuj+1 + unj−12(19.1.14)This turns (19.1.11) intoun+1=j v∆t n1 nuu+ unj−1 −− unj−12 j+12∆x j+1(19.1.15)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).x or j838Chapter 19.Partial Differential Equationsstableunstable∆t∆t∆x∆xx or j(a)( b)Figure 19.1.3. Courant condition for stability of a differencing scheme. The solution of a hyperbolicproblem 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 choiceof points on one time slice (shown as connected solid dots) whose values are used in determining a newpoint (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 bedetermined simply by the outermost points.Substituting equation (19.1.12), we find for the amplification factorξ = cos k∆x − iv∆tsin k∆x∆x(19.1.16)The stability condition |ξ|2 ≤ 1 leads to the requirement|v|∆t≤1∆x(19.1.17)This is the famous Courant-Friedrichs-Lewy stability criterion, oftencalled simply the Courant condition.
Intuitively, the stability condition can bein equation (19.1.15) isunderstood as follows (Figure 19.1.3): The quantity un+1jcomputed 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. Now recall that in the continuum wave equation, informationinformation to un+1jactually propagates with a maximum velocity v. If the point un+1is outside ofjthe 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 toan instability.
Therefore, ∆t cannot be made too large.The surprising result, that the simple replacement (19.1.14) stabilizes the FTCSscheme, is our first encounter with the fact that differencing PDEs is an art as muchas a science. To see if we can demystify the art somewhat, let us compare theFTCS and Lax schemes by rewriting equation (19.1.15) so that it is in the form ofequation (19.1.11) with a remainder term: n− unjun+1uj+1 − unj−11 unj+1 − 2unj + unj−1j= −v+(19.1.18)∆t2∆x2∆tBut this is exactly the FTCS representation of the equation∂u (∆x)2 2∂u= −v+∇ u∂t∂x2∆t(19.1.19)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).t or n19.1 Flux-Conservative Initial Value Problems839The Lax method for this equation is1 nn(r+ rj−1)+2 j+11= (snj+1 + snj−1 ) +2rjn+1 =sn+1jv∆t n(s− snj−1 )2∆x j+1v∆t nn(r− rj−1)2∆x j+1(19.1.21)The von Neumann stability analysis now proceeds by assuming that the eigenmodeis of the following (vector) form, 0 nrjn ikj∆x r=ξ(19.1.22)esnjs0Here 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), and dividing by the power ξ n , gives the homogeneous vector equation(cos k∆x) − ξv∆tsin k∆xi∆x v∆tr0sin k∆x0∆x· = (cos k∆x) − ξs00i(19.1.23)This admits a solution only if the determinant of the matrix on the left vanishes, acondition easily shown to yield the two roots ξξ = cos k∆x ± iv∆tsin k∆x∆x(19.1.24)The stability condition is that both roots satisfy |ξ| ≤ 1. This again turns out to besimply the Courant condition (19.1.17).Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).where ∇2 = ∂ 2 /∂x2 in one dimension. We have, in effect, added a diffusion term tothe equation, or, if you recall the form of the Navier-Stokes equation for viscous fluidflow, a dissipative term. The Lax scheme is thus said to have numerical dissipation,or numerical viscosity. We can see this also in the amplification factor.
Unless |v|∆tis exactly equal to ∆x, |ξ| < 1 and the amplitude of the wave decreases spuriously.Isn’t a spurious decrease as bad as a spurious increase? No. The scales that wehope to study accurately are those that encompass many grid points, so that they havek∆x 1. (The spatial wave number k is defined by equation 19.1.12.) For thesescales, the amplification factor can be seen to be very close to one, in both the stableand unstable schemes. The stable and unstable schemes are therefore about equallyaccurate.
For the unstable scheme, however, short scales with k∆x ∼ 1, which weare 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. Both the stable and the unstable schemes are inaccurate for these shortwavelengths, but the inaccuracy is of a tolerable character when the scheme is stable.When the independent variable u is a vector, then the von Neumann analysisis slightly more complicated.
For example, we can consider equation (19.1.3),rewritten as ∂ vs∂ r=(19.1.20)∂t s∂x vr840Chapter 19.Partial Differential EquationsOther Varieties of Errorv∆tξ = e−ik∆x + i 1 −sin k∆x∆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),depending on their value of k. If ∆t = ∆x/v, then the exact solution for each modeof a wave packet f(x − vt) is obtained if each mode gets multiplied by exp(−ik∆x).For this value of ∆t, equation (19.1.25) shows that the finite-difference solutiongives the exact analytic result. However, if v∆t/∆x is not exactly 1, the phaserelations of the modes can become hopelessly garbled and the wave packet disperses.Note from (19.1.25) that the dispersion becomes large as soon as the wavelengthbecomes comparable to the grid spacing ∆x.A third type of error is one associated with nonlinear hyperbolic equations andis therefore sometimes called nonlinear instability.
For example, a piece of the Euleror Navier-Stokes equations for fluid flow looks like∂v∂v= −v+...∂t∂x(19.1.26)The nonlinear term in v can cause a transfer of energy in Fourier space fromlong wavelengths to short wavelengths. This results in a wave profile steepeninguntil a vertical profile or “shock” develops. Since the von Neumann analysissuggests that the stability can depend on k∆x, a scheme that was stable for shallowprofiles can become unstable for steep profiles. This kind of difficulty arises ina differencing scheme where the cascade in Fourier space is halted at the shortestwavelength representable on the grid, that is, at k ∼ 1/∆x. If energy simplyaccumulates in these modes, it eventually swamps the energy in the long wavelengthmodes of interest.Nonlinear instability and shock formation is thus somewhat controlled bynumerical viscosity such as that discussed in connection with equation (19.1.18)above.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
