On dispersive difference schemes

Physica 18D (1986) 250-254 North-Holland, Amsterdam

O N DISPERSIVE DIFFERENCE S C H E M E S Peter D. LAX* New York University, Courant Institute of Mathematical Sciences, New York, N Y 10012, USA

1. Introduction

as e tends to zero, while the initial value of u,

In 1944 [4], von Neumann wrote a report outlining a method for calculating one-dimensional compressible flows with shocks. The method used Lagrangian description, the space and time derivatives were discretized symmetrically, and viscosity and heat conduction were set equal to zero. Calculations carried out at the Ballistics Research Laboratory at Aberdeen, using punched card equipment, produced approximate solutions with post-shock oscillations on mesh scale, von N e u m a n n boldly conjectured that the mesh-scale oscillations in velocity can be interpreted as heat energy produced by the irreversible action of the shock wave, and that as A x ~ 0, the approximations converge in the weak sense to the appropriate discontinuous solution of the equations. I believe that von Neumann was wrong in this assertion; the purpose of this note is to present my reasons why I believe the assertion to be wrong, and to analyze them as far as possible at this writing.

u(x,O)=uo(X )

(2.2)

is kept fixed. A far more detailed study of this limit was carried out by S. Venakides [7,8]. These investigations show that, as foretold by numerical calculations, the solutions u(x, t, e) of (2.1), (2.2) are oscillatory, and tend to a limit as e tends to zero in the weak sense only; we denote this weak limit as fi(x, t). Now rewrite (2.1) in conservation form,

utq-(~u )x"l-£2Uxxx=O, 1

2

(2.3)

and take the limit of each term in the sense of distributions. The first term tends to f,; the second term tends to ½(U2)x, where u 2 is the weak limit of u2(x, t, e). The third term tends to 0. So we obtain the relation fit + 12~[~-~]l x = 0.

(2.4)

Since f is the weak, but not strong, limit of

u(x,t,e), 2. Dispersive limits

m

U2=~f2.

In a series of papers [2], the author and C.D. Levermore have investigated the limit of solutions of the KdV equation

This shows that the weak limit fi does not satisfy eq. (2.3) with e = 0, i.e. that u~+ (½fi2)x ~ 0.

U t + UU x + e 2 U x x x = 0

(2.5)

(2.6)

(2.1)

*Work supported by Dept. of Energy under Contract DEAC02-76ER03077.

This is indeed borne out by the explicit formula obtained in [2] for ft. A difference scheme that is entirely centered and without viscous terms, genuine or artificial, is anal-

0167-2789/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

251

P.D. L a x ~ O n dispersive difference schemes

ogous to a modification of the original equations by a dispersive term. It was this analogy that led me to disbelieve von Neumann's conjecture; but I have. not been able to push through the argument, the finite difference case being, as so often, trickier than the differential case. Here is the story: The underlying equations governing o n e dimensional compressible flow are, in Lagrange coordinates, Vt + p~ = 0,

(2.7)

Vt - u~ = 0,

(2.8)

et + pu~ = 0,

(2.9)

where u is velocity, p is pressure, V is specific volume and e is specific internal energy. This system is made determined by adjoining the equation of state

p =p(e,V).

(2.10)

Eq. (2.7) is conservation of momentum, (2.8) is conservation of mass. The law of conservation of energy is obtained by multiplying (2.7) by u and adding it to (2.9):

(e+ r ) t + ( u p ) ~ = O ,

(2.11)

where K denotes kinetic energy:

K = u2/2.

*Actually this is the difference scheme used in [5].

u,+1/2_ k - - ,,,-1/2_X(p,k+l/2_pT,_l/2)" ~k

(2.13)

The difference approximations to (2.8) and (2.9) are centered at (k + 1/2) a~,(n + 1/2) At,

Vk,+l +1/2

u~+X/2) = V k,+ l / 2 + " \ ~ kx[,,,+1/2 +l --

(2.14)

and

en+l

n

)~,n+l/2(,,n+1/2

k+l/2= ek+l/2--"Fk+l/2~,~k+l

u~+l/2), (2.15)

where n+l/2

Pk+l/2

=

1(

n+X-

n

Pk+l/z-rPk+l/2

)

,

(2.16)

and

= AtlAS.

(2.17)

Extend the lattice functions u~,+x/2 etc. to all real 4, real positive t as piecewise constant over rectangles of size (A~, At). Denote by w(x, t) a smooth test function with compact support in t > 0. Multiply eq. (2.13) by w(kA~, nSt)A~ and sum over k, n. Assuming that u and p, as extended, have weak limits fi and fi as A~, At ---, 0, we conclude, after summation by parts, see [3], that

(2.12)

von N e u m a n n ' s difference scheme*, roughly speaking, can be described as follows: The approximate values of velocity are associated with positions that are integer multiples of Af and i n t e g e r + h a l f multiples of At, while the thermodynamical variables are associated with positions that are integer+ half multiples of Af, integer multiples of At,

u~,+1/2, (e,p,V)~+l/2.

Eq. (2.7) is differenced by centering at kAy, nAt,

f f (-fiwt+ pw¢)d~dt=O. This implies that fit + fi~ = 0

(2.18)

in the distribution sense. Similarly, we deduce from (2.14) that V, - fi~ = O.

(2.19)

We cannot apply this procedure to eq. (2.15) since, unlike (2.13) and (2.14), it is not in conservation

P.D. Lax/ On dispersivedifferenceschemes

252

form. However, Trulio and Trigger [6], have shown how to derive from (2.13)-(2.15) an energy equation in conservation*. Since [6] is not readily available, we reproduce here their argument: Add (2.13), and (2.13),+1, and multiply by ~..+1/2. using the abbreviation (2.16) and 2U k 1 n--1/2 n+1/2 K k = ~u k uk , n

Similarly, we can combine (2.16) and (2.25) to write

p=p~+a/2 =

T+T-I ) 2 P"

T(e+K)-T

l(e+K)

+ XS(Pu) - X S - I ( P u ) = 0.

K ; +1 = K ; -- ½)~( p,~ + 1/2 _ ,,, +,/2 ~u~,+ 1/2 1/2 k ' k - 1/2 ] (2.21) Averaging (2.21) over k and k + 1, and using the abbreviation n

1

n

n

(2.27)

We then rewrite (2.23) as

(2.20)

we get

(S+S-1)( 2

(2.28)

Suppose the sequence of approximations u, p and e remain uniformly bounded. Then so are K and P, and therefore a subsequence of e + K and Pu will have a weak limit. We denote these weak limits by e + K and Pu; since eq. (2.28) is in conservation form, we can deduce as before that

(2.22)

Kk+l/2 = ~ ( K k + Kk+l/2),

(e~)t

+ (P--u) ~ = 0.

(2.29)

we get

gn+l

__r~n+l/2"~

n 1 ( nn+l/2 k+l/2=Kk+l/2---~k~l-'k+3/2

n+l/2

l-'k+l/2}Uk+l

lX[ _n+1/2

nn+l/2"~un+l/2 --2"\tlk+l/2--Uk-1/2] k

•

(2.22')

Adding this to (2.15) gives En+l

k+l/2 ~

( ~ + ½fiz)t+ (fifi)~= 0;

1fi2 = K,

n+1/2_ n + l / 2

Pk+l /~k+l

(2.30)

this is identical with (2.29) if

n Ek+l/2

--~k(

The law of conservation of energy requires on the other hand that

n+1/2

--Pk

(2.31)

n+l/2"~

Uk

),

and

(2.23) tiff =

where E is total energy, E ~ + I / 2 = e~,+1/2 + K~+1/2,

(2.24)

P--~.

(2.32)

Suppose that the sequence of approximations u converges to fi weakly but not strongly. Set

and

p~+1/2

1[ , , + 1 / 2 - ,+1/2~ = ~L t'k+l/2 -t'Pk-1/2)-

u=Ft+m. (2.25) By (2.26),

Denote by $ translation in ~ by A~/2, and by T translation in t by At~2. We can combine formulas (2.20) and (2.22) to write

K= ½(T + T-X)(Sfi + Sm)(S-~fi + S-~m). Take the weak limit of both sides; we get

I(T+T KF'+I/2 = 2 2

1)(

St/)(S- in)"

(2.26)

*I thank Bill Noh for callingmy attention to this work.

K = ½~2 + ½(T + T - 1 ) ( S m ) ( S - a m ) .

(2.33)

If m consisted mostly of mesh-scale oscillations,

P.D. Lax/On dispersive difference schemes

253

then the second term on the right in (2.33) could indeed be zero, and relation (2.31) would hold. Similarly, if

odic in k, with period L / A - - K . Sum (3.3) with respect to k over such a period; we obtain that

p-p

Y~ log u k

K

(3.4)

1

consisted mostly of mesh scale oscillations, then P given by (2.27) could tend strongly to ~, and so relation (2.32) would hold. Thus it is very hard, without a more detailed understanding of the nature of the oscillations, to exclude the possibility that von N e u m a n n was right, after all, and the weak limits satisfy the law of conservation of energy.

3. The model equation It is traditional in computational fluid dynamics, when one is unable to prove anything about the behavior of a scheme for the real equations, to turn to the model equation

ut dp ( lu2 ) ~ = O.

(3.1)

Accordingly, the author and Jonathan Goodman have undertaken to study the following semidiscrete approximations to (3.1), in conservation form: 1

Otldk q- -'~(Uk+lU k -- UkUk_l) = 0,

is independent of t. On the other hand, summing (3.2) over a period shows that K

~.u k

(3.5)

1

is independent of t. At any time when all the Uk(t ) are positive, (3.5) yields an upper bound for each uk(t ) that is independent of t. Using this, (3.4) yields a lower bound for log Uk(t), which shows that uk(t ) > d > 0 for all t. Extend now the lattice function uk(t ) to all real by setting

U a ( ~ , t ) = u k ( t ),

I~-kAI
(3.6)

Denote by S translation in ~ by A. If uk(t ) and Ukl(t) were bounded uniformly for all t as A ---, 0, we could select a subsequence such that u a, ( $ u a ) ( S - lua) and log ua converge weakly to limits that we denote as fi u (2) and log u. Since both (3.2) and (3.3) are in conservation form, it follows by the standard argument used above that weakly

(3.2) fit + ½ ( ~ ) , =

0

(3.7)

(i-ff~), + fi~ = 0.

(3.8)

where we use the abbreviation and

a =a~. Assuming that u k ~ 0 we divide (3.2) by u k to obtain 0flog uk + ~1- (uk+ 1 - uk-1 ) = 0.

(3.3)

Note that (3.3) also is in conservation form. Suppose the solutions u we are approximating are periodic, with period L, and that the prescribed initial values are positive. Then u~ is peri-

If u a converged strongly to f, then U(2) = f 2 and

log u = log f ;

then eq. (3.7) and (3.8) become

(log ~), + ~ = o.

254

P.D. Lax / On dispersive differenceschemes

T h e s e e q u a t i o n s are incompatible if ~ c o n t a i n s a d i s c o n t i n u i t y ; therefore they can hold only if the initial d a t a p e r m i t a c o n t i n u o u s solution. C o n versely, if the initial data are such that (3.1) has n o c o n t i n u o u s solution, then u a does not converge strongly. E v e n if ua t e n d s to ~ only weakly, ( S u ~ ) ( S - l u a ) c o u l d converge to u 2, as the analysis presented at the e n d of section 2 shows. T o decide this matter, o n e n e e d s m o r e detailed knowledge of the oscillations of ua. Here it is of help to note that the system (3.3) is completely integrable, as shown in [1] b y K a c a n d v a n Moerbeke. The author a n d J o n a t h a n G o o d m a n hope to use the explicit solut i o n p r o v i d e d i n [1] to study the limit A ~ 0. This p a p e r is dedicated to M a r t i n Kruskal, with affection a n d a d m i r a t i o n , a n d in trust of his tolera n c e for speculative studies, as long as the hypotheses are clearly stated.

References [1] M. Kac and P. van Moerbeke, On an Explicitly Soluble System of Nonlinear Differential Equations related to Certain Toda Lattices, Adv. in Math. 16 (1975) 160-169. [2] P.D. Lax and C.D. Levermore, The Small Dispersion Limit of the KdV Equation, Comm. Pure Appl. Math., vol. XXXVI (1983) I p. 253-290, II p. 571-593, III p. 809-830. [3] P.D. Lax and B. Wendroff, Systems of Conservation Laws, Comm. Pure Appl. Math. III (1960) 217-237. [4] J. v. Neumann, Proposal and Analysis of a New Numerical Method for the Treatment of Hydrodynamical Shock Problems, vol. VI, Collected Works (Pergamon, London). [5] J. v. Neumann and R.D. Richtmyer, A Method for the Numerical Calculation of hydrodynamicalShocks, J. Appl. Phys. 21 (1950) 380-385. [6] J.G. Trulio and K.R. Trigger, Numerical Solution of One Dimensional Hydrodynamical Shock Problem, UCRL Report 6522, 1961. [7] S. Venakides, The Zero Dispersion Limit of the KdV Equation with Nontrivial Reflection Coefficient, Comm. Pure Appl. Math. XXXVIII (1985) 125-155. [8] S. Venakides, The Generation of Modulated Wavetrains in the Solution of the KdV Equation, to appear in Comm. Pure Appl. Math.