Stability, Constraints, and CFL Limits

A numerical evolution can remain bounded while drifting away from the constraint surface of the continuum formulation. Conversely, a convergent spatial operator can become unstable when the timestep is too large. Stability, constraint preservation, and accuracy are therefore separate questions. A useful test measures field energy, high-frequency content, and constraint norms while scanning the Courant factor and resolution.

Required background. Ordinary Differential Equations supplies amplification factors, eigenvalues, and time-integration stability.

Helpful background. Formulations of General Relativity provides continuum constraint propagation and characteristic speeds, while Partial Differential Equations and Green Functions supplies energy estimates and domains of dependence.

The discrete domain of dependence

For a hyperbolic equation with characteristic speed $a$, an explicit local scheme usually requires

$$\Delta t\le C_{\max}\frac{\Delta x}{\lvert a\rvert},$$

with a method-dependent $C_{\max}$. In several dimensions, on mapped meshes, or with refinement, the relevant bound uses the largest local coordinate characteristic speed and the smallest effective spacing. Courant, Friedrichs, and Lewy 1928, § 2 derives the necessary mesh restriction by comparing the difference equation’s domain of dependence with the hyperbolic equation’s characteristic domain. That condition is necessary for convergence; it is not sufficient to make every discretization stable.

For a linear constant-coefficient periodic problem, von Neumann analysis diagonalizes translation-invariant stencils into Fourier modes and tests their amplification factors. The relation among stability, consistency, and convergence for linear initial-value problems is developed in Lax and Richtmyer 1956, Part I §§ 5–8, pp. 271–275. Boundaries, variable coefficients, and nonnormal systems require an energy estimate or stronger analysis beyond the periodic symbol.

Wave-and-constraint benchmark

Write the one-dimensional wave equation as

$$\partial_tu=\Pi, \qquad \partial_t\Pi=\partial_x\Phi, \qquad \partial_t\Phi=\partial_x\Pi,$$

with auxiliary constraint

$$\mathcal C=\Phi-\partial_xu=0.$$

In the continuum, consistent initial data preserve $\mathcal C=0$. A discrete scheme may approximate the three evolution equations stably while boundary closures, interpolation, or inconsistent derivative operators inject constraint error. Monitor

$$E_h=\frac12\sum_jh\left(\Pi_j^2+\Phi_j^2\right), \qquad \lVert\mathcal C\rVert_{2,h},$$

along with field error and a high-frequency norm.

For a completely specified baseline, take $0\leq x\lt 2\pi$, $N=64,128,256$ uniform points, periodic centered differences

$$(D_0f)_j=\frac{f_{j+1}-f_{j-1}}{2h},$$

and classical four-stage RK4. Initialize

$$u_j=\sin x_j, \qquad \Pi_j=-(D_0u)_j+10^{-12}q_j, \qquad \Phi_j=(D_0u)_j+10^{-12}r_j,$$

where $q_j$ and $r_j$ are independent uniform samples on $[-1,1]$ produced by NumPy’s PCG64 generator with seed $20260818$. The unperturbed part is a right-moving semidiscrete wave and satisfies the discrete constraint exactly; the perturbation gives a known $10^{-12}r_j$ constraint seed.

Use exactly the same $D_0$ in the two derivative equations and in $\mathcal C_h=\Phi-D_0u$. Then

$$\dot{\mathcal C}_h=D_0\Pi-D_0\Pi=0.$$

A common Runge–Kutta update preserves this linear invariant up to roundoff at every step, including steps above the stability threshold. The energy can therefore blow up while the constraint remains near its seed until roundoff is amplified. This baseline tests derivative consistency and shows why a small constraint is insufficient evidence of stability. A separate negative test can replace one occurrence of $D_0$ by a different stencil; the resulting nonzero $\dot{\mathcal C}_h$ must then be explained by that operator difference.

For Fourier angle $\theta$, $D_0$ has symbol $i\sin\theta/h$, while RK4 is stable on the imaginary axis for $\lvert z\rvert\leq2\sqrt2$. Thus the predicted Courant limit for the unit-speed periodic problem is $C_{\max}=2\sqrt2$. Evolve to $t=40\pi$ at $C=\Delta t/h=2.0,2.7,2.82,2.84,3.0$, shortening only the last step. Record the maximum energy, the discrete constraint relative to its initial seed, and the Fourier norm of modes satisfying $\lvert\sin\theta\rvert\ge0.9$, together with the individual modes $\theta=\pm\pi/2$. These have the largest derivative eigenvalues; the Nyquist mode $\theta=\pi$ is actually annihilated by $D_0$, so measuring only the largest wavenumbers can miss the instability. Below threshold, the random component should remain bounded and the smooth error should converge; above it, an excited mode must grow more strongly as $N$ increases. Repeat $C=2.7$ at $N=512$ to check the apparently stable side of the boundary rather than inferring it from three coarse runs.

Constraint damping and projection

Constraint damping modifies the evolution so small violations decay, while projection maps numerical data back onto the constraint surface. Both can be useful, but neither licenses hiding the undamped behavior. Report the damping coefficient in physical units or relative to the grid, derive how it changes characteristic or stiff timescales, and rerun the CFL scan. A coefficient proportional to $1/h$ can impose a new timestep restriction.

Projection may reduce $\lVert\mathcal C\rVert$ discontinuously while altering conserved quantities or physical modes. Measure that change. A stable field norm with a growing constraint is a failed formulation–discretization pairing for long evolutions, even if the waveform looks smooth.

The role of boundary stability and the extension from periodic symbols to initial-boundary problems is treated in Lax and Richtmyer 1956, Part II §§ 12–14, pp. 279–281. Calabrese, Hinder, and Husa 2006, §§ 5–5.1, Eq. 119, and Figs. 2–5 uses random-data numerical-stability and noisy-convergence tests to expose high-frequency failures that smooth data can miss. That paper explicitly distinguishes this procedure from the Apples-with-Apples robust-stability test; the labels should not be conflated.

Choosing an operational timestep

The analytic bound should be checked empirically, not estimated from the first run that survives. For each grid, record maximum field energy, constraint norm, and highest-mode amplitude over a fixed number of crossing times. Define failure thresholds before scanning. The production Courant factor should include a safety margin for variable characteristic speeds, moving meshes, nonlinear states, and different boundary patches.

Convergence above a threshold cannot be restored by adding enough artificial dissipation if the underlying amplification factor remains outside the time integrator’s stable region. Dissipation can move eigenvalues, but then the combined operator must be analyzed and its accuracy measured. A short run that ends before the unstable mode grows is not evidence of stability.

Reproducible environment and failure criteria

Pin the implementation revision, precision, derivative and dissipation operators, boundary closure, integrator tableau, timestep rule, characteristic-speed estimate, damping or projection parameters, random seed, mesh sequence, and norm definitions. A clean rerun exports the Courant scan, energy and constraint histories, Fourier-mode amplitudes, refinement errors, and the predicted stability boundary.

The validity regime is the declared mesh, boundary treatment, coefficient range, and time interval. Failure modes include a CFL violation, resolution-amplified noise, constraint growth hidden by a field norm, an incorrect characteristic speed, and damping that introduces unrecognized stiffness. Reject stability inferred from one resolution, one smooth input, or a final-time snapshot without growth history.

Check your understanding

A run has bounded field energy and convergent $u$, but $\lVert\mathcal C\rVert_{2,h}$ grows exponentially at a rate that increases with resolution. Is the discretization stable for the constrained system?

Solution

No. The physical-looking variables do not exhaust the state’s required properties. Resolution-amplified constraint growth indicates a discrete constraint or boundary instability. Inspect the constraint-propagation operator, boundary closure, derivative consistency, and Courant factor; then repeat the random-noise and convergence tests with the constraint norm as an acceptance condition.

References

• Calabrese, G., Hinder, I., and Husa, S. (2006), “Numerical Stability for Finite Difference Approximations of Einstein’s Equations,” Journal of Computational Physics 218, 607–634. arXiv.
• Courant, R., Friedrichs, K., and Lewy, H. (1928), “On the Partial Difference Equations of Mathematical Physics,” Mathematische Annalen 100, 32–74. DOI.
• Lax, P. D., and Richtmyer, R. D. (1956), “Survey of the Stability of Linear Finite Difference Equations,” Communications on Pure and Applied Mathematics 9, 267–293. DOI.