Origin of Numerical Dispersion
Let us consider the solution of the simple convective wave equation:
d U d U
+ c = 0 dt dx
























and t. They are as follows:
N
K(x, t) =——– ^2 ajU(x + jAx, t) (4.17a)
AX j=N
3
u(x, t + At) = u(x, t) + At ^2 bjK(x, t — jAt), (4.17b)
j=0
with initial condition
( Ф(x), 0 < t < At
u(x, t) = (4.17c)
I 0, t < 0.
By taking the FourierLaplace transform of problem (4.17) and using the Shifting Theorem (see Appendix A), it is straightforward to find, after some algebra, that the solution is
ІФ (a)
2п (ш — ca)
where U is the FourierLaplace transform of the solution of the finite difference equation (4.17). Eq. (4.18) is identical to the FourierLaplace transform of the solution of the original convective wave equation; i. e., Eq. (4.13), if ш and a in Eq. (4.18) are replaced by ш and a.
The exact solution of the finite difference solution of (4.11) is given by the inverse transform of Eq. (4.18) as follows:
The double integral may be evaluated by first calculating the contributions of the poles in the ш plane. The poles are given by the zeros of the denominator of the integrand, that is, the roots of
ш(ш) = ca(a). (4.20)
Eq. (4.20) is the dispersion relation. This is formally the same as the dispersion relation of the original partial differential equation; i. e., Eq. (4.14), provided ш is replaced by ш and a is replaced by a. For this reason, the finite difference scheme is referred to as dispersionrelationpreserving (DRP).
For a fourlevel time marching scheme, there would be four poles or Eq. (4.20) has four roots. Let the poles be ш = шj (a), j = 1, 2, 3, 4. One of the poles has Im(«) nearly equal to zero. This is the pole that corresponds to the physical solution. The other poles lead to solutions that are heavily damped.
Numerical dispersion can arise from spatial as well as temporal discretization. Numerical dispersion due to spatial discretization alone will be considered first. For this purpose, the limit At ^ 0 is applied, that is, time is treated exactly. This is
The integrand has a pole at ш = ca(a). On deforming the inverse contour Г to pick up this pole in the ш plane, it is easy to find through the use of the Residue Theorem
TO
u(x, t) = J Ф(a)ei(ax —ш) da. (4.23)
— то ш=са( a)
Now for large t (with x/t fixed), the aintegral of Eq. (4.23) may be evaluated by the method of stationary phase (see Appendix A). The stationary phase point a = as is given by the zero of the derivative of the phase function Ф = (ax/t — ca(a)) with respect to a. This leads to
which yields as as a function of x/t. The asymptotic solution is
where sgn(z) is the sign of z.
Now, ш and a are related by Eq. (4.22). In general, da/da is not a constant equal to 1, so that the group velocity varies with wave number. Different wave numbers of the initial disturbances will propagate with a slightly different speed leading to numerical dispersion.
In the more general case of solution (4.19), when both spatial and temporal discretization are used, the group velocity may be calculated by implicit differentiation of Eq. (4.20). This gives
da
dm = C da da dm
dш
Eq. (4.26) indicates that numerical dispersion can be the result of imperfect spatial discretization, i. e., da/da = 1, or imperfect temporal discretization, i. e., drn/drn = 1. Consideration of imperfect time discretization will be deferred to later because it is quite involved. The effect of dispersion due to spatial discretization will now be illustrated by a numerical example.
For the convective wave equation (4.11), Eq. (4.24) reveals that the variation in group velocity is caused by the variation of the slope of the a Ax versus a Ax curve. Figure 4.2 shows the da/da versus a Ax curves for the 7point and 15point
Figure 4.2. The da/da vs. aAx curves………. , sixthorder central difference scheme; 7point DRP scheme——– , 15point DRP scheme. 
stencil DRP finite difference schemes. Notice that the da/da curve for the 7point stencil DRP scheme has a peak at a Ax = 0.67 at which da/da = 1.00276. That is, the wave component with a Ax = 0.67 (wave with 9.38 mesh points per wavelength) will propagate faster than the exact wave speed da/da = 1.0 by about 0.3 percent. Suppose, in a computation, the waves propagate over a distance of 400 mesh points. In this case, when the main part of the wave packet reaches x = 400, the wave component with a Ax = 0.67 will reach x = 401.1. Under these circumstances, numerical dispersion would just be noticeable as wave dispersion exceeds one mesh point. On the other hand, if the standard sixthorder scheme is used instead, severe numerical dispersion will result when the main wave reaches x = 400. Waves in the wave number range of a Ax > 0.6, having a group velocity less than 1.0, would form trailing waves.
A relevant question to ask is how to reduce numerical dispersion due to spatial discretization. This can be done by using a scheme with a larger stencil. Figure 4.2 shows the da/da curve for the 15point stencil DRP scheme. For aAx < 1.73, the group velocity differs from 1.0 by no more than 0.1 percent. Thus, if the wave packet of the solution contains only waves with wave number aAx < 1.73, then there will be no observable numerical dispersion even after the wave packet propagates over a distance of up to 1000 mesh points. As a concrete example, Figure 4.3 shows the
Figure 4.3. Comparisons between computed solutions and exact solution for the onedimensional convective wave equation. (a) sixthorder central difference scheme; (b) 7point DRP scheme; (c) 15point DRP scheme. 
computed solution of the convective wave equation (4.11) with initial condition in the form of a Gaussian with a halfwidth of 3 mesh spacings,
u(x, 0) = e(ln2)()2,
using the standard sixthorder scheme, the 7point stencil, and the 15point stencil DRP scheme. It is easy to see that there is significant numerical dispersion (with extensive trailing waves) when computed by the standard sixthorder scheme. The 7point stencil DRP scheme is better. Both schemes, having the same stencil size, require the same amount of computation. The solution by the 15point stencil is even better. For all intents and purposes, the solution is identical numerically to the exact solution.
Leave a reply