c19-1 (779616), страница 3

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 3)

In some fluid problems, however, shock formation is not merely an annoyance,but an actual physical behavior of the fluid whose detailed study is a goal. Then,numerical viscosity alone may not be adequate or sufficiently controllable. This is acomplicated subject which we discuss further in the subsection on fluid dynamics,below.For wave equations, propagation errors (amplitude or phase) are usually mostworrisome. For advective equations, on the other hand, transport errors are usuallyof greater concern. In the Lax scheme, equation (19.1.15), a disturbance in theadvected quantity u at mesh point j propagates to mesh points j + 1 and j − 1 atthe next timestep. In reality, however, if the velocity v is positive then only meshpoint j + 1 should be affected.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).Thus far we have been concerned with amplitude error, because of its intimateconnection with the stability or instability of a differencing scheme. Other varietiesof error are relevant when we shift our concern to accuracy, rather than stability.Finite-difference schemes for hyperbolic equations can exhibit dispersion, orphase errors. For example, equation (19.1.16) can be rewritten as19.1 Flux-Conservative Initial Value Problems841vt or nvx or jFigure 19.1.4.

Representation of upwind differencing schemes. The upper scheme is stable when theadvection constant v is negative, as shown; the lower scheme is stable when the advection constant v ispositive, also as shown. The Courant condition must, of course, also be satisfied.The simplest way to model the transport properties “better” is to use upwinddifferencing (see Figure 19.1.4): nuj − unj−1n+1,nuj − uj∆x= −vjn un − un∆tj+1j,∆xvjn > 0(19.1.27)vjn < 0Note that this scheme is only first-order, not second-order, accurate in thecalculation of the spatial derivatives. How can it be “better”? The answer isone that annoys the mathematicians: The goal of numerical simulations is notalways “accuracy” in a strictly mathematical sense, but sometimes “fidelity” to theunderlying physics in a sense that is looser and more pragmatic.

In such contexts,some kinds of error are much more tolerable than others. Upwind differencinggenerally adds fidelity to problems where the advected variables are liable to undergosudden changes of state, e.g., as they pass through shocks or other discontinuities.You will have to be guided by the specific nature of your own problem.For the differencing scheme (19.1.27), the amplification factor (for constant v) is v∆t (1 − cos k∆x) − i v∆t sin k∆xξ =1−∆x ∆xv∆t v∆t 1 − (1 − cos k∆x)|ξ|2 = 1 − 2 ∆x ∆x (19.1.28)(19.1.29)So the stability criterion |ξ|2 ≤ 1 is (again) simply the Courant condition (19.1.17).There are various ways of improving the accuracy of first-order upwinddifferencing.

In the continuum equation, material originally a distance v∆t awaySample 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).upwind842Chapter 19.Partial Differential Equationsstaggeredleapfrogx or jFigure 19.1.5.

Representation of the staggered leapfrog differencing scheme. Note that informationfrom two previous time slices is used in obtaining the desired point. This scheme is second-orderaccurate in both space and time.arrives at a given point after a time interval ∆t. In the first-order method, thematerial always arrives from ∆x away. If v∆t ∆x (to insure accuracy), this cancause a large error.

One way of reducing this error is to interpolate u between j − 1and j before transporting it. This gives effectively a second-order method. Variousschemes for second-order upwind differencing are discussed and compared in [2-3].Second-Order Accuracy in TimeWhen using a method that is first-order accurate in time but second-orderaccurate in space, one generally has to take v∆t significantly smaller than ∆x toachieve desired accuracy, say, by at least a factor of 5.

Thus the Courant conditionis not actually the limiting factor with such schemes in practice. However, there areschemes that are second-order accurate in both space and time, and these can often bepushed right to their stability limit, with correspondingly smaller computation times.For example, the staggered leapfrog method for the conservation equation(19.1.1) is defined as follows (Figure 19.1.5): Using the values of un at time tn ,compute the fluxes Fjn .

Then compute new values un+1 using the time-centeredvalues of the fluxes:− un−1=−un+1jj∆t nn(F− Fj−1)∆x 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 requiresthat un−1 and un be stored to compute un+1 .For our simple model equation (19.1.6), staggered leapfrog takes the formun+1− un−1=−jjv∆t n(u− unj−1 )∆x j+1(19.1.31)The von Neumann stability analysis now gives a quadratic equation for ξ, rather thana linear one, because of the occurrence of three consecutive powers of ξ when theSample 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 Problems843form (19.1.12) for an eigenmode is substituted into equation (19.1.31),ξ 2 − 1 = −2iξv∆tsin k∆x∆x(19.1.32)whose solution is1−v∆tsin k∆x∆x2(19.1.33)Thus the Courant condition is again required for stability.

In fact, in equation(19.1.33), |ξ|2 = 1 for any v∆t ≤ ∆x. This is the great advantage of the staggeredleapfrog method: There is no amplitude dissipation.Staggered leapfrog differencing of equations like (19.1.20) is most transparentif the variables are centered on appropriate half-mesh points:nrj+1/2n+1/2sjnunj+1 − unj∂u ≡v=v∂x j+1/2∆xn+1/2n+1uj − unj∂u ≡=∂t j∆t(19.1.34)This is purely a notational convenience: we can think of the mesh on which r ands are defined as being twice as fine as the mesh on which the original variable u isdefined.

The leapfrog differencing of equation (19.1.20) isn+1n− rj+1/2rj+1/2n+1/2sjn+1/2=∆tn−1/2− sj=v∆tn+1/2− sj∆xnnrj+1/2− rj−1/2sj+1(19.1.35)∆xIf you substitute equation (19.1.22) in equation (19.1.35), you will find that onceagain the Courant condition is required for stability, and that there is no amplitudedissipation when it is satisfied.If we substitute equation (19.1.34) in equation (19.1.35), we find that equation(19.1.35) is equivalent to− 2unj + un−1un+1un − 2unj + unj−1jj2 j+1=v(∆t)2(∆x)2(19.1.36)This is just the “usual” second-order differencing of the wave equation (19.1.2). Wesee that it is a two-level scheme, requiring both un and un−1 to obtain un+1 .

Inequation (19.1.35) this shows up as both sn−1/2 and r n being needed to advancethe solution.For equations more complicated than our simple model equation, especiallynonlinear equations, the leapfrog method usually becomes unstable when the gradients get large. The instability is related to the fact that odd and even mesh points arecompletely decoupled, like the black and white squares of a chess board, as shownSample 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.

Характеристики

Список файлов книги