Conventions and normalizations

These definitions are normative for the whole library.

Notation

  • \(\alpha\) is a fractional time order, \(0<\alpha<1\).

  • \(s\) is a fractional power of the Laplacian, \(0<s<1\).

  • \(d\) is the spatial dimension.

  • \((-\Delta)^s\) is always written together with its realization, because the Dirichlet spectral, zero-exterior Riesz/restricted, and periodic Fourier operators are different operators.

Orders are immutable and must lie in the open unit interval.

Boundary terminology

Fourier transform

For integrable \(u:\mathbb R^d\to\mathbb C\),

\[ \widehat u(\xi)=\int_{\mathbb R^d}e^{-i x\cdot\xi}u(x)\,dx, \qquad u(x)=\frac{1}{(2\pi)^d} \int_{\mathbb R^d}e^{i x\cdot\xi}\widehat u(\xi)\,d\xi , \]

so the whole-space fractional Laplacian has multiplier

\[ \widehat{(-\Delta)^s u}(\xi)=|\xi|^{2s}\widehat u(\xi). \]

On a periodic box of lengths \(L_j\), the continuous frequencies become \(2\pi k_j/L_j\). PeriodicFractionalLaplacian applies the corresponding discrete Fourier multiplier to the unique periodic nodal values and assigns zero to the constant mode.

Riesz normalization

Consistency with that multiplier fixes

\[ C_{d,s} =\frac{4^s\Gamma(d/2+s)} {\pi^{d/2}|\Gamma(-s)|} =\frac{2^{2s}s\Gamma(d/2+s)} {\pi^{d/2}\Gamma(1-s)}, \qquad C_{2,s}=\frac{2^{2s}s\Gamma(1+s)}{\pi\Gamma(1-s)} . \]

All Riesz backends share this constant.

Finite-element action

The Riesz implementation assembles the operator through its integral against every finite-element test function \(v\),

\[ \langle v,Lu\rangle=\int_\Omega v(x)Lu(x)\,dx , \]

instead of evaluating the pointwise formula for \(Lu\) at mesh nodes and then interpolating those values. The dense, matrix-free, and H-matrix backends all compute this same Galerkin action. Since the continuous operator is symmetric, the discrete identity \(\langle v,Lu\rangle=\langle u,Lv\rangle\) is used to check the implementation. The assembled action is left unchanged so that any discrete asymmetry remains visible instead of being hidden by post-processing.

Numerical precision and requested accuracy

Quadrature rates and weights used in simulations are computed with standard 64-bit floating-point arithmetic (float64). Independently evaluated higher-precision values are used for validation, not during simulations.

For spectral sinc quadrature, a requested target below meaningful float64 precision is raised to machine epsilon and emits a warning. A request that would need more than 100,000 sinc nodes is rejected with an explanation rather than silently using a different rule. Diethelm2022 instead reports whether its fixed mode count meets target_error through the target_achievable metadata field.