Numerical geometry and constraint monitoring

Use the first circle factor of the fixed periodic three-torus from the earlier fixed-gauge data definition. This item is a separate nondimensionalized scalar model: both the spatial coordinate $x$ and time $t$ are dimensionless, the period is $2\pi$, and the model wave speed is $1$. Let $N\ge 3$ be an integer, $h:=2\pi \mathrm{/}N$, and $x_j:=jh(\mathrm{m}\mathrm{o}\mathrm{d}2\pi )$ for $j\in \mathbb{Z}\mathrm{/}N\mathbb{Z}$. A periodic grid function is an $N$-tuple $u=(u_j{)}_{j\in \mathbb{Z}\mathrm{/}N\mathbb{Z}}$ of real numbers. Define the centered periodic second-difference operator

where subscripts are reduced modulo $N$.

The algorithm takes as input such arrays $u^0,v^0$, a real time step $k>0$, and an integer $n_{\max }\ge 1$. It sets

then, for each integer $n=1,\dots ,n_{\max }-1$, computes

It outputs $(u^0,u^1,\dots ,u^{n_{\max }})$. All operations here are exact real arithmetic. No rounding-error model or claim about a floating-point implementation is part of the definition.

Fixed-gauge finite-slab data and solution topology

The recurrence is the centered finite-difference model for the scalar equation $u_{tt}=u_{xx}$. It is a one-dimensional diagnostic benchmark only; its stability properties do not establish convergence of a discretization of the Einstein equations. The initialization uses the second-order time expansion with the discrete spatial operator. The time step and mesh spacing are nondimensional, so their ratio has no unit conversion hidden in it.

For every admissible integer $N\ge 3$, time-step count $n_{\max }\ge 1$, real step $k>0$, and pair of real periodic input arrays, the algorithm terminates after exactly $n_{\max }-1$ update passes and outputs all levels through $u^{n_{\max }}$. In the unit-cost exact-real-arithmetic model, where each scalar arithmetic operation and each array read or write has bounded cost independent of $N$ and $n_{\max }$, its operation count is $O(N{n}_{\max })$.

A finite periodic scalar-wave benchmark

Let $N\ge 3$, $h=2\pi \mathrm{/}N$, $k>0$, and $\lambda =k\mathrm{/}h$. Complexify the real linear recurrence of the periodic wave algorithm. For each discrete Fourier index $\mathrm{\ell }\in \{0,1,\dots ,N-1\}$ define the mode

If $0<\lambda <1$, every nonconstant mode ($\mathrm{\ell }\mathrm{\ne }0$) has recurrence coefficient

and two distinct amplification roots

Consequently, each fixed nonconstant mode has bounded amplitude on every finite slab. For arbitrary real initial arrays define their spatial means by

The constant mode has the repeated root $1$ and evolves exactly as ${\bar{u}}^0+nk{\bar{v}}^0$; for any fixed $T\ge 0$ it is bounded on the slab $0\le nk\le T$ but has no uniform-in-time displacement bound when ${\bar{v}}^0\mathrm{\ne }0$. This statement does not assert that the condition $0<\lambda <1$ alone proves convergence of a numerical relativity scheme. Define the nonzero discrete frequencies ${\omega }_{\mathrm{\ell }}:=2h^{-1}\mathrm{\mid }\sin (\pi \mathrm{\ell }\mathrm{/}N)\mathrm{\mid }$, $1\le \mathrm{\ell }<N$. On this fixed torus they satisfy the mesh-independent lower bound ${\omega }_{\mathrm{\ell }}\ge \frac{2}{\pi }(N\ge 3,1\le \mathrm{\ell }<N)\mathrm{.}$ In particular each fixed grid has a positive minimum. This frequency bound is not a mesh-independent amplitude bound for arbitrary unscaled two-step data $(a_0,a_1)$: the distinct-root coefficients in the mode solution can depend on the grid and the initial data. The constant-mode qualifications above are unchanged. For even $N$, the alternating mode with Fourier index $\mathrm{\ell }=N\mathrm{/}2$ is called the Nyquist mode. At the endpoint $\lambda =1$, if $N$ is even and $u_j^0=a_0(-1{)}^j$, $v_j^0={\beta }_0(-1{)}^j$ for real amplitudes $a_0,{\beta }_0$, then each output level satisfies

The Nyquist amplification root is repeated at $-1$, and a nonzero Nyquist velocity amplitude ${\beta }_0$ gives linear growth in the step number.

A finite periodic scalar-wave benchmark

The complex exponential by its power series

The complex numbers as $\mathbb{R}[x]\mathrm{/}(x^2+1)$, with the real embedding and imaginary unit $i$

$\mathbb{C}=\mathbb{R}[x]\mathrm{/}(x^2+1)$ is a field, every element is uniquely $a+bi$, and every nonzero element has inverse $(a-bi)\mathrm{/}(a^2+b^2)$

Real and imaginary parts, complex conjugation, and modulus

Conjugation is an involutive real-field automorphism, $z\bar{z}=\mathrm{\mid }z{\mathrm{\mid }}^2$, and modulus is definite, multiplicative, and subadditive

Absolute value in an ordered field

Basic properties of the absolute value

The triangle inequality

Squares of nonzero elements are positive

Every complete ordered field is Archimedean

Euler's formula: $\exp (i\theta )=\cos \theta +i\sin \theta$ for every real $\theta$

The addition formulas for sine and cosine

The derivatives of sine and cosine are cosine and minus sine

Quarter-turn values and shifts by pi/2 and pi

Parity and the Pythagorean identity for sine and cosine

The zero sets of sine and cosine and the least positive common period 2 pi

Existence and uniqueness of $n$-th roots: a unique $a^{1\mathrm{/}n}\ge 0$ with $(a^{1\mathrm{/}n}{)}^n=a$

Finite sums and finite products, by recursion

Finite commutative-monoid sums are invariant under bijective reindexing, split over disjoint unions, and satisfy the finite Fubini rule

The mean value theorem, as the case $g(x)=x$ of Cauchy's: for $f$ continuous on $[a,b]$ with $a<b$ and differentiable on $(a,b)$ there is $c\in (a,b)$ with $f(b)-f(a)=f^{\mathrm{\prime }}(c)(b-a)$

Pi is the first positive zero of sine

The root calculation is mode-by-mode. A convergence theorem for a discretized Einstein system would also need an error estimate for the initialization and a discrete energy estimate for the actual reduced equations, including its boundary treatment. The analytic finite-slab estimate on Finite-slab stability in a fixed hyperbolic gauge concerns continuum solutions of a fixed gauge system and does not supply those discrete estimates.

Let $I\subseteq \mathbb{R}$ be open, and let $u:\mathbb{R}\times I\to \mathbb{R}$ be six times continuously differentiable and $2\pi$-periodic in its spatial variable $x$. Assume the nondimensional scalar wave equation $u_{tt}=u_{xx}$ on $\mathbb{R}\times I$. For $N\ge 3$, $h=2\pi \mathrm{/}N$, $k>0$, $x_j=jh$, and an integer $n$ with $t_n=nk$ and $[t_n-k,t_n+k]\subset I$, put

Set $M_n$ to be the supremum

Then for every grid index $j$,

Thus the normalized local truncation residual is $O(k^2+h^2)$ on each compact time window where the displayed fourth-derivative bound is uniform. This consistency estimate alone is not a convergence theorem.

A finite periodic scalar-wave benchmark

A uniform derivative bound gives a uniform Taylor remainder bound

Clairaut--Schwarz theorem for continuous second partial derivatives

Heine-Borel in ${\mathbb{R}}^n$: with the Euclidean metric a subset of ${\mathbb{R}}^n$ is compact if and only if it is closed and bounded, and the proof by bisection uses no choice principle; the same holds on the real line

A continuous real-valued function on a nonempty compact metric space is bounded and attains a greatest and a least value

Absolute value in an ordered field

The triangle inequality

Use the dimensionless angular coordinates $x^i$ and the periodic spatial three-torus from the earlier fixed-gauge data definition. Choose integers $N_i\ge 3$, coordinate spacings ${\mathrm{\Delta }}_i:=2\pi \mathrm{/}{N}_i$, and the finite periodic grid

At each grid node $p$, supply a real symmetric positive-definite matrix $h_{ij}(p)$ and its matrix inverse $h^{ij}(p)$ (so $h^{ik}(p)h_{kj}(p)={\delta }^i{}_j$), a real symmetric array $K_{ij}(p)$, a scalar $\rho (p)$, a covector array $j_i(p)$, and real constants $\mathrm{\Lambda }$ and $\kappa >0$. These are the sampled spatial metric, second fundamental form, energy density and momentum density. For a future unit normal $n$ and the coordinate basis ${\mathrm{\partial }}_i$, use the convention $K(X,Y)=g({\mathrm{\nabla }}_{X}n,Y)$, $\rho =T(n,n)$ and $j_i=-T(n,{\mathrm{\partial }}_i)$, where $T$ is a smooth symmetric covariant two-tensor. The array definition itself needs only the listed samples and does not assume the continuum field equation. The vacuum case is $\rho =j_i=\mathrm{\Lambda }=0$. In the SI angular-coordinate convention used here, $[x^i]=1$, $[h_{ij}]={\mathrm{m}}^2$, $[K_{ij}]=\mathrm{m}$, $[\rho ]=\mathrm{J}{\mathrm{m}}^{-3}$, $[j_i]=\mathrm{J}{\mathrm{m}}^{-2}$, $[\mathrm{\Lambda }]={\mathrm{m}}^{-2}$, and $[\kappa ]={\mathrm{m}}^{-1}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{s}}^2$. Thus the scalar residual below has units ${\mathrm{m}}^{-2}$ and its covector components in this angular basis have units ${\mathrm{m}}^{-1}$.

Let ${\mathbb{R}}^{G_{\mathbf{\Delta }}}$ denote the real arrays indexed by the grid. Choose and specify linear maps ${\delta }_i:{\mathbb{R}}^{G_{\mathbf{\Delta }}}\to {\mathbb{R}}^{G_{\mathbf{\Delta }}}$ with ${\delta }_{i}1=0$; apply them componentwise to tensor arrays. Periodic index wrapping is used here. Define the discrete connection components, curvature, Ricci array and scalar curvature by nodewise operations, with every index $i,j,k,\mathrm{\ell },a,b$ ranging over $1,2,3$ and repeated upper/lower indices summed:

Set

The signed discrete Hamiltonian and momentum residuals are

These signs match the Hamiltonian and momentum constraint equations when the pointwise equation $G+\mathrm{\Lambda }g=\kappa T$ is assumed. The displayed discrete formulas are array operations: they do not define a coordinate-covariant discrete connection or assume a discrete Leibniz rule. A reconstruction map has the exact type

where ${\mathcal{D}}_{\mathbf{\Delta }}$ is the set of all grid arrays $(h_{ij},K_{ij},\rho ,j_i)$ with the pointwise symmetry and positive-definite conditions above, and $\mathcal{C}$ is the set of pointwise assignments on ${\mathbb{T}}_{2\pi }^3$ of a positive-definite symmetric covariant two-tensor $h$, a symmetric covariant two-tensor $K$, a scalar $\rho$ and a covector $j$, with the same units. No continuity or differentiability of a member of $\mathcal{C}$ is implied; if continuum derivatives are used, the needed regularity must be stated separately. Specify ${\mathcal{I}}_{\mathbf{\Delta }}$ as part of each diagnostic, since it records which continuum field the grid is intended to represent. For every reported diagnostic, also record the grid norm and its weights, the mesh hierarchy, the units, the derivative maps, the periodic boundary treatment, and ${\mathcal{I}}_{\mathbf{\Delta }}$. A small residual by itself does not assert closeness to an exact constraint-satisfying development.

Fixed-gauge finite-slab data and solution topology

Gauss-Codazzi decomposition in Gaussian normal coordinates and the constraint equations in terms of the induced metric and second fundamental form

The geometric convention in the constraint theorem fixes the signs above. Gourgoulhon's 3+1 lecture notes use the opposite sign for the second fundamental form in their convention; their displayed equations are therefore context and cross-checks only, not a source for silently replacing the local sign convention.