Finite-difference stencils and numerical dispersion

On a periodic uniform grid, insert the mode exp(ij theta) into a linear translation-invariant stencil. A centered first difference multiplies it by i sin(theta)/h, whereas the exact derivative on a resolved continuum mode multiplies it by ik with theta=kh. Expanding sine shows second-order consistency when k is fixed and h tends to zero. That low-frequency expansion is not uniform over modes near the grid’s highest representable frequencies, where theta remains order one.

The three-point second difference has symbol minus four sin squared(theta/2) over h squared. A fourth-order centered first derivative instead has symbol i times[8sin(theta)−sin(2theta)] divided by 6h. These formulas let one compare stencils without confusing their coefficient patterns with an error estimate. A nominally high-order stencil can still have substantial error near the resolution limit, and wide boundary closures need their own consistency and stability checks.

For u_t+a u_x=0 with a>0, centered spatial differencing gives a semidiscrete Fourier mode with angular frequency a sin(theta)/h. Dividing by physical wave number k gives phase speed a sin(theta)/theta, with its continuous limit a at theta zero. Differentiating angular frequency with respect to k gives group velocity a cos(theta), describing the movement of a narrow wave packet. A high-frequency packet can have a negative group velocity even though the PDE transports every profile to the right.

This behavior is numerical dispersion: different wavelengths acquire different phase shifts. It is distinct from dissipation, which changes amplitudes. The periodic centered operator is skew-adjoint in the h-weighted L2 inner product, so its exact semidiscrete flow preserves that norm. Norm preservation does not imply that a wave packet is in the correct location. A computation can have almost no amplitude loss and still produce an inaccurate signal through accumulated phase error.

Forward Euler applied to the centered semidiscrete advection equation has multiplier one minus i times the Courant number times sin(theta). Its squared magnitude exceeds one for every nonzero Courant number and every mode with nonzero sine. Thus decreasing the step to satisfy the geometric advection CFL condition alone does not make this scheme stable on a fixed-ratio refinement. The spatial operator’s nondissipation and the temporal method’s instability are separate facts that happen to coexist.

Other time integrators can have stable intervals on the imaginary axis, and upwind spatial differencing adds dissipation rather than preserving the centered symbol. The right choice depends on the PDE, the target accuracy and the resolved spectrum. To interpret a dispersion study, report the spatial symbol, the time integration rule and the boundary setting. Otherwise an observed phase lag may be attributed to the wrong component of the algorithm, or an unstable amplitude trend may be mistaken for a physical wave phenomenon.