Group Velocity Consideration

There is no question that a good finite difference scheme should have a Ax equal to a Ax over as wide a bandwidth as possible. However, from the standpoint of wave propagation, it is also important to ensure that the group velocity of the finite difference scheme is a good approximation of that of the original equation. It will be shown later that the group velocity of a finite difference scheme is determined by da /da. da /da, the slope of the a Ax — a Ax curve, should be equal to 1 if the scheme is to reproduce the same group velocity or speed of propagation of the original partial differential equation.

Figure 2.2 shows the da/da curve of the optimized 7-point stencil scheme as a function of a Ax. It is seen that for 0.4 < aAx < 1.15, da/da is larger than 1. This will lead to numerical dispersion, because some wave components would propagate faster than others and faster than the wave speed of the original partial differential equation. Around aAx = 0.8 to 1.0, the slope is greater than 1.02. Thus, the wave component in this wave number range will propagate at a speed 2 percent faster (assuming time is computed exactly). That is to say, after propagating over 100 mesh spacings, the group of fast waves has reached the 102 mesh point. This is too much dispersion for long-range propagation.

One way to reduce numerical dispersion is to reduce the range of optimization. Instead of the range of aAx = 0 to n/2, one may use the range of 0 to n. That is, the integrated error is

n

E = j |aAx — aAx|2 d (aAx).

0

By setting n to a specified set of values, the coefficients ax, a2, and a3 of the 7-point central difference stencil may be found as before. Upon examining the corresponding da/da curve, it is recommended that the case n = 1.1 has the best overall wave

Figure 2.3. aAx versus a Ax.————- ,

optimized fourth-order scheme;——- ,

standard sixth-order scheme;.. , stan­dard fourth-order scheme; , stan­

dard second-order scheme.

propagation characteristics. The stencil coefficients are

a0 = 0 a1 = – a-1 = 0.77088238051822552

a2 = – a-2 = -0.166705904414580469 a3 = – a-3 = 0.02084314277031176.

Figure 2.3 shows the a Ax versus a Ax relations for the standard sixth-order, fourth-order, and second-order central difference scheme and the optimized scheme (n = 1.1). The values of aAx at approximately the point of deviation from the a Ax – aAx straight-line relation for these schemes, are

Scheme	acAx	Resolution, A.
Sixth-order	0.85	7.4
Fourth-order	0.60	10.5
Second-order	0.30	21.0
Optimized	0.95	6.6

where X is the wavelength. Thus, a high-order finite difference scheme can better resolve small-scale features.

Figure 2.4 is a comparison of the curves. Figure 2.5 is the same plot but with an enlarged vertical scale. The increase in group velocity near aAx = 0.7 is around 0.3 percent. This means that the scheme has very small numerical dispersion. The resolved bandwidth for | dal da -1.01 < 0.003 is about 15 percent wider than that of the standard sixth-order scheme. This improvement is achieved at no cost because both schemes have a 7-point stencil.

example. Consider the initial value problem governed by the one-dimensional convective wave equation. Dimensionless variables with Ax as the length scale,

Ax/c (c = sound speed) as the timescale will be used. The mathematical problem is

d u	A = 0, d X	— ^ < X < ^.	(2.15)
	t = 0	u = f( x ).	(2.16)

The exact solution is

u = f (x —1).

Figure 2.5. Enlarged d(fAX) versus aAx

curves.—— , optimized scheme and the

…… , standard sixth-order central differ­ence scheme.

In other words, the solution consists of the initial disturbance propagating to the right at a dimensionless speed of 1 without any distortion of the waveform. In particular, if the initial disturbance is in the form of a Gaussian function,

with a half-width of three mesh spacings, then the exact solution is

(2.18)

On using a (2N + 1)-point stencil finite difference approximation to d/dx in Eq. (2.15), the original initial value problem is converted to the following system of ordinary differential equations:

t = 0 ul = f (l ).

l is the spatial index.

For the Gaussian pulse, the initial condition is

The system of ordinary differential equations (2.19) and initial condition (2.21) can be integrated in time numerically by using the Runge-Kutta or a multistep method. The solutions using the standard central second-order (N = 1), fourth-order (N = 2), sixth-order (N = 3), and the optimized 7-point (N = 3, n = 1.1) finite difference schemes at t = 400 are shown in Figure 2.6. The exact solution in the form of a Gaussian pulse is shown as the dotted curves. The second-order solution is in the form of an oscillatory wave train. There is no resemblance to the exact solution. The numerical solution is totally dispersed. The fourth-order solution is better. The sixth – order scheme is even better, but there are still some dispersive waves trailing behind the main pulse. The wavelength of the trailing wave is about 7, which corresponds to a wave number of 0.9, which is beyond the resolved range of the sixth-order scheme. The optimized scheme gives very acceptable numerical results.

To understand why the different finite difference scheme performs in this way, it is to be noted that the Fourier transform of the initial disturbance is

It is easy to verify that a significant fraction of the initial wave spectrum lies in the short-wave range (a Ax > ac Ax) of the a Ax versus a Ax curve of the standard second – and fourth-order finite difference schemes. (Note: Ax = 1 in this problem.) These wave components propagate at a slower speed (smaller than the group velocity). This causes the pulse solution to disperse in space and exhibits a wavy tail.

x/Ax Figure 2.6. Comparison between the computed and the exact solution of the simple one­dimensional convective wave equation. (a) second-order, (b) fourth-order, and (c) sixth-order standard central difference scheme, (d) 7-point optimized scheme.—- , numerical solution; ………… , exact solution.

Another way to understand the dispersion effect is to superimpose the f (a) curve on the dal da versus a curve. The standard sixth-order scheme maintains a wave speed nearly equal to 1 (da/da = 1.0) up to a Ax = 0.85. Beyond this value the wave speed decreases. The optimized scheme maintains an accurate wave speed up to a Ax = 1.0. Now, for the given initial condition, there is a small quantity of wave energy around a Ax = 0.9. This component, having a slower wave speed, appears as trailing waves in the solution of the standard sixth-order scheme. The wavelength of the trailing wave is X = (2n|a) ~ 7.0. This is consistent with the numerical result shown in Figure 2.6. For the optimized scheme, there are trailing waves with wave number ~ 1.0 or X ~ 6.3, but the amplitude is much smaller and is, therefore, not so easily detectable. Basically, none of these schemes can resolve waves with wavelengths less than six mesh spacings really well.