- Email: [email protected]
Finite difference solutions of dispersive partial differential equations
268
Mathematics
FINITE DIFFERENCE EQUATIONS
SOLUTIONS
and Computers
OF DISPERSIVE
in Simulation
XXV (1983) 268-278 North-Holland
PARTIAL DIFFERENTIAL
William L. BRIGGS Department
of Mathematics
and Computer Science, Clarkson College of Technology, Potsdam,
NY 13676, U.S.A.
Talib SARIE Department
of Mathematics,
Jordan University, Amman,
Jordan
The second order leapfrog method is used to discretize the linearized KdV equation which is itself a dispersive partial differential equation. The resulting difference equation is solved and analyzed in terms of its dispersion relation and propagation properties. Numerical experiments are included to illustrate some of the results obtained. Finally, a novel aliasing error which affects the propagation of wave packets is presented
1.
INTRODUCTION
It might be thought that the analysis of finite difference schemes with constant coefficients is Inan exercise which has long been exhausted. deed the earliest investigations, which were concerned primarily with the stability of Cerschemes, date back at least thirty years. tainly at that time it was recognized that even stable finite difference schemes, applied to time dependent problems, necessarily introduce errors in resolving the amplitude and phase of The orooerties of these damoinq and solutions. dispersion errors has been the subject of &IStained research ever since. Despite this level of activity, only recently has the role of group velocity in the analysis of difference schemes been recognized. This study is undertaken in the belief that, while the analysis of difference schemes is an area that has been wellcovered, there are still some details that could use further illumination.
In the followino oaoes, we will be concerned with the numerica! solution of a dispersive partial differential equation. Many of the conclusions and results would apply to any equation However the choice of a disperof propagation. sive equation does provide some new effects and The importance of propagation properinsights. ties (e.g. dispersion and group velocity) in understanding the behavior of difference schemes seems under-emphasized in much of the numerical Although some of the analysis literature. results which follow are not new or unexpected, they do appear in a new setting and it serves at least a pedagogic purpose to present them again. Much of the work which directly precedes and anticipates this present study-was done to develop schemes for fluid flow problems. Among this early work, the papers which deal specifically with phase error and its control seem In the 1963 work of most appropriate here. 0378-4754/83/$3.00 0 1983, IMACS/Elsevier
Science Publishers
Stone and Brian [6], a family of schemes for a one-dimensional convection-diffusion problem was presented. Although the terminology was not present, the dispersion relation for the schemes was derived and an optimal scheme which minimizes damping and phase errors was developed. In 1966, Roberts and Weiss [5] pointed out the role of damping and dispersion in pure advection problems and developed a fourth order scheme (in space) which considerably reduces phase error. Crowley [l] compared the damping and dispersion properties of several existing schemes, particularly for advection problems in two dimensions. At about the same time, Fromm [4] analyzed split time step methods and devised his well-known method of zero average phase error, again for two dimensional advection problems. All of these studies clearly recognized the fact that while damoinq errors can be avoided by usinq time-centered (leapfrog) difference schemes; all difference schemes are dispersive. That is, when any finite difference scheme is applied to an equation of propagation, each spatial mode of the difference scheme has a phase speed which depends upon the wavelength of the mode. This is true whether the original partial differential equation has this property or not. The earliest recognition that group velocity plays a role in the analysis of finite difference schemes seems to have occurred in 1971 when Vichnevetsky and his colleagues [9,10] considered the one-dimensional advection equation and examined semi-discrete (time continuous) Using Fourier analysis, propagation schemes. properties of these schemes were derived, the important notions of group velocity and energy propagation of the finite difference scheme were introduced. and the role of soatial comoutational modes was illustrated through numerical experiments. Shortlv thereafter Grotjahn [2,3] studied group velocity in finite difference schemes arising from one and two dimensional
B.V. (North-Holland)
269
W.L. Briggs, T. Sarie / Finite difference solutions of dispersive PDEs
More recently, Trefethen geophysical problems. [7,8] has studied the one-dimensional advection equation modeled by fully discrete schemes. He derived the relevant dispersion relations and in a series of imaginative numerical experiments, clearly showed the behavior of the computational mode, pointinq out that a fully discrete scheme can support computational modes in both space and time. In addition, known stability criteria were interpreted in terms of the group velocity. In the present order leapfrog equation
study, we consider the second method applied to the dispersive
ut + Uux + buxxx = 0.
U,be
R
This can be regarded as the linear analog of the Korteweg-deVries (KdV) equation. A brief discussion in the next section presents the propagation properties of this equation. Most important for our purposes is the fact that the dispersion relation, which relates the frequency w of a spatial mode to its wave number u, is a nonlinear function of u. This coupled with the fact that the difference scheme has its own nonlinear dispersion relation means that discrepancies between the propagation properties of the continuous and discrete problems are even greater than with a non-dispersive partial differFor example, the continuous ential equation. problem can support wave trains with zero phase speed and wave packets with zero group velocity. We will see that the difference scheme approximates these stationary modes, but in the process, produces spurious stationary modes of its own.
2.
THE CONTINUOUS
In this section we will investigate finite difference,solutions to the dispersive partial differential equation ut + lJux + bu
In Section 4, several numerical experiments are performed to illustrate some of the ideas of the previous section. The presence of computational modes in a calculation is explored and an argument is given to show that the choice of initial conditions at the first two time levels ultimately determines whether computational modes appear. The paper concludes by pointing out an aliasing error that can occur when wave packets are propagated by finite difference Whether this error is a serious danger schemes. or merely a curiosity has yet to be determined.
xxx
= 0,
U,b
E
R
(2.1)
on the intervals 0 < x < ZIT, t > 0 subject to periodic boundary conditions and an initial condition u(x,O) = f(x). Without loss of generality it is assumed that U > 0, however consideration b > 0 and b < 0 leads to some interesting situations. A brief look at the propagation properties of continuous solutions to this problem will provide some background and also a point of comparison for the discrete solutions. Equation (2.1) admits normal mode solutions of the form u(x,t) = ei(ux-mt)
where the wave number u and the frequency w are related through the dispersion relation w = Uu - bu3 = -bu(u'
- ;).
(2.2)
This expression, in which u assumes integer values, gives the frequency w with which any The graphs of the dispersion wave travels. relation for the cases b > 0 and b < 0 are shown Note that in in Figures la and lb respectively. both cases w is an odd function of u, a property which will be preserved by the analogous discrete dispersion relations. The phase speed or speed of propagation of a pure mode with wave number u is given by c(u) = ti
In the third section of the paper, the difference equations arising from the Leapfrog scheme are solved subject to periodic boundary conditions using separation of variables. This approach seems neglected, but it has the advantage that it produces a general solution in terms of corresponding spatial and temporal modes and also gives stability conditions and the dispersion relation in a straightforward way. The phase speed and group velocity curves for the difference scheme are presented and discussed, particularly in relation to the corresponding curves for the continuous problem. The stability condition for the scheme proves to have some interesting structure which is also discussed.
PROBLEM
= U - bu2.
(2.3)
The graphs o/the phase speed for b > 0 and b < 0 are shown in Figures 2a and 2b respectiveIn both cases the phase speed is a nonly. constant function of p. which is to say that different modes travel at different speeds. Of special interest is the case b > 0 in which the phase speed is of both signs with a stationary mode at u= Jms. Finally, the group velocity or the speed at which a packet of waves of wave number u propagates is given by C(u) = 4
= U - 3bu2.
( 2 :4)
Again the group velocity is a non-trivial function of the wave number (Figures 3a and 3b) which gives rise to the term dispersive in describing this problem. In the case b > 0, the group velocity is of both signs. Energy in the long waves (small u) propagates in the positive direction. Energy in the p = Jm56 mode is stationary and energy in the short waves (large Both u) propagates in the negative direction. phase speed and group velocity are even functions of the wave number. 3.
ANALYSIS
OF THE DIFFERENCE
a.
Solution
by separation
EQUATIONS
of variables.
With this understanding of the continuous problem, we now turn to finite difference approxima-
270
W.L. Briggs, T. Sarie / Finite difference solutions of dispersive PDEs
tions to (2.1).
A uniform
grid of width
h = 'K1_
is established on the interval 0 < x < Zn, while the interval t > 0 is partitioned-at intervals of length k. In this study we shall be concerned only with finite difference representations of equation (2.1). The Leapfrog method is obtained by replacing all derivatives in (2.1) by second order finite differences. Letting u(m,n) be the discrete approximation to the exact solution at the point x = nh, t = mk, the following difference equations result. u(m+l,n)
- u(m-1,n)
+ a[u(m,n+l)
- u(m,n-1)]
+ B[u(m,n+2)-2u(m,n+l)+2u(m,n-1)-u(m,n-2)]=0 for 0 5 n 5 K-l, m > 1 with f(x,,) and u(l,n)
specified.
u(m,O)=u(m,K),u(O,n)= The two parameters
a = !$ and @ = k will have important parts to h3 play in all that follows. There are several approaches to solving linear, constant coefficient partial difference equations. We proceed by separation of variables. Letting u(m,n) = M(m)N(n) it follows that M(m+l) Y
1C
lb
20
(a)
25
- M(m-1)
= XM(m) and
N(n+2)+(a-2B)N(n+l)+~N(n)-(a-24)N(n-l)-N(n-2) -0 B 5 B where A is the separation constant which will be determined. Turning to the N equation, we seek solutions of the form N(n) = g". The periodic boundary conditions of the original equation appear as the condition N(0) = N(K). This requires that < be of the form
ii!??%! c=e
K
JJ = 0, -+ 1, _ + 2, . . . . $
.
We see immediately at this point that the spatial part of the discrete solution is quantized and consists
Figure
oispersion la)
relation
for
b = + 1.52 I lo-'
the
is constant
n = Fmode
has a wavelength
modes.
The
over the grid while
the
of 2h‘. This is, of
course, in contrast to the continuous problem in which modes of all wavenumbers are available For the discrete problem, only to the system. those initial conditions which are a linear com-
1 co”tin”o”s (b)
of only $ distinct
n = 0 mode
problem,
b = -1.52 x 10
” = 1 -2
bination
of these t modes
can be represented
An initial condition which lies outexactly. side of this space will not be represented or propagated without error. We can then express the spatial part of the solution as a linear combination of modes of the form N(,,)=e
?_.
K =elpnh
= 0, -+ 1, _ + 2,..., t;
.
Substituting this into the original equation for N(n) allows the separation constant X to be determined for each spatial mode. We find that A(p) = -Zi[(a-2B)sinuh
+ Bin
2uh].
271
W. L. Briggs, T. Sarie / Finite difference solutions of dispersme PDEs Using these values assuming
of X in the M(m) equation
a solution x&z
v+(u) -
= -is(u) +
= 2
and
of the form M(m) = ym gives l-g(u)
2
The fact where g(u) = (a-2a)sinuh + @sin 2ph. that there are two values of y associated with each mode arises, of course, from the necessity of solving a second order difference equation to determine 14. This will have important implications in all that follows. The temporal
Y+(V) which
behavior of mode u is determined is turn depends upon g(u). Three
cases
must be considered.
(i)
If lg(lJ) > 1 for then IY u) 1 > 1
by
at least one value of u. and hence solutions can qrow with increasing m (time). These are instable solutions which are not consistent with the continuous solutions.
I
(ii)
If ig(u)I = 1 for at least one value of u, then y, = y_ = i or y, = y_ = -i. In either case the difference equation for M has a double root and mode u has temporal behavior
of the form M(m) = clim + c2mim.
This slow linear growth, while not catastropic, is also inconsistent with the continuous solutions and must be avoided. (iii) iG+lz\l;'l Y Bn~~lf'lo~~~"~~v~'a~~o~~~~nt The conamFlitude, oscillatory behavior. dition ig(p)I < 1 imposes restrictions on the parameters CL and B which amount to stability conditions for each mode and for the entire scheme. We will return to the precise form of these stability conditions shortly. Clearly a useful difference scheme results only in this last case. b.
Discrete
dispersion
If we assume that (g(p)1 < 1 for all wave numbers u, then the temporal part of the solution for mode p takes the form
M(m) = cle
-iw2mk +ce
2
where o,k = sin- '{(a-2a)sinuh 0 (~,k
~II
and
+ Bsin 2u.h)
IJ = 0, + 1,
+
(2.5)
w2k = n-w,k.
The full solution of (2.1) associated u finally appears as i(unh-wlmk) i(unh-m2mk) u(m,n)=c,e +ce 2 2,...,
+- . 2
Rewriting
with mode
(2.5) in the form
sinwk = (a-2b)sinph
+ f3sin2ph
(2.6)
gives the dispersion relation for the difference As shown in Figure 0, this relation equations. is double valued with respect to both w and u and hence both w, and w2 are represented. The figure also shows the continuous dispersion relation. For small wave numbers, the continuous dispersion relation is matched very closely by the w, branch of the discrete dispersion relation.
For this reason,
the o, component
of the
discrete solution is called the physical mode since it corresponds and converges to the continous solution. The w2 component of the discrete solution is entirely spurious and is called the computational mode. To be completely accurate, these modes should be called the physical mode and the computational mode in time since they are two temporal modes corresponding to the same wave number. However, the dispersion relation could be read in the other direction also. If the frequency is fixed (say be time forcing at a boundary)then two spatial modes arise, one physical and one computational. Therefore for one (u,w) pair in the physical problem, the difference scheme is capable of producing four (u,w) pairs, only one of which approximates the physical pair. c.
Phase speed and group
velocity.
The propagation properties of various discrete modes may be obtained from the dispersion relation (2.5) or (2.6). The phase speed is given by c(u) = &
relation.
-iw,mk
These may be regarded as the normal modes of Any linear combination of the discrete problem. these modes can be represented by the difference scheme.
sin-'[(a-26)sinph
+ Bsin2uh].
Expanding this expression for small wave numbers (uh << 1) we find that for the physical mode c(,,) = "_bu2 _ +
(uh)2 + O(uh)4.
The difference scheme gives a second order approximation to the phase speed of long waves. Figure 5 is a graph of the phase speed of the physical mode as it varies with u. The phase speed curve for the corresponding continuous problem is also shown for comparison. As is well known, any difference scheme applied to a propagation equation will be dispersive. In this study, in which the continuous problem is also dispersive, the departures of the numerical phase speed from the continuous phase speed are even more pronounced. Of particular interest are the stationary modes. For b > 0, the continuous phase speed has a zero at n = 0 and uO = + m in the u-m plane
W. L. Briggs, T. Sarie / Finiie difference solutions
212
mare speeafor (a)
b-
the
+ 1.52 x 10
(bl
(bl
Figure 2
Figure 3
contin”o”s -2
prm,en, (b)
”
=
(a)
c
2o
IO _
for the b - t1.52 x 10-2
Group velocity
1
b = -1.52 x IO-'
2-
-I
of dispersive PDEs
continuousproblem. 0 = 1 (b)
b = -1.52 x 10-2
W.L. Briggs, T. Sarie / Finite difference solutions (Both the
yields
numerical and continuous phase speeds are even functions of u and negative wave numbers can be omitted from the discussion.) Now if Usb%O(l),
u0
oratuOh
= +m
in the ph-mk
plane.
then u0 % O(1) << $ and the stationary
.
This puts the
stationary wave near the short wave end of the Assuming that U and b are discrete spectrum. given and assuming that the zero of the phase speed is of physical importance, then h and k should be chosen so that u0 lies near the center of the discrete
spectrum.
Of course the difference scheme will not reproduce the stationary mode of the continuous problem even if it is centered in the discrete The discrete phase speed has zeros spectrum. at u=O, u = k and at
"
0(h2).
The above remarks pertain to the phase speed of By the symmetry of the the physical mode. numerical dispersion relation, it can be seen that the phase speed of the computational mode is equal in magnitude but opposite in sign to that of the physical mode. In many physical problems and numerical experiments the group velocity may be of more interReturning to the est than the phase speed. discrete dispersion relation the numerical group velocity may be computed to be C(u) = ;
(a-2!3)cosnh + 2BcosZiJh 1-[(a-2B)sinuh
=
3 + $ h2
+ 4sig
U-3bp2 - $(‘-a2)(ph2) + O(j&,)4_
Again the group velocity of the continuous lem is reproduced by the physical mode to
+ 0(h2)
The first value of u3
(2.7) prob-
The group velocity of the computational O(uh)2. mode is equal in magnitude but of opposite sign Figure 6 is a to that of the physical mode. graph of the exact group velocity and the numerical group velocity as it varies with u. As with phase speed, the group velocity of the continuous problem has much more structure than it For b > 0 would in a non-dispersive problem. the mode p. = Jm36 has zero group velocity and
corresponds
to the sta-
tionary mode of the continuous problem, while the second value is a spurious stationary mode at a higher wave number. The difference scheme clearly has the capability of moving energy in the wrong direction even in the intermediate wave numbers. When the computational mode is introduced the situation becomes even more comWe will return to these results in plicated. later numerical experiments. d.
Stability
conditions.
We now return to the question of stability and investigate the conditions which determine a stable difference scheme. It was seen that if /g(p)1 = I(a-2g)sinph for all wave numbers
+
273
$ + O(h2)
2=
wave may
not be resolved within the spectrum of the If U > b then the situation difference scheme. is even worse with the stationary wave longer than the fundamental wave of the difference At the other extreme, if U s O(1) and scheme. b 'irh*U then a%fi and 'JO Q $
of dispersive PDEs
+ BsinPuhl
< 1
0 2 u 5 !j then solutions
of
the difference scheme remain bounded as m (time) These inequalities generate regions increases. of stability in the (@,a) plane which are bounded by the pairs of parallel lines (a-2e)sinph
- BsinPuh
= + 1
or a = 2(1-cosph)b
+ cscuh
(2.8)
where ~1 > 0 and -m < B cm. Figure 7 shows the stability strips of several of the higher and more critical modes. For a given choice of o. and B to yield a stable scheme the point (B,a) must lie in the intersection of the stability regions of all modes. Looking at the slopes and intercepts of the boundary lines in (2.8) it is evident that the stability regions enlarge with decreasing wave number. It is the higher modes (uh = F,
z, 5, $) that determine
of the stability The boundaries given below:
the boundaries
region.
for these critical
modes
are
ZCc('+2_+3B 3 -J5 E:a=t,t2fi 2 -
as in the case of phase speed, h and k must be chosen carefully in order to resolve VU within
;:o=lZ+@
the spectrum of the difference scheme. However is resolved the difference scheme only appr even if "8 ximates it and furthermore produces a The zeros second stationary mode of its own. of the discrete group velocity are not simply expressed, however expanding in powers of h
;
2
: a = + VT +
(e-n)@.
The outline of the common region of stability of these modes is shown in the figure. The stability strips for all other modes contain this comnon region. We see that if u 2 1, then
W. L. Briggs, T. Sarie / Finite difference solutions
214
k
5-
-
a-
3-
Z-
I-
L) 5.8
h 5.2
6.6
0.0
I.6
of dispersive
PDEs
215
W. L. Briggs, T. Sarie / Finiie difference solutions of dispersive PDEs
then a stable 161 >$
scheme
is not possible,
= .39, then stability
while
if
is not possible.
A different picture of the stability conditions is obtained if the conditions Ig(p)I < 1 are expressed in terms of k and h. Some calculation eventually leads to the stability condition k<
h3
where
? ? UsinnhlRL-hLI
a2 =
2b l-cosph) (,, "
If vh is fixed to consider for the case b > 0. the stability of a single mode, then the stability region in the (h,k) plane will have the general appearance shown in Figure 8. The location of the asymptote at h = R depends upon the particular mode under consideration and on the relative values of the parameters b and U. For longer (more stable) waves the asymptote It is also interesting that shifts to the left. the choice of h = R, which would allow arbitrarily large time steps, is the combination of parameters for which the mode has zero phase The stability region for a single mode speed. in the case b < 0 is given by k<._
h3
1
0
* A I’
-1
Usinnh(R2
+ h2)
*t
.
II
-11
II
’
Again, the details of the curve will depend upon the choice of b and U and upon the mode under A general curve for this case is consideration. shown in Figure 8.
4.
EXPERIMENTS
AND CONSEQUENCES
In this section several numerical experiments are presented to illustrate the ideas of the In addition some examples are previous section. given to point out further consequences of using a discrete model to represent a dispersive partial differential equation.
"=l.b
- 1.52 x 10-2, k - 8.7 x 10-3
(a)
Initial cmditions.
(b)
u(nh. 1Ok)
"("h.0) = cor(o"h)cos(mh).
(1 = 2. B - 30
As in previous studies, we present some wave packet experiments to demonstrate the propagaIn the tion properties of the Leapfrog method. first two experiments, we use a laboratory with K=72 grid points, U=l, b=l.52x 10e2, k=8.7x10m3. In the first case, we consider an initial condition consisting of a wave with wave number B=24 (wavelength 3h) modulated by a long wave with wave number a=1 (wavelength 72h), i.e. This forms a series of u(x,O) = cosux cosox. wave packets which should propagate at the group velcoity C(24) = -25.5 according to the group velocity relation (2.4). Figure 9 shows the discrete solution at m=O and m=lO time steps. Clearly the packet is almost stationary which is consistent with the discrete group velocity (2.7) which qives C(24) s 0. In the second case.initial-conditions of the form u(x,O) = The cosax COSBX with a=2 and 0=30 are taken. continuous equation should propagate these wave packets with a group velocity of C(30) = Figure 10 shows the result of solving -41.0.
I -Nunrlcal
solution
0
-1
n
-2 0
100
200
300 Figure 11
400
500
216
W.L. Briggs, T. Sarie / Finite differmce
the problem by the Leapfrog method. Between time steps m=O and m=lO, the packet moves with a group velocity of C 2 5 which agrees well with the predicted discrete group velocity. In this case we see that the physical mode of the difference scheme actually transports energy in the wrong direction. These experiments are illustrative, but not surprising when one understands the dispersive nature of the difference scheme. A slightly more subtle effect is exposed in the next experiment. A single wave packet formed by a Gaussian envelope over a wave with wavelength 8h is taken as an initial condition, i.e. u(x,O) = e-8(x-5)2sinPx for 0 < x < 10. In this case K = 500 grid points have been-used with b=l, U=l and k = 2.4 x 10m6. Initially the packet is centered on the grid. Figure 11 shows the wave packet m=ZOO time steps later, having propagated to the left at the group velocity predicted by the discrete dispersion relation. The figure also includes the exact solution to the problem propagating to the left at the group velocity of the continuous dispersion relation, but leading the numerical solution. Notice that the numerical solution also consists of a smaller packet moving to the right at the same speed with which the initial packet moves to the left. This is the computational mode with a group velocity equal and opposite to the group velocity of the physical mode. This raises the important question of how the comoutational mode qets introduced into a calculation in the first place. The following argument shows how it happens and how difficult it-is to avoid introducing the computational It was seen in the previous section that mode. the general solution of the difference equations arising from the Leapfrog method is i(pnh-m@,(p)) i(pnh-m@,(p)) ‘ I u(m,n) = c A e +Bpe (3.1) P p where I$ = w,k cal mode
is the phase shift of the physi-
in one time step and $2 = w2k = TI - $I,
is the phase shift of the computational mode in The frequencies W, and ~2 are one time step. given by the discrete dispersion relation The sums range over all wave numbersofthe crete system
(0 5 p 5 $).
(Throughout
(2.5). dis-
f(x) = e'p', 0 5 P 5 5
used as an initial
condition
.
i(nPh-6) second initial condition as u(l,n) = e the initial profile shifted by a phase 6. If ’ these two initial conditions are substituted into the general solution (3.1), the coefficients Ap and BP can be determined. We find Ap=Bp=O for p # P and Ap + Bu = 1
*ve -i@,
Solving
BPe
i$,
=e
-ii5
gives
-i (69, ) Au = 1+e_2ie l+e
.
’
and B = lJ
l-e l+e
-i(Wl) .
2+i,
(3.2
It is evident that only if the second initial condition has a phase shift equal to the phase shift of the physical mode of wave number-u (i.e. 6 = $,) then A = 1, B = 0 and the computational mode will be sup:ressed. If the exact solution is used to obtain the second initial condition (a common practice in many numerical experiments), the computation mode is introduced imnediately since @l(u) # w(u)k (numerical phase shift # exact phase shift). One might hope that by specifying the exact solution at the first time step (m=l) and thereby introducing the computational mode, it might be possible to reproduce the exact solution for all times as a combination of the physical and computational mode. That this does not happen can be seen by looking at the evolution of the discrete solution. As before, let w be the frequency of the exact solution (mode P ) given by the continuous dispersion relation. Let o1 be the frequency
of the physical
mode and w2 be the
frequency of the computational mode. Assume that the second initial condition is given by i(pnh-ok) the exact solution: u(l,n) = e . Then for m 2 1 we find that i(Pnh-w,mk) i(Pnh-w2mk) u(m,n) = Alle + Bpe = ei(M-wmk)lA
e lJ
i(w-w,)mk
i(w-w2)mk +BPe
)
this arg-
ument, we assume that the appropriate conjugate modes are included to give a real solution). The continuous problem is accompanied by an initial condition u(x,O) = f(x). For now, consider an initial profile consisting of a single mode:
solutions of dispersive PDEs
This will be
for the difference
The difference scheme scheme: u(O,x) = elunh. also reouires a second initial condition which Since the continuous solumust be'specified. tion after one time step is the initial profile shifted by a phase w(p)k, we will give the
where Al, and Bu are given by (3.2). The term Q(m) represents the total modification of the exact solution after m time steps. It is important to notice that Q(m) amounts to more than a simple phase shift of the exact solV Although the scheme is neutral (nonution. dissipative), the presence of the computational mode (Bu # 0) makes the amplitude of the wave This re-iterates the fact train non-constant. that using the exact solution or a very good approximation to it introduces the computational mode and may not improve the discrete solution as an approximation to the exact solution.
W.L. Bnggs, T. Sarie / Finite difference solutions of dispersive PDEs
This discussion has concerned the situation in which the initial condition is a single spatial If the initial condition is a superposimode. tion of modes of the discrete problem, the above analysis pertains to each mode. This means that the computational mode of a particular spatial mode will be introduced if that mode does not have the phase shift given by the discrete dispersion relation at the first time step (m=l). If the initial condition cannot be represented by the spatial modes of the discrete problem, this additional error will inevitably excite the computational mode. Returning to the experiment of Figure 11, the presence of the computational mode can be explained. The initial condition cannot be represented exactly by the discrete spatial modes. Furthermore, the second initial condition for the Leapfrog method was given by the exact solution to the problem which, as we have seen, guarantees the appearance of the computational mode. 5.
A LINEAR ALIASING
We now investigate how the Leapfrog scheme will handle this initial profile. First note that
which
cosgx
where
of linear modes.
K is the number
of grid points,
w,(B+a)mkl + (*I + ei[(O-cdnh + w,(8-ahkl + (*) +
where wl(f3Aa) is the frequency
u(m,n)
Series
about
the wave
; ei(Bnh+wl(@)mk)
+ ei(Bnh+wl(8)mk)
solution
traveling
is a wave
train
at phase speed c(8)
by an envelope group velocity
of wave number a traveling C(B) = w;(B).
This case can now be compared in which
B < !$ ci < !$ but
at
to the situation
B +c,;K
2
.
Since
wave number B + c1 can no longer be resolved the grid, it is misrepresented by the wave number r where 6 + c1 = K - r since
by
ei(O+a)nh ,,i(K-r)nh,,i2neTirnh=,-irnh The two wave numbers seen by the difference scheme are now r and B - ~1. We define the average wave number of r and 8 - a to be
of the physical
number
l3 to
eic((nh+wi(6)mk) + (*)
eia(nh-wi(B)mk)
Then for some s, it follows r=i+s The relationship shown below.
The solution given by
+ (*)
among
these wave
of the difference
i[rnh-wl(r)mkl
,i[(B-a)nh where r.
that
B-cl=;-s.
-
Introducing
numbers
equations
is
is
+ (*) +
wl(B-ahkl
+ (*)
l? + a has been replaced
by the alias mode
ii and s, we have
u(m n) _ ,i[(ii+s)nh - wl(!i+s)mkl
,
ei[(bs?nh
mode associated with wave numbers @ +a given by (2.5). We assume no computational mode is present. Assuming @XX and using a standard group velocity argument, we expand W,(B + KY) in Taylor obtain
This
If
then the discrete solution will be a train of packets moving with the numerical group velocity C(o) given by (2.7). This can be seen as follows. The discrete solution is given by u(m,n)=ei[(6+a)nh
mk)]cos[a(nh-wi(B)mk]
where the approximate equality reflects the truncation of the Taylor Series for w,(B +a).
u(m,n)=e
= ei(0+a)x+(*)+ei(6-c)x+(*)
is a superposition
l? + cx ~5,
+ -J--
ERROR
Aliasing or the misrepresentation of high wave number components is usually associated with nonlinear problems in which interaction of modes produces higher harmonics. An interesting analogous effect can occur in a discrete representation of propagating wave packets. Consider an initial condition of the form u(x,O) = 4coscrx cosax where B >XX. This corresponds to a train of wave -packets which will be propagated by the continuous differential equation at the group velocity C(B) given by (2.4).
u(x,0)=4cos~x
$9
=4cos[B(nh
277
- w,(khkl
+ (*) +
+ (*).
Noting that ii >> s, we can repeat the same group velocity argument and expand wl(p + s) about After collecting terms and simwave number i. plifying, the discrete solution appears as u(m,n)=4cos[i(nh
$9 -yk)]cos[s(nh-wi(G)mk)].
This can be interprEted as a wave train with wave number lJ (short wave length) moving with w,(i;) phase speed c(li) = -modulated by an envel11 ope with wave number s (long waves) moving with
278
W.L. Briggs, T. Sarie / Finite difference solutions of dispersiue PDEs group velocity
C(i) = wi(i).
The two wave num-
bers of the discrete solution i and s are entirely spurious in that they appear nowhere in the initial conditions. They are introduced due to aliasing of the B + a mode. On final curiosity occurs
if B = !$ in which
case the short wave
within the packet has wavelength 2h. In this case s=O and the solution of the difference equations becomes u(m,n) = 4cos[$ih
- ol(u)mk]
which is a traveling monochromatic wave. In this case the difference scheme does not even "see" the wave packets of the initial condition and produces a uniform wave train. It might be argued that these effects are anomalies which are confined to the high wave number end of the spectrum. However if the energy contained even in the medium wave numbers is of any physical significance, it seems evident that this energy and its propagation can be grossly misrepresented by the difference scheme. Furthermore, these is no obvious argument to assert that similar behavior could not occur in nonlinear problems in which there s a mechanism to transfer energy originally in the low wave numbers to the high wave numbers. REFERENCES:
iI11
Crowley, vection, 484.
PI
Grotjahn, R. and O'Brien, J.J., Some inaccuracies in finite differencing hyperbolic equations, Mon. Wea. Rev.,104 (1976) pp. 180-194.
c31
Grotiahn, R.. Some consequences of finite differencing~revealed by 'group velocity errors, Beitrage zur Physik der Atmosphare, 50 (1977) pp. 231-238.
[41
A method for reducing disperFrom, J.E., sion in convective difference schemes, J. Comp. Phys. 3 (1968) pp. 176-189.
151
Roberts, K.V., and Weiss, difference schemes, Math. pp. 272-299.
[‘31
Stone, H.L. and Brian, P.L.T., Numerical solution of convective transport problems, A.1.Ch.E. Journal, g (1963) pp. 682-688.
[71
Trefethen, L., Group velocity in finite difference schemes, SIAM Review, 24 (1982) pp. 113-136.
181
Trefethen, L., Wave propagation and stability for finite difference schemes, Ph.D. Thesis, Stanford University, 1982.
W.P., Second order numerical adJ. Comp. Phys. 1 (1967) pp. 471-
N.O., Convective Comp. 20 (1966)
[9]
Vichnevetsky, R., Tomalesky, A.W., Spurious wave phenomena in numerical approximations of hyperbolic equations, Proceedings of the Fifth Annual Princeton Conf. on Information Sciences and Systems, 1971.
[lo] Vichnevetsky, R., Peiffer, B., Error waves in finite element and finite difference methods for hyperbolic equations, Advances in Computer Methods for Partial Differential Equations, AICA 1975.
Mathematics
FINITE DIFFERENCE EQUATIONS
SOLUTIONS
and Computers
OF DISPERSIVE
in Simulation
XXV (1983) 268-278 North-Holland
PARTIAL DIFFERENTIAL
William L. BRIGGS Department
of Mathematics
and Computer Science, Clarkson College of Technology, Potsdam,
NY 13676, U.S.A.
Talib SARIE Department
of Mathematics,
Jordan University, Amman,
Jordan
The second order leapfrog method is used to discretize the linearized KdV equation which is itself a dispersive partial differential equation. The resulting difference equation is solved and analyzed in terms of its dispersion relation and propagation properties. Numerical experiments are included to illustrate some of the results obtained. Finally, a novel aliasing error which affects the propagation of wave packets is presented
1.
INTRODUCTION
It might be thought that the analysis of finite difference schemes with constant coefficients is Inan exercise which has long been exhausted. deed the earliest investigations, which were concerned primarily with the stability of Cerschemes, date back at least thirty years. tainly at that time it was recognized that even stable finite difference schemes, applied to time dependent problems, necessarily introduce errors in resolving the amplitude and phase of The orooerties of these damoinq and solutions. dispersion errors has been the subject of &IStained research ever since. Despite this level of activity, only recently has the role of group velocity in the analysis of difference schemes been recognized. This study is undertaken in the belief that, while the analysis of difference schemes is an area that has been wellcovered, there are still some details that could use further illumination.
In the followino oaoes, we will be concerned with the numerica! solution of a dispersive partial differential equation. Many of the conclusions and results would apply to any equation However the choice of a disperof propagation. sive equation does provide some new effects and The importance of propagation properinsights. ties (e.g. dispersion and group velocity) in understanding the behavior of difference schemes seems under-emphasized in much of the numerical Although some of the analysis literature. results which follow are not new or unexpected, they do appear in a new setting and it serves at least a pedagogic purpose to present them again. Much of the work which directly precedes and anticipates this present study-was done to develop schemes for fluid flow problems. Among this early work, the papers which deal specifically with phase error and its control seem In the 1963 work of most appropriate here. 0378-4754/83/$3.00 0 1983, IMACS/Elsevier
Science Publishers
Stone and Brian [6], a family of schemes for a one-dimensional convection-diffusion problem was presented. Although the terminology was not present, the dispersion relation for the schemes was derived and an optimal scheme which minimizes damping and phase errors was developed. In 1966, Roberts and Weiss [5] pointed out the role of damping and dispersion in pure advection problems and developed a fourth order scheme (in space) which considerably reduces phase error. Crowley [l] compared the damping and dispersion properties of several existing schemes, particularly for advection problems in two dimensions. At about the same time, Fromm [4] analyzed split time step methods and devised his well-known method of zero average phase error, again for two dimensional advection problems. All of these studies clearly recognized the fact that while damoinq errors can be avoided by usinq time-centered (leapfrog) difference schemes; all difference schemes are dispersive. That is, when any finite difference scheme is applied to an equation of propagation, each spatial mode of the difference scheme has a phase speed which depends upon the wavelength of the mode. This is true whether the original partial differential equation has this property or not. The earliest recognition that group velocity plays a role in the analysis of finite difference schemes seems to have occurred in 1971 when Vichnevetsky and his colleagues [9,10] considered the one-dimensional advection equation and examined semi-discrete (time continuous) Using Fourier analysis, propagation schemes. properties of these schemes were derived, the important notions of group velocity and energy propagation of the finite difference scheme were introduced. and the role of soatial comoutational modes was illustrated through numerical experiments. Shortlv thereafter Grotjahn [2,3] studied group velocity in finite difference schemes arising from one and two dimensional
B.V. (North-Holland)
269
W.L. Briggs, T. Sarie / Finite difference solutions of dispersive PDEs
More recently, Trefethen geophysical problems. [7,8] has studied the one-dimensional advection equation modeled by fully discrete schemes. He derived the relevant dispersion relations and in a series of imaginative numerical experiments, clearly showed the behavior of the computational mode, pointinq out that a fully discrete scheme can support computational modes in both space and time. In addition, known stability criteria were interpreted in terms of the group velocity. In the present order leapfrog equation
study, we consider the second method applied to the dispersive
ut + Uux + buxxx = 0.
U,be
R
This can be regarded as the linear analog of the Korteweg-deVries (KdV) equation. A brief discussion in the next section presents the propagation properties of this equation. Most important for our purposes is the fact that the dispersion relation, which relates the frequency w of a spatial mode to its wave number u, is a nonlinear function of u. This coupled with the fact that the difference scheme has its own nonlinear dispersion relation means that discrepancies between the propagation properties of the continuous and discrete problems are even greater than with a non-dispersive partial differFor example, the continuous ential equation. problem can support wave trains with zero phase speed and wave packets with zero group velocity. We will see that the difference scheme approximates these stationary modes, but in the process, produces spurious stationary modes of its own.
2.
THE CONTINUOUS
In this section we will investigate finite difference,solutions to the dispersive partial differential equation ut + lJux + bu
In Section 4, several numerical experiments are performed to illustrate some of the ideas of the previous section. The presence of computational modes in a calculation is explored and an argument is given to show that the choice of initial conditions at the first two time levels ultimately determines whether computational modes appear. The paper concludes by pointing out an aliasing error that can occur when wave packets are propagated by finite difference Whether this error is a serious danger schemes. or merely a curiosity has yet to be determined.
xxx
= 0,
U,b
E
R
(2.1)
on the intervals 0 < x < ZIT, t > 0 subject to periodic boundary conditions and an initial condition u(x,O) = f(x). Without loss of generality it is assumed that U > 0, however consideration b > 0 and b < 0 leads to some interesting situations. A brief look at the propagation properties of continuous solutions to this problem will provide some background and also a point of comparison for the discrete solutions. Equation (2.1) admits normal mode solutions of the form u(x,t) = ei(ux-mt)
where the wave number u and the frequency w are related through the dispersion relation w = Uu - bu3 = -bu(u'
- ;).
(2.2)
This expression, in which u assumes integer values, gives the frequency w with which any The graphs of the dispersion wave travels. relation for the cases b > 0 and b < 0 are shown Note that in in Figures la and lb respectively. both cases w is an odd function of u, a property which will be preserved by the analogous discrete dispersion relations. The phase speed or speed of propagation of a pure mode with wave number u is given by c(u) = ti
In the third section of the paper, the difference equations arising from the Leapfrog scheme are solved subject to periodic boundary conditions using separation of variables. This approach seems neglected, but it has the advantage that it produces a general solution in terms of corresponding spatial and temporal modes and also gives stability conditions and the dispersion relation in a straightforward way. The phase speed and group velocity curves for the difference scheme are presented and discussed, particularly in relation to the corresponding curves for the continuous problem. The stability condition for the scheme proves to have some interesting structure which is also discussed.
PROBLEM
= U - bu2.
(2.3)
The graphs o/the phase speed for b > 0 and b < 0 are shown in Figures 2a and 2b respectiveIn both cases the phase speed is a nonly. constant function of p. which is to say that different modes travel at different speeds. Of special interest is the case b > 0 in which the phase speed is of both signs with a stationary mode at u= Jms. Finally, the group velocity or the speed at which a packet of waves of wave number u propagates is given by C(u) = 4
= U - 3bu2.
( 2 :4)
Again the group velocity is a non-trivial function of the wave number (Figures 3a and 3b) which gives rise to the term dispersive in describing this problem. In the case b > 0, the group velocity is of both signs. Energy in the long waves (small u) propagates in the positive direction. Energy in the p = Jm56 mode is stationary and energy in the short waves (large Both u) propagates in the negative direction. phase speed and group velocity are even functions of the wave number. 3.
ANALYSIS
OF THE DIFFERENCE
a.
Solution
by separation
EQUATIONS
of variables.
With this understanding of the continuous problem, we now turn to finite difference approxima-
270
W.L. Briggs, T. Sarie / Finite difference solutions of dispersive PDEs
tions to (2.1).
A uniform
grid of width
h = 'K1_
is established on the interval 0 < x < Zn, while the interval t > 0 is partitioned-at intervals of length k. In this study we shall be concerned only with finite difference representations of equation (2.1). The Leapfrog method is obtained by replacing all derivatives in (2.1) by second order finite differences. Letting u(m,n) be the discrete approximation to the exact solution at the point x = nh, t = mk, the following difference equations result. u(m+l,n)
- u(m-1,n)
+ a[u(m,n+l)
- u(m,n-1)]
+ B[u(m,n+2)-2u(m,n+l)+2u(m,n-1)-u(m,n-2)]=0 for 0 5 n 5 K-l, m > 1 with f(x,,) and u(l,n)
specified.
u(m,O)=u(m,K),u(O,n)= The two parameters
a = !$ and @ = k will have important parts to h3 play in all that follows. There are several approaches to solving linear, constant coefficient partial difference equations. We proceed by separation of variables. Letting u(m,n) = M(m)N(n) it follows that M(m+l) Y
1C
lb
20
(a)
25
- M(m-1)
= XM(m) and
N(n+2)+(a-2B)N(n+l)+~N(n)-(a-24)N(n-l)-N(n-2) -0 B 5 B where A is the separation constant which will be determined. Turning to the N equation, we seek solutions of the form N(n) = g". The periodic boundary conditions of the original equation appear as the condition N(0) = N(K). This requires that < be of the form
ii!??%! c=e
K
JJ = 0, -+ 1, _ + 2, . . . . $
.
We see immediately at this point that the spatial part of the discrete solution is quantized and consists
Figure
oispersion la)
relation
for
b = + 1.52 I lo-'
the
is constant
n = Fmode
has a wavelength
modes.
The
over the grid while
the
of 2h‘. This is, of
course, in contrast to the continuous problem in which modes of all wavenumbers are available For the discrete problem, only to the system. those initial conditions which are a linear com-
1 co”tin”o”s (b)
of only $ distinct
n = 0 mode
problem,
b = -1.52 x 10
” = 1 -2
bination
of these t modes
can be represented
An initial condition which lies outexactly. side of this space will not be represented or propagated without error. We can then express the spatial part of the solution as a linear combination of modes of the form N(,,)=e
?_.
K =elpnh
= 0, -+ 1, _ + 2,..., t;
.
Substituting this into the original equation for N(n) allows the separation constant X to be determined for each spatial mode. We find that A(p) = -Zi[(a-2B)sinuh
+ Bin
2uh].
271
W. L. Briggs, T. Sarie / Finite difference solutions of dispersme PDEs Using these values assuming
of X in the M(m) equation
a solution x&z
v+(u) -
= -is(u) +
= 2
and
of the form M(m) = ym gives l-g(u)
2
The fact where g(u) = (a-2a)sinuh + @sin 2ph. that there are two values of y associated with each mode arises, of course, from the necessity of solving a second order difference equation to determine 14. This will have important implications in all that follows. The temporal
Y+(V) which
behavior of mode u is determined is turn depends upon g(u). Three
cases
must be considered.
(i)
If lg(lJ) > 1 for then IY u) 1 > 1
by
at least one value of u. and hence solutions can qrow with increasing m (time). These are instable solutions which are not consistent with the continuous solutions.
I
(ii)
If ig(u)I = 1 for at least one value of u, then y, = y_ = i or y, = y_ = -i. In either case the difference equation for M has a double root and mode u has temporal behavior
of the form M(m) = clim + c2mim.
This slow linear growth, while not catastropic, is also inconsistent with the continuous solutions and must be avoided. (iii) iG+lz\l;'l Y Bn~~lf'lo~~~"~~v~'a~~o~~~~nt The conamFlitude, oscillatory behavior. dition ig(p)I < 1 imposes restrictions on the parameters CL and B which amount to stability conditions for each mode and for the entire scheme. We will return to the precise form of these stability conditions shortly. Clearly a useful difference scheme results only in this last case. b.
Discrete
dispersion
If we assume that (g(p)1 < 1 for all wave numbers u, then the temporal part of the solution for mode p takes the form
M(m) = cle
-iw2mk +ce
2
where o,k = sin- '{(a-2a)sinuh 0 (~,k
~II
and
+ Bsin 2u.h)
IJ = 0, + 1,
+
(2.5)
w2k = n-w,k.
The full solution of (2.1) associated u finally appears as i(unh-wlmk) i(unh-m2mk) u(m,n)=c,e +ce 2 2,...,
+- . 2
Rewriting
with mode
(2.5) in the form
sinwk = (a-2b)sinph
+ f3sin2ph
(2.6)
gives the dispersion relation for the difference As shown in Figure 0, this relation equations. is double valued with respect to both w and u and hence both w, and w2 are represented. The figure also shows the continuous dispersion relation. For small wave numbers, the continuous dispersion relation is matched very closely by the w, branch of the discrete dispersion relation.
For this reason,
the o, component
of the
discrete solution is called the physical mode since it corresponds and converges to the continous solution. The w2 component of the discrete solution is entirely spurious and is called the computational mode. To be completely accurate, these modes should be called the physical mode and the computational mode in time since they are two temporal modes corresponding to the same wave number. However, the dispersion relation could be read in the other direction also. If the frequency is fixed (say be time forcing at a boundary)then two spatial modes arise, one physical and one computational. Therefore for one (u,w) pair in the physical problem, the difference scheme is capable of producing four (u,w) pairs, only one of which approximates the physical pair. c.
Phase speed and group
velocity.
The propagation properties of various discrete modes may be obtained from the dispersion relation (2.5) or (2.6). The phase speed is given by c(u) = &
relation.
-iw,mk
These may be regarded as the normal modes of Any linear combination of the discrete problem. these modes can be represented by the difference scheme.
sin-'[(a-26)sinph
+ Bsin2uh].
Expanding this expression for small wave numbers (uh << 1) we find that for the physical mode c(,,) = "_bu2 _ +
(uh)2 + O(uh)4.
The difference scheme gives a second order approximation to the phase speed of long waves. Figure 5 is a graph of the phase speed of the physical mode as it varies with u. The phase speed curve for the corresponding continuous problem is also shown for comparison. As is well known, any difference scheme applied to a propagation equation will be dispersive. In this study, in which the continuous problem is also dispersive, the departures of the numerical phase speed from the continuous phase speed are even more pronounced. Of particular interest are the stationary modes. For b > 0, the continuous phase speed has a zero at n = 0 and uO = + m in the u-m plane
W. L. Briggs, T. Sarie / Finiie difference solutions
212
mare speeafor (a)
b-
the
+ 1.52 x 10
(bl
(bl
Figure 2
Figure 3
contin”o”s -2
prm,en, (b)
”
=
(a)
c
2o
IO _
for the b - t1.52 x 10-2
Group velocity
1
b = -1.52 x IO-'
2-
-I
of dispersive PDEs
continuousproblem. 0 = 1 (b)
b = -1.52 x 10-2
W.L. Briggs, T. Sarie / Finite difference solutions (Both the
yields
numerical and continuous phase speeds are even functions of u and negative wave numbers can be omitted from the discussion.) Now if Usb%O(l),
u0
oratuOh
= +m
in the ph-mk
plane.
then u0 % O(1) << $ and the stationary
.
This puts the
stationary wave near the short wave end of the Assuming that U and b are discrete spectrum. given and assuming that the zero of the phase speed is of physical importance, then h and k should be chosen so that u0 lies near the center of the discrete
spectrum.
Of course the difference scheme will not reproduce the stationary mode of the continuous problem even if it is centered in the discrete The discrete phase speed has zeros spectrum. at u=O, u = k and at
"
0(h2).
The above remarks pertain to the phase speed of By the symmetry of the the physical mode. numerical dispersion relation, it can be seen that the phase speed of the computational mode is equal in magnitude but opposite in sign to that of the physical mode. In many physical problems and numerical experiments the group velocity may be of more interReturning to the est than the phase speed. discrete dispersion relation the numerical group velocity may be computed to be C(u) = ;
(a-2!3)cosnh + 2BcosZiJh 1-[(a-2B)sinuh
=
3 + $ h2
+ 4sig
U-3bp2 - $(‘-a2)(ph2) + O(j&,)4_
Again the group velocity of the continuous lem is reproduced by the physical mode to
+ 0(h2)
The first value of u3
(2.7) prob-
The group velocity of the computational O(uh)2. mode is equal in magnitude but of opposite sign Figure 6 is a to that of the physical mode. graph of the exact group velocity and the numerical group velocity as it varies with u. As with phase speed, the group velocity of the continuous problem has much more structure than it For b > 0 would in a non-dispersive problem. the mode p. = Jm36 has zero group velocity and
corresponds
to the sta-
tionary mode of the continuous problem, while the second value is a spurious stationary mode at a higher wave number. The difference scheme clearly has the capability of moving energy in the wrong direction even in the intermediate wave numbers. When the computational mode is introduced the situation becomes even more comWe will return to these results in plicated. later numerical experiments. d.
Stability
conditions.
We now return to the question of stability and investigate the conditions which determine a stable difference scheme. It was seen that if /g(p)1 = I(a-2g)sinph for all wave numbers
+
273
$ + O(h2)
2=
wave may
not be resolved within the spectrum of the If U > b then the situation difference scheme. is even worse with the stationary wave longer than the fundamental wave of the difference At the other extreme, if U s O(1) and scheme. b 'irh*U then a%fi and 'JO Q $
of dispersive PDEs
+ BsinPuhl
< 1
0 2 u 5 !j then solutions
of
the difference scheme remain bounded as m (time) These inequalities generate regions increases. of stability in the (@,a) plane which are bounded by the pairs of parallel lines (a-2e)sinph
- BsinPuh
= + 1
or a = 2(1-cosph)b
+ cscuh
(2.8)
where ~1 > 0 and -m < B cm. Figure 7 shows the stability strips of several of the higher and more critical modes. For a given choice of o. and B to yield a stable scheme the point (B,a) must lie in the intersection of the stability regions of all modes. Looking at the slopes and intercepts of the boundary lines in (2.8) it is evident that the stability regions enlarge with decreasing wave number. It is the higher modes (uh = F,
z, 5, $) that determine
of the stability The boundaries given below:
the boundaries
region.
for these critical
modes
are
ZCc('+2_+3B 3 -J5 E:a=t,t2fi 2 -
as in the case of phase speed, h and k must be chosen carefully in order to resolve VU within
;:o=lZ+@
the spectrum of the difference scheme. However is resolved the difference scheme only appr even if "8 ximates it and furthermore produces a The zeros second stationary mode of its own. of the discrete group velocity are not simply expressed, however expanding in powers of h
;
2
: a = + VT +
(e-n)@.
The outline of the common region of stability of these modes is shown in the figure. The stability strips for all other modes contain this comnon region. We see that if u 2 1, then
W. L. Briggs, T. Sarie / Finite difference solutions
214
k
5-
-
a-
3-
Z-
I-
L) 5.8
h 5.2
6.6
0.0
I.6
of dispersive
PDEs
215
W. L. Briggs, T. Sarie / Finiie difference solutions of dispersive PDEs
then a stable 161 >$
scheme
is not possible,
= .39, then stability
while
if
is not possible.
A different picture of the stability conditions is obtained if the conditions Ig(p)I < 1 are expressed in terms of k and h. Some calculation eventually leads to the stability condition k<
h3
where
? ? UsinnhlRL-hLI
a2 =
2b l-cosph) (,, "
If vh is fixed to consider for the case b > 0. the stability of a single mode, then the stability region in the (h,k) plane will have the general appearance shown in Figure 8. The location of the asymptote at h = R depends upon the particular mode under consideration and on the relative values of the parameters b and U. For longer (more stable) waves the asymptote It is also interesting that shifts to the left. the choice of h = R, which would allow arbitrarily large time steps, is the combination of parameters for which the mode has zero phase The stability region for a single mode speed. in the case b < 0 is given by k<._
h3
1
0
* A I’
-1
Usinnh(R2
+ h2)
*t
.
II
-11
II
’
Again, the details of the curve will depend upon the choice of b and U and upon the mode under A general curve for this case is consideration. shown in Figure 8.
4.
EXPERIMENTS
AND CONSEQUENCES
In this section several numerical experiments are presented to illustrate the ideas of the In addition some examples are previous section. given to point out further consequences of using a discrete model to represent a dispersive partial differential equation.
"=l.b
- 1.52 x 10-2, k - 8.7 x 10-3
(a)
Initial cmditions.
(b)
u(nh. 1Ok)
"("h.0) = cor(o"h)cos(mh).
(1 = 2. B - 30
As in previous studies, we present some wave packet experiments to demonstrate the propagaIn the tion properties of the Leapfrog method. first two experiments, we use a laboratory with K=72 grid points, U=l, b=l.52x 10e2, k=8.7x10m3. In the first case, we consider an initial condition consisting of a wave with wave number B=24 (wavelength 3h) modulated by a long wave with wave number a=1 (wavelength 72h), i.e. This forms a series of u(x,O) = cosux cosox. wave packets which should propagate at the group velcoity C(24) = -25.5 according to the group velocity relation (2.4). Figure 9 shows the discrete solution at m=O and m=lO time steps. Clearly the packet is almost stationary which is consistent with the discrete group velocity (2.7) which qives C(24) s 0. In the second case.initial-conditions of the form u(x,O) = The cosax COSBX with a=2 and 0=30 are taken. continuous equation should propagate these wave packets with a group velocity of C(30) = Figure 10 shows the result of solving -41.0.
I -Nunrlcal
solution
0
-1
n
-2 0
100
200
300 Figure 11
400
500
216
W.L. Briggs, T. Sarie / Finite differmce
the problem by the Leapfrog method. Between time steps m=O and m=lO, the packet moves with a group velocity of C 2 5 which agrees well with the predicted discrete group velocity. In this case we see that the physical mode of the difference scheme actually transports energy in the wrong direction. These experiments are illustrative, but not surprising when one understands the dispersive nature of the difference scheme. A slightly more subtle effect is exposed in the next experiment. A single wave packet formed by a Gaussian envelope over a wave with wavelength 8h is taken as an initial condition, i.e. u(x,O) = e-8(x-5)2sinPx for 0 < x < 10. In this case K = 500 grid points have been-used with b=l, U=l and k = 2.4 x 10m6. Initially the packet is centered on the grid. Figure 11 shows the wave packet m=ZOO time steps later, having propagated to the left at the group velocity predicted by the discrete dispersion relation. The figure also includes the exact solution to the problem propagating to the left at the group velocity of the continuous dispersion relation, but leading the numerical solution. Notice that the numerical solution also consists of a smaller packet moving to the right at the same speed with which the initial packet moves to the left. This is the computational mode with a group velocity equal and opposite to the group velocity of the physical mode. This raises the important question of how the comoutational mode qets introduced into a calculation in the first place. The following argument shows how it happens and how difficult it-is to avoid introducing the computational It was seen in the previous section that mode. the general solution of the difference equations arising from the Leapfrog method is i(pnh-m@,(p)) i(pnh-m@,(p)) ‘ I u(m,n) = c A e +Bpe (3.1) P p where I$ = w,k cal mode
is the phase shift of the physi-
in one time step and $2 = w2k = TI - $I,
is the phase shift of the computational mode in The frequencies W, and ~2 are one time step. given by the discrete dispersion relation The sums range over all wave numbersofthe crete system
(0 5 p 5 $).
(Throughout
(2.5). dis-
f(x) = e'p', 0 5 P 5 5
used as an initial
condition
.
i(nPh-6) second initial condition as u(l,n) = e the initial profile shifted by a phase 6. If ’ these two initial conditions are substituted into the general solution (3.1), the coefficients Ap and BP can be determined. We find Ap=Bp=O for p # P and Ap + Bu = 1
*ve -i@,
Solving
BPe
i$,
=e
-ii5
gives
-i (69, ) Au = 1+e_2ie l+e
.
’
and B = lJ
l-e l+e
-i(Wl) .
2+i,
(3.2
It is evident that only if the second initial condition has a phase shift equal to the phase shift of the physical mode of wave number-u (i.e. 6 = $,) then A = 1, B = 0 and the computational mode will be sup:ressed. If the exact solution is used to obtain the second initial condition (a common practice in many numerical experiments), the computation mode is introduced imnediately since @l(u) # w(u)k (numerical phase shift # exact phase shift). One might hope that by specifying the exact solution at the first time step (m=l) and thereby introducing the computational mode, it might be possible to reproduce the exact solution for all times as a combination of the physical and computational mode. That this does not happen can be seen by looking at the evolution of the discrete solution. As before, let w be the frequency of the exact solution (mode P ) given by the continuous dispersion relation. Let o1 be the frequency
of the physical
mode and w2 be the
frequency of the computational mode. Assume that the second initial condition is given by i(pnh-ok) the exact solution: u(l,n) = e . Then for m 2 1 we find that i(Pnh-w,mk) i(Pnh-w2mk) u(m,n) = Alle + Bpe = ei(M-wmk)lA
e lJ
i(w-w,)mk
i(w-w2)mk +BPe
)
this arg-
ument, we assume that the appropriate conjugate modes are included to give a real solution). The continuous problem is accompanied by an initial condition u(x,O) = f(x). For now, consider an initial profile consisting of a single mode:
solutions of dispersive PDEs
This will be
for the difference
The difference scheme scheme: u(O,x) = elunh. also reouires a second initial condition which Since the continuous solumust be'specified. tion after one time step is the initial profile shifted by a phase w(p)k, we will give the
where Al, and Bu are given by (3.2). The term Q(m) represents the total modification of the exact solution after m time steps. It is important to notice that Q(m) amounts to more than a simple phase shift of the exact solV Although the scheme is neutral (nonution. dissipative), the presence of the computational mode (Bu # 0) makes the amplitude of the wave This re-iterates the fact train non-constant. that using the exact solution or a very good approximation to it introduces the computational mode and may not improve the discrete solution as an approximation to the exact solution.
W.L. Bnggs, T. Sarie / Finite difference solutions of dispersive PDEs
This discussion has concerned the situation in which the initial condition is a single spatial If the initial condition is a superposimode. tion of modes of the discrete problem, the above analysis pertains to each mode. This means that the computational mode of a particular spatial mode will be introduced if that mode does not have the phase shift given by the discrete dispersion relation at the first time step (m=l). If the initial condition cannot be represented by the spatial modes of the discrete problem, this additional error will inevitably excite the computational mode. Returning to the experiment of Figure 11, the presence of the computational mode can be explained. The initial condition cannot be represented exactly by the discrete spatial modes. Furthermore, the second initial condition for the Leapfrog method was given by the exact solution to the problem which, as we have seen, guarantees the appearance of the computational mode. 5.
A LINEAR ALIASING
We now investigate how the Leapfrog scheme will handle this initial profile. First note that
which
cosgx
where
of linear modes.
K is the number
of grid points,
w,(B+a)mkl + (*I + ei[(O-cdnh + w,(8-ahkl + (*) +
where wl(f3Aa) is the frequency
u(m,n)
Series
about
the wave
; ei(Bnh+wl(@)mk)
+ ei(Bnh+wl(8)mk)
solution
traveling
is a wave
train
at phase speed c(8)
by an envelope group velocity
of wave number a traveling C(B) = w;(B).
This case can now be compared in which
B < !$ ci < !$ but
at
to the situation
B +c,;K
2
.
Since
wave number B + c1 can no longer be resolved the grid, it is misrepresented by the wave number r where 6 + c1 = K - r since
by
ei(O+a)nh ,,i(K-r)nh,,i2neTirnh=,-irnh The two wave numbers seen by the difference scheme are now r and B - ~1. We define the average wave number of r and 8 - a to be
of the physical
number
l3 to
eic((nh+wi(6)mk) + (*)
eia(nh-wi(B)mk)
Then for some s, it follows r=i+s The relationship shown below.
The solution given by
+ (*)
among
these wave
of the difference
i[rnh-wl(r)mkl
,i[(B-a)nh where r.
that
B-cl=;-s.
-
Introducing
numbers
equations
is
is
+ (*) +
wl(B-ahkl
+ (*)
l? + a has been replaced
by the alias mode
ii and s, we have
u(m n) _ ,i[(ii+s)nh - wl(!i+s)mkl
,
ei[(bs?nh
mode associated with wave numbers @ +a given by (2.5). We assume no computational mode is present. Assuming @XX and using a standard group velocity argument, we expand W,(B + KY) in Taylor obtain
This
If
then the discrete solution will be a train of packets moving with the numerical group velocity C(o) given by (2.7). This can be seen as follows. The discrete solution is given by u(m,n)=ei[(6+a)nh
mk)]cos[a(nh-wi(B)mk]
where the approximate equality reflects the truncation of the Taylor Series for w,(B +a).
u(m,n)=e
= ei(0+a)x+(*)+ei(6-c)x+(*)
is a superposition
l? + cx ~5,
+ -J--
ERROR
Aliasing or the misrepresentation of high wave number components is usually associated with nonlinear problems in which interaction of modes produces higher harmonics. An interesting analogous effect can occur in a discrete representation of propagating wave packets. Consider an initial condition of the form u(x,O) = 4coscrx cosax where B >XX. This corresponds to a train of wave -packets which will be propagated by the continuous differential equation at the group velocity C(B) given by (2.4).
u(x,0)=4cos~x
$9
=4cos[B(nh
277
- w,(khkl
+ (*) +
+ (*).
Noting that ii >> s, we can repeat the same group velocity argument and expand wl(p + s) about After collecting terms and simwave number i. plifying, the discrete solution appears as u(m,n)=4cos[i(nh
$9 -yk)]cos[s(nh-wi(G)mk)].
This can be interprEted as a wave train with wave number lJ (short wave length) moving with w,(i;) phase speed c(li) = -modulated by an envel11 ope with wave number s (long waves) moving with
278
W.L. Briggs, T. Sarie / Finite difference solutions of dispersiue PDEs group velocity
C(i) = wi(i).
The two wave num-
bers of the discrete solution i and s are entirely spurious in that they appear nowhere in the initial conditions. They are introduced due to aliasing of the B + a mode. On final curiosity occurs
if B = !$ in which
case the short wave
within the packet has wavelength 2h. In this case s=O and the solution of the difference equations becomes u(m,n) = 4cos[$ih
- ol(u)mk]
which is a traveling monochromatic wave. In this case the difference scheme does not even "see" the wave packets of the initial condition and produces a uniform wave train. It might be argued that these effects are anomalies which are confined to the high wave number end of the spectrum. However if the energy contained even in the medium wave numbers is of any physical significance, it seems evident that this energy and its propagation can be grossly misrepresented by the difference scheme. Furthermore, these is no obvious argument to assert that similar behavior could not occur in nonlinear problems in which there s a mechanism to transfer energy originally in the low wave numbers to the high wave numbers. REFERENCES:
iI11
Crowley, vection, 484.
PI
Grotjahn, R. and O'Brien, J.J., Some inaccuracies in finite differencing hyperbolic equations, Mon. Wea. Rev.,104 (1976) pp. 180-194.
c31
Grotiahn, R.. Some consequences of finite differencing~revealed by 'group velocity errors, Beitrage zur Physik der Atmosphare, 50 (1977) pp. 231-238.
[41
A method for reducing disperFrom, J.E., sion in convective difference schemes, J. Comp. Phys. 3 (1968) pp. 176-189.
151
Roberts, K.V., and Weiss, difference schemes, Math. pp. 272-299.
[‘31
Stone, H.L. and Brian, P.L.T., Numerical solution of convective transport problems, A.1.Ch.E. Journal, g (1963) pp. 682-688.
[71
Trefethen, L., Group velocity in finite difference schemes, SIAM Review, 24 (1982) pp. 113-136.
181
Trefethen, L., Wave propagation and stability for finite difference schemes, Ph.D. Thesis, Stanford University, 1982.
W.P., Second order numerical adJ. Comp. Phys. 1 (1967) pp. 471-
N.O., Convective Comp. 20 (1966)
[9]
Vichnevetsky, R., Tomalesky, A.W., Spurious wave phenomena in numerical approximations of hyperbolic equations, Proceedings of the Fifth Annual Princeton Conf. on Information Sciences and Systems, 1971.
[lo] Vichnevetsky, R., Peiffer, B., Error waves in finite element and finite difference methods for hyperbolic equations, Advances in Computer Methods for Partial Differential Equations, AICA 1975.