Finite-Difference Operators

A finite-difference stencil is more than a Taylor-series formula. In a relativistic evolution it must combine with boundaries, dissipation, time integration, and constraints to produce a stable discrete system. Interior order alone is therefore an incomplete specification. This page builds the operator, derives its discrete energy structure, and attaches convergence and high-frequency failure tests.

Required background. Ordinary Differential Equations supplies semidiscrete systems, stability functions, and convergence concepts.

Helpful background. Formulations of General Relativity provides the continuum evolution and constraint systems being discretized, while Partial Differential Equations and Green Functions supplies hyperbolicity, boundary data, and energy estimates.

Consistency in the interior

On a uniform grid $x_j=jh$, the centered first derivative

$$(D_0u)_j=\frac{u_{j+1}-u_{j-1}}{2h} =u'(x_j)+\frac{h^2}{6}u'''(x_j)+O(h^4)$$

is second-order consistent for a smooth function. Taylor expansion determines local truncation error, not stability. The stencil also has a blind Nyquist mode: for $u_j=(-1)^j$, $D_0u=0$. Unresolved grid-scale oscillations may therefore survive unless the system, boundary closure, or explicit high-order dissipation controls them.

On curvilinear or nonuniform grids, metric factors and Jacobians belong in the operator definition. Applying a uniform-grid formula after a nonlinear map can reduce order or violate a conservation identity. The distinction among consistency, stability, and convergence is developed in Gustafsson, Kreiss, and Oliger 2013, Chapters 2–4. The implementation record must state grid locations, ghost zones, stencil coefficients, mapping, derivative order, and whether variables are nodal, cell centered, or staggered.

Summation by parts and boundary closure

A summation-by-parts derivative has a positive norm matrix $H$ and satisfies

$$HD+D^TH=B,$$

where $B$ represents the two boundary terms. This is the discrete analogue of integration by parts. For scalar advection $u_t+a u_x=0$, multiplying the semidiscrete equation by $u^TH$ exposes the boundary energy flux. An incoming datum can be imposed weakly through a simultaneous-approximation-term penalty whose coefficient is chosen so the discrete energy does not grow spuriously.

For $a\gt0$ on nodes $0,\ldots,N$, take $B=e_Ne_N^T-e_0e_0^T$ and impose left data $g(t)$ through

$$\dot{\mathbf u}=-aD\mathbf u+\tau H^{-1}e_0(u_0-g).$$

Here $e_0$ and $e_N$ select the endpoint values. With $E_h=\mathbf u^TH\mathbf u/2$ and homogeneous data,

$$\dot E_h=-\frac a2u_N^2+\left(\frac a2+\tau\right)u_0^2.$$

Thus $\tau\le-a/2$ gives a nonincreasing discrete energy. The common choice $\tau=-a$ also gives, for nonzero data,

$$\dot E_h=-\frac a2u_N^2-\frac a2(u_0-g)^2+\frac a2g^2.$$

The outflow removes energy, mismatch at the inflow is damped, and the prescribed data supply the allowed energy. This derivation fixes the penalty sign before any numerical pulse is evolved. The choice $\tau=-1$ below corresponds to $a=1$.

Boundary rows usually have lower local order than the centered interior, yet the global convergence order can be higher than that row order. Judge the complete initial-boundary scheme rather than quoting the interior stencil. The energy and convergence framework is developed in Gustafsson, Kreiss, and Oliger 2013, Chapters 4–5 and 8–11. Strand 1994, §§ 2–4 and Tables 1–4 constructs diagonal-norm SBP families through eighth-order interior accuracy; their derivative coefficients, norm, and boundary closure form one operator and must not be mixed across tables.

Periodic and finite-interval benchmarks

First solve periodic advection on $0\le x\lt 2\pi$,

$$u_t+u_x=0, \qquad u(0,x)=\sin x, \qquad u(t,x)=\sin(x-t).$$

Use the displayed $D_0$ operator with periodic indexing, classical RK4, $N=64,128,256$ equally spaced points, Courant number $\Delta t/h=1$, and final time $t=2\pi$. Shorten the last step to land exactly at $t=2\pi$. Report $L^2$ and $L^\infty$ errors, the Fourier phase error of the $k=1$ mode, and the observed order

$$p_{\mathrm{obs}}=\log_2\frac{E_h}{E_{h/2}}.$$

Then repeat on $0\leq x\leq1$ with the fourth-order-interior, second-order-boundary diagonal-norm SBP operator from Strand and a left simultaneous-approximation-term penalty $\tau H^{-1}e_0(u_0-g)$ with $\tau=-1$. Take $g=0$ and

$$u(0,x)=B(x),\qquad B(x)= \begin{cases} \exp\!\left[1-\dfrac{1}{1-s^2}\right],&\lvert s\rvert\lt 1,\\ 0,&\lvert s\rvert\geq1, \end{cases} \qquad s=\frac{x-0.25}{0.10}.$$

Use $N=64,128,256$ equal subintervals, RK4, and $\Delta t/h=0.2$. Compare with the zero-extended exact translate $B(x-t)$ at $t=0.5$, then continue to $t=1.2$, after the exact pulse has left, and quote the reflected $L^2$ and $L^\infty$ amplitudes. Monitor $u^THu$, the inflow residual $u_0-g$, and the outflow flux. The periodic test isolates the interior; this finite-interval test makes the closure and penalty sign falsifiable. A second run with $\Delta t/h=0.1$ separates temporal from spatial error.

High-frequency control and constraints

Kreiss–Oliger-type dissipation can damp unresolved modes while leaving smooth convergence order intact when scaled consistently. It is not a cure for an unstable boundary penalty or a violated CFL limit. Scan the dissipation strength, report its effect on phase and amplitude, and include a grid-scale initial perturbation in addition to the smooth sine wave.

For constrained relativistic systems, monitor both field error and constraint norm. A stencil can propagate the physical field stably while exciting a constraint mode. Calabrese et al. 2006, §§ 5–5.1, Eq. 119, and Figs. 2–5 illustrates numerical-stability and noisy-convergence tests for formulations of Einstein’s equations. The authors explicitly distinguish their test from the Apples-with-Apples robust-stability test, so the two names are not interchangeable. Translate the source’s continuum variables, normalization, and signs before adapting either test.

Reproducible environment and failure criteria

Pin the implementation revision, language and compiler, precision, grid generator, stencil and SBP coefficients, norm matrix, penalty parameters, dissipation, time integrator, Courant factor, boundaries, and output norms. A clean rerun emits inputs, errors on at least three grids, energy and constraint histories, high-mode amplitudes, and the exact comparison. This specification does not imply that the calculation has been executed.

The validity regime assumes the regularity needed by the designed order and a well-posed continuum problem. Reject unstable boundary closures, a wrong penalty sign, order reduction hidden by field norms, a CFL violation masked by excessive dissipation, or apparent convergence at only two resolutions. Also reject a smooth-wave success if a Nyquist perturbation grows without resolution convergence. Spatial and temporal errors must be separated by repeating with a smaller timestep at fixed grid.

Check your understanding

A fourth-order interior stencil with a second-order boundary closure produces second-order convergence in the measured norm. Is the implementation necessarily wrong?

Solution

No. The global order depends on the boundary closure, norm, stability estimate, and how boundary error propagates. Compare with the theoretical order of the complete SBP operator, refine at fixed Courant number, and reduce the timestep independently. It is wrong only if the observed order disagrees with the documented full scheme after pre-asymptotic and temporal errors are excluded.

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.
• Gustafsson, B., Kreiss, H.-O., and Oliger, J. (2013), Time-Dependent Problems and Difference Methods, 2nd edition, Wiley. DOI.
• Strand, B. (1994), “Summation by Parts for Finite Difference Approximations for d/dx,” Journal of Computational Physics 110, 47–67. DOI.