Fast Spectral Simulation of One-Dimensional Nonlinear Fractional Waves in Attenuating Media

1. Introduction

Wave propagation in media with memory, such as porous and fluid-saturated rocks, biological tissue and polymers, is commonly described by integro-differential equations with fractional time derivatives [1,2,3,4,5]. In seismology, fractional models are used to represent intrinsic attenuation and the accompanying velocity dispersion [6,7,8,9,10,11]. Two computational difficulties are characteristic of such models: the fractional operator is non-local in time, so that a direct discretisation requires storing and re-using the whole history, and nonlinear terms complicate the stability analysis.

Spectral methods for linear time-fractional wave and diffusion-wave problems have been studied, e.g., in [12,13,14,15], and a considerable body of work addresses nonlinear time-fractional Klein–Gordon-type equations [16,17,18,19,20,21,22,23,24], and high-order discretisations of related nonlinear fractional models are reported in [25,26]. The memory cost of the Caputo derivative can be reduced by sum-of-exponentials (SOE) or kernel-compression techniques [27,28,29,30], which rest on exponential-sum approximations of the power kernel [31,32].

In this paper we consider a one-dimensional wave equation with Caputo fractional damping of order $\alpha \in (0,1)$ and a cubic restoring term, discretised by Legendre–Gauss–Lobatto (LGL) collocation in space and by a linearly implicit three-level scheme with the L1 formula in time. The contributions are as follows.

A fast L1–SOE implementation in which the kernel $t^{-\alpha }$ is approximated by an exponentially convergent trapezoidal rule with strictly positive nodes and weights, together with an explicit recursion for the history term.

A stability and convergence analysis based on a discrete energy identity and a bound on the L1 memory operator. The argument does not use a sign of the discrete memory form; we show numerically that, for the pairing with the centred difference that arises in the energy method, this form is indefinite, so that it cannot simply be discarded. The resulting time-step restriction is independent of the polynomial degree $N$.

• A machine-checked formalisation in Lean 4 [33] with Mathlib [34] of the algebraic core of the analysis. All assumptions that are not formalised are stated explicitly.
• A seismic application in which the damping coefficient is calibrated from a prescribed quality factor, the discretisation is verified against an exact solution and against a second-order finite-difference discretisation, and the attenuation law recovered from the simulated seismograms is compared with the dispersion relation of the continuous model.
• An extension of the scheme to Caputo’s constitutive model [4], in which the fractional derivative acts on the strain; this model yields a nearly frequency-independent quality factor, as is typical of rocks.

Individually, LGL collocation, the L1 formula and SOE compression are known; the novelty of this work lies in the stability argument that avoids the sign of the discrete memory form (and in the observation that this sign is not available), in its machine-checked formalisation, and in the quantitative validation of the simulated attenuation law against the continuous dispersion relation for two different fractional attenuation models. The cubic term is used as a phenomenological Klein–Gordon-type nonlinearity, as in [16,17,22]; its physical interpretation is discussed in Section 2.2.

The paper is organised as follows. Section 2 states the problem and analyses the attenuation law of the model. Section 3 describes the scheme and the SOE acceleration. Section 4 contains the stability and convergence analysis and its formalisation. Section 5 reports numerical tests, Section 6 the seismic application, and Section 7 and Section 8 contain a discussion and conclusions.

2. Problem Formulation

Let $\mathsf{\Omega }=(0,L)$ and $T>0$. We consider

$${\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}=c^2{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+\gamma \hspace{0.17em}{}^{C}D_t^{\alpha }u+f\left(u\right)+g(x,t), x\in \mathsf{\Omega },\hspace{0.25em}t\in (0,T],$$

(1)

$$u(x,0)=u_0\left(x\right), u_t(x,0)=v_0\left(x\right), u(0,t)=u(L,t)=0,$$

(2)

where $c>0$ is the wave speed, $\gamma <0$ is the damping coefficient, $g$ is a source term and ${}^{C}D_t^{\alpha }$ is the Caputo derivative of order $\alpha \in (0,1)$,

$${}^{C}D_t^{\alpha }u(x,t)={\displaystyle \frac{1}{\mathsf{\Gamma }(1-\alpha )}}{\int }_0^t(t-s{)}^{-\alpha }\hspace{0.17em}{\displaystyle \frac{\partial u}{\partial s}}(x,s)\hspace{0.17em}ds.$$

(3)

The nonlinearity is the cubic restoring force $f\left(u\right)=-{\beta }^2{u}^3$ with $\beta \ge 0$. In the analysis we only use that $f$ is locally Lipschitz with $f\left(0\right)=0$.

2.1. Attenuation Law of the Linear Model

For $f=g=0$ and a plane wave $u=\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathrm{i}\right(\omega t-kx\left)\right)$, the (steady-state) Caputo derivative acts as multiplication by $(\mathrm{i}\omega {)}^{\alpha }$, and (1) yields the dispersion relation

$$c^2{k}^2={\omega }^2-\left|\gamma \right|\hspace{0.17em}{\omega }^{\alpha }\hspace{0.17em}{e}^{\mathrm{i}\alpha \pi /2}.$$

(4)

Defining the quality factor by $Q^{-1}=|\mathrm{I}\mathrm{m}\hspace{0.17em}{k}^2|/\mathrm{R}\mathrm{e}\hspace{0.17em}{k}^2$ we obtain

$$Q^{-1}\left(\omega \right)={\displaystyle \frac{\left|\gamma \right|\hspace{0.17em}{\omega }^{\alpha -2}\mathrm{s}\mathrm{i}\mathrm{n}(\alpha \pi /2)}{1-\left|\gamma \right|\hspace{0.17em}{\omega }^{\alpha -2}\mathrm{c}\mathrm{o}\mathrm{s}(\alpha \pi /2)}}.$$

(5)

Three consequences are important for the interpretation of the model. (i) For $\alpha \to 1$ the damping term becomes the classical viscous term $\gamma u_t$ and $Q=\omega /\left|\gamma \right|$; for $\alpha \to 0$ it becomes $\gamma (u-u_0)$, which does not dissipate energy. For fixed $\gamma$ and $\omega \ge 1\hspace{0.17em}{\mathrm{s}}^{-1}$ (in particular in the seismic frequency band), the attenuation therefore increases with $\alpha$; because the unit of $\gamma$ depends on $\alpha$ (see (iii)), such comparisons at equal $\gamma$ depend on the units and are only indicative. (ii) For weak attenuation, $Q\propto {\omega }^{2-\alpha }$, i.e., the model describes a power-law, frequency-dependent attenuation and not a constant-$Q$ medium; constant-$Q$ models are obtained by placing a fractional derivative in the stress–strain relation [6,7,35]. (iii) The coefficient $\gamma$ is not dimensionless ($\left[\gamma \right]={\mathrm{s}}^{\alpha -2}$) and must be calibrated. Given a target quality factor $Q_0$ at the dominant angular frequency ${\omega }_0$, (5) gives exactly

$$\left|\gamma \right|={\displaystyle \frac{{\omega }_0^{\hspace{0.17em}2-\alpha }}{Q_0\mathrm{s}\mathrm{i}\mathrm{n}(\alpha \pi /2)+\mathrm{c}\mathrm{o}\mathrm{s}(\alpha \pi /2)}},$$

(6)

so that $Q\left({\omega }_0\right)=Q_0$. The order $\alpha$ then controls the frequency dependence of $Q$ around ${\omega }_0$. This calibration is used in Section 6.

2.2. Role of the Cubic Term

The cubic term $-{\beta }^2{u}^3$ is an on-site (Klein–Gordon or Duffing-type) nonlinearity. It should not be confused with the nonlinearity of finite-strain elasticity: the Murnaghan–Landau expansion of the elastic energy [36] leads to a nonlinear flux ${\partial }_x\sigma \left({\partial }_{x}u\right)$ whose leading correction is quadratic in the strain. We therefore regard $f\left(u\right)=-{\beta }^2{u}^3$ as a phenomenological model and as a standard benchmark for numerical methods [16,17,22]. For a wave of amplitude $A_0$ and angular frequency ${\omega }_0$, the relative size of the cubic term is measured by

$$\nu ={\displaystyle \frac{{\beta }^2{A}_0^2}{{\omega }_0^2}}.$$

(7)

For the strain amplitudes of seismic waves away from the source ($\lesssim {10}^{-5}$), $\nu$ is negligible unless $\beta$ is extremely large; the nonlinear seismic computations in Section 6.5 are therefore a parametric study rather than a prediction for a specific rock.

3. Numerical Scheme

3.1. LGL Collocation in Space

Let ${\xi }_0=-1<{\xi }_1<\dots <{\xi }_N=1$ be the LGL nodes, i.e., the zeros of $(1-{\xi }^2)P_{N^{\mathrm{\text{'}}}}\left(\xi \right)$, with quadrature weights $w_j$, and let $x_j=L({\xi }_j+1)/2$. Let $D$ be the first-order Legendre collocation differentiation matrix [37,38] and let

$$A=\left({\displaystyle \frac{2}{L}}{)}^2\right(D^2{)}_{\mathrm{i}\mathrm{n}\mathrm{t}}\in {\mathbb{R}}^{(N-1)\times (N-1)}$$

(8)

be the second-derivative matrix restricted to the interior nodes, which incorporates the homogeneous Dirichlet conditions. For $\mathbf{v}\in {\mathbb{R}}^{N-1}$ we denote by $p_{\mathbf{v}}\in {\mathbb{P}}_N$ the polynomial interpolating $\mathbf{v}$ at the interior nodes and vanishing at $x=0,L$; then $(A\mathbf{v}{)}_j=p_{\mathbf{v}}^{\mathrm{\text{'}}\mathrm{\text{'}}}(x_j)$. We use the discrete inner product and norm

$$⟨\mathbf{v},\mathbf{w}{⟩}_{\omega }={\displaystyle \sum _{j=1}^{N-1}}{\omega }_j{v}_j{w}_j, \parallel \mathbf{v}{\parallel }_{\omega }^2=⟨\mathbf{v},\mathbf{v}{⟩}_{\omega }, {\omega }_j={\displaystyle \frac{L}{2}}{w}_j,$$

(9)

and the stiffness form $a\left(\mathbf{v}\right)=⟨-A\mathbf{v},\mathbf{v}{⟩}_{\omega }$.

3.2. L1 Approximation of the Caputo Derivative

On the uniform grid $t_n=n\tau$, $\tau =T/N_t$, the L1 formula reads [39]

$${}^{C}D_t^{\alpha }u\left(t_n\right)\approx H^n={\displaystyle \frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}}{\displaystyle \sum _{k=0}^{n-1}}{b}_k\left(U^{n-k}-U^{n-k-1}\right), b_k=(k+1{)}^{1-\alpha }-k^{1-\alpha }.$$

(10)

The weights satisfy $b_k>0$ and ${\displaystyle \sum _{k=0}^{n-1}}{b}_k=n^{1-\alpha }$. For $u\in C^2[0,T]$ the truncation error is $O\left({\tau }^{2-\alpha }\right)$. A direct evaluation of (10) costs $O\left(n\right)$ operations per node at step $n$, hence $O\left(N_t^2\right)$ in total, and requires storing all time levels.

3.3. Fast Evaluation by a Sum of Exponentials

The kernel admits the exact representation

$$t^{-\alpha }={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{e}\mathrm{x}\mathrm{p}\left(-t\hspace{0.17em}{e}^y+\alpha y\right)\hspace{0.17em}dy, t>0,$$

(11)

which follows from ${\int }_0^{\mathrm{\infty }}{e}^{-ts}{s}^{\alpha -1}ds=\mathsf{\Gamma }\left(\alpha \right)t^{-\alpha }$ with $s=e^y$. Truncating the integral to $[y_{\mathrm{m}\mathrm{i}\mathrm{n}},y_{\mathrm{m}\mathrm{a}\mathrm{x}}]$ and applying the trapezoidal rule with step $h$ gives [31,32]

$$t^{-\alpha }\approx {\displaystyle \sum _{l=1}^{N_e}}{w}_l\hspace{0.17em}{e}^{-s_{l}t}, s_l=e^{y_l}, w_l={\displaystyle \frac{h\hspace{0.17em}{e}^{\alpha y_l}}{\mathsf{\Gamma }\left(\alpha \right)}}, t\in [\tau ,T].$$

(12)

All nodes and weights are strictly positive by construction. The step $h$ is reduced until the relative error on $[\tau ,T]$, measured on a dense logarithmic grid, is below a prescribed tolerance ${\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}$, and exponentials that are negligible on $[\tau ,T]$ are removed. For ${\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}={10}^{-8}$ and $\alpha =0.5$ this gives $N_e=116$ exponentials for $N_t=64$ and $N_e=131$ for $N_t=32,768$ (on $T=1$), in agreement with the logarithmic growth of $N_e$ with $T/\tau$ [32].

Following [27], the most recent interval is treated exactly by the L1 formula and the history part by (12). With the piecewise linear interpolant of $U$ this gives

$$H_{\mathrm{S}\mathrm{O}\mathrm{E}}^n={\displaystyle \frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}}\left(U^n-U^{n-1}\right)+{\displaystyle \frac{1}{\mathsf{\Gamma }(1-\alpha )}}{\displaystyle \sum _{l=1}^{N_e}}{w}_l\hspace{0.17em}{Z}_l^n,$$

(13)

where the history variables satisfy the recursion

$$Z_l^n=e^{-s_l\tau }{Z}_l^{n-1}+c_l\left(U^{n-1}-U^{n-2}\right), Z_l^1=0, c_l={\displaystyle \frac{e^{-s_l\tau }-e^{-2s_l\tau }}{s_l\tau }}.$$

(14)

Indeed, $Z_l^n={\int }_0^{t_{n-1}}{e}^{-s_l(t_n-s)}\hspace{0.17em}{\mathsf{\Pi }}^{\mathrm{\text{'}}}\left(s\right)\hspace{0.17em}ds$ for the piecewise linear interpolant $\mathsf{\Pi }$ of the time levels, and splitting off the last interval yields (14). The cost per time step is $O\left(N_{e}N\right)$ and the memory $O\left(N_{e}N\right)$, independent of $n$.

3.4. Fully Discrete Scheme

Let ${\mathbf{U}}^n\in {\mathbb{R}}^{N-1}$ approximate $u(x_j,t_n)$ at the interior nodes. The fully discrete scheme is

$${\displaystyle \frac{{\mathbf{U}}^{n+1}-2{\mathbf{U}}^n+{\mathbf{U}}^{n-1}}{{\tau }^2}}=c^{2}A\hspace{0.17em}{\displaystyle \frac{{\mathbf{U}}^{n+1}+{\mathbf{U}}^{n-1}}{2}}+\gamma \hspace{0.17em}{H}^n+f\left({\mathbf{U}}^n\right)+{\mathbf{g}}^n, n=1,\dots ,N_t-1,$$

(15)

where $H^n$ is given by (10) or (13), ${\mathbf{U}}^0=u_0\left(x_j\right)$ and the first step is computed from a Taylor expansion,

$${\mathbf{U}}^1={\mathbf{U}}^0+\tau {\mathbf{v}}_0+{\displaystyle \frac{{\tau }^2}{2}}\left(c^{2}A{\mathbf{U}}^0+f\left({\mathbf{U}}^0\right)+{\mathbf{g}}^0\right),$$

(16)

using that the Caputo term vanishes at $t=0$. Since $H^n$ and $f\left({\mathbf{U}}^n\right)$ involve only known levels, each step requires the solution of

$$\left(I-{\displaystyle \frac{c^2{\tau }^2}{2}}A\right){\mathbf{U}}^{n+1}=2{\mathbf{U}}^n-{\mathbf{U}}^{n-1}+{\displaystyle \frac{c^2{\tau }^2}{2}}A{\mathbf{U}}^{n-1}+{\tau }^2\left(\gamma H^n+f\left({\mathbf{U}}^n\right)+{\mathbf{g}}^n\right).$$

(17)

By Lemma 1 below, the matrix $I-{\displaystyle \frac{c^2{\tau }^2}{2}}A$ is symmetric positive definite with respect to $⟨\cdot ,\cdot {⟩}_{\omega }$ and hence nonsingular for every $\tau >0$; it is factorised once by LU decomposition.

4. Stability and Convergence

We write ${\delta }_t{\mathbf{U}}^n=({\mathbf{U}}^{n+1}-{\mathbf{U}}^n)/\tau$, ${\widehat{\delta }}_t{\mathbf{U}}^n=({\mathbf{U}}^{n+1}-{\mathbf{U}}^{n-1})/\left(2\tau \right)$ and define the discrete energy

$${\mathcal{E}}^n={\displaystyle \frac{1}{2}}\parallel {\delta }_t{\mathbf{U}}^n{\parallel }_{\omega }^2+{\displaystyle \frac{c^2}{4}}\left(a\left({\mathbf{U}}^{n+1}\right)+a\left({\mathbf{U}}^n\right)\right), n\ge 0.$$

(18)

Lemma 1 (Properties of the LGL stiffness matrix). For all $\mathit{v},\mathit{w}\in {\mathbb{R}}^{N-1}$: (a) $⟨-A\mathit{v},\mathit{w}{⟩}_{\omega }={\int }_0^L{p}_{\mathit{v}}^{\text{'}}{p}_{\mathit{w}}^{\text{'}}\hspace{0.17em}dx$; in particular $A$ is symmetric with respect to $⟨\cdot ,\cdot {⟩}_{\omega }$and $a\left(\mathit{v}\right)=\parallel p_{\mathit{v}}^{\text{'}}{\parallel }_{L^2}^2\ge 0$with equality only for $\mathit{v}=0$; (b) $\parallel \mathit{v}{\parallel }_{\omega }^2\le C_P^2\hspace{0.17em}a\left(\mathit{v}\right)$with $C_P^2=3L^2/{\pi }^2$; (c) ${max}_j\left|v_j\right|\le \sqrt{L}\hspace{0.17em}a(\mathit{v}{)}^{1/2}.$

Proof. (a) Since $p_{\mathbf{v}}^{\mathrm{\text{'}}\mathrm{\text{'}}}{p}_{\mathbf{w}}$ has degree $2N-2$ and the LGL rule is exact for polynomials of degree $2N-1$, and $p_{\mathbf{w}}$ vanishes at the end points, $⟨-A\mathbf{v},\mathbf{w}{⟩}_{\omega }=-{\int }_0^L{p}_{\mathbf{v}}^{\mathrm{\text{'}}\mathrm{\text{'}}}{p}_{\mathbf{w}}\hspace{0.17em}dx={\int }_0^L{p}_{\mathbf{v}}^{\mathrm{\text{'}}}{p}_{\mathbf{w}}^{\mathrm{\text{'}}}\hspace{0.17em}dx$. (b) For $p\in {\mathbb{P}}_N$ the discrete and continuous $L^2$ norms satisfy $\parallel p{\parallel }_{L^2}^2\le \sum _j{\omega }_{j}p(x_j{)}^2\le (2+1/N)\parallel p{\parallel }_{L^2}^2$ [37]; combined with the Poincaré inequality $\parallel p{\parallel }_{L^2}\le (L/\pi )\parallel p^{\mathrm{\text{'}}}{\parallel }_{L^2}$ this gives (b). (c) $\left|p_{\mathbf{v}}\right(x\left)\right|=|{\int }_0^x{p}_{\mathbf{v}}^{\mathrm{\text{'}}}\hspace{0.17em}ds|\le \sqrt{L}\hspace{0.17em}\parallel p_{\mathbf{v}}^{\mathrm{\text{'}}}{\parallel }_{L^2}$. □

In our implementation the relative deviation of $WA$ from symmetry, $W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\omega }_j\right)$, is $1.4\times {10}^{-14}$ for $N=48$.

Lemma 2 (Discrete energy identity). Let ${\mathit{U}}^n$satisfy (15) and set ${\mathit{F}}^n=\gamma H^n+f\left({\mathit{U}}^n\right)+{\mathit{g}}^n$. Then, for $n\ge 1,$

$${\mathcal{E}}^n-{\mathcal{E}}^{n-1}=\tau \hspace{0.17em}⟨{\mathbf{F}}^n,{\widehat{\delta }}_t{\mathbf{U}}^n{⟩}_{\omega }.$$

Proof. Take the inner product of (15) with ${\mathbf{U}}^{n+1}-{\mathbf{U}}^{n-1}$. The identities $⟨\mathbf{a}-2\mathbf{b}+\mathbf{d},\mathbf{a}-\mathbf{d}{⟩}_{\omega }=\parallel \mathbf{a}-\mathbf{b}{\parallel }_{\omega }^2-\parallel \mathbf{b}-\mathbf{d}{\parallel }_{\omega }^2$ and, by the symmetry of $A$, $⟨A(\mathbf{a}+\mathbf{d}),\mathbf{a}-\mathbf{d}{⟩}_{\omega }=⟨A\mathbf{a},\mathbf{a}{⟩}_{\omega }-⟨A\mathbf{d},\mathbf{d}{⟩}_{\omega }$ telescope the inertial and stiffness terms; dividing by $2{\tau }^2$ gives the claim. □

Lemma 3 (Bound on the memory operator). Let$n\tau \le T$ and $K_m=T^{1-\alpha }/\Gamma (2-\alpha )$. Then

$$\parallel H^n{\parallel }_{\omega }\le K_m{\mathrm{m}\mathrm{a}\mathrm{x}}_{0\le j<n}\parallel {\delta }_t{\mathbf{U}}^j{\parallel }_{\omega },$$

$$\parallel H_{\mathrm{S}\mathrm{O}\mathrm{E}}^n{\parallel }_{\omega }\le (1+{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}})\hspace{0.17em}{K}_m{\mathrm{m}\mathrm{a}\mathrm{x}}_{0\le j<n}\parallel {\delta }_t{\mathbf{U}}^j{\parallel }_{\omega }.$$

Proof. Since $b_k\ge 0$ and ${\mathbf{U}}^{n-k}-{\mathbf{U}}^{n-k-1}=\tau \hspace{0.17em}{\delta }_t{\mathbf{U}}^{n-k-1}$, the triangle inequality gives $\parallel H^n{\parallel }_{\omega }\le {\tau }^{1-\alpha }\mathsf{\Gamma }(2-\alpha {)}^{-1}\sum _{k<n}{b}_k{\mathrm{m}\mathrm{a}\mathrm{x}}_j\parallel {\delta }_t{\mathbf{U}}^j{\parallel }_{\omega }=t_n^{1-\alpha }\mathsf{\Gamma }(2-\alpha {)}^{-1}{\mathrm{m}\mathrm{a}\mathrm{x}}_j\parallel {\delta }_t{\mathbf{U}}^j{\parallel }_{\omega }$. The fast operator (13) has the same structure with the history weights $b_k$, $k\ge 1$, replaced by integrals of the SOE kernel over the corresponding intervals. Since the SOE kernel is positive and approximates $t^{-\alpha }$ with relative accuracy ${\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}$ on $[\tau ,T]$, these weights lie in $\left[\right(1-{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}})b_k,(1+{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}\left)b_k\right]$. □

Remark 1 (On the sign of the memory term). For$\gamma <0$one would like to discard the memory contribution$\tau \gamma \sum _n⟨H^n,{\widehat{\delta }}_t{\mathit{U}}^n{⟩}_{\omega }$in the summed energy identity as non-positive. With${\mathit{x}}_j={\mathit{U}}^j-{\mathit{U}}^{j-1}$, this sum is a quadratic form${\displaystyle \sum _{n=1}^m}{\displaystyle \sum _{j=1}^n}{b}_{n-j}⟨{\mathit{x}}_j,{\mathit{x}}_n+{\mathit{x}}_{n+1}{⟩}_{\omega }$. The positivity results for the L1 kernel [39,40] concern the pairing with the backward difference ${\mathit{x}}_n$ only; for the centred pairing the form is indefinite: the smallest eigenvalue of its symmetric part is between $-0.24$ and $-0.21$ for all $\alpha \in \{0.1,\dots ,0.9\}$ and $m\in \left\{\mathrm{100,1000,3000}\right\}$ (script exp_memory_form.jl). The analysis below therefore bounds the memory term by Lemma 3 instead of using its sign.

Theorem 1 (Stability). Let ${\mathit{g}}^n=0$, $\alpha \in \left(\mathrm{0,1}\right)$, $N_t\tau \le T$, and assume that $\parallel f\left({\mathit{U}}^n\right){\parallel }_{\omega }\le L_f\parallel {\mathit{U}}^n{\parallel }_{\omega }$for all $n$. Let

$$K={\displaystyle \frac{1}{2}}\left(3\left(2{\gamma }^2{K}_m^2+{\displaystyle \frac{4L_f^2{C}_P^2}{c^2}}\right)+1\right).$$

If $\tau K\le 1/2$, then the solution of (15) with the L1 memory (10) satisfies

$${\mathcal{E}}^m\le (1+4\tau K{)}^m\hspace{0.17em}{\mathcal{E}}^0\le e^{4K{t}_m}\hspace{0.17em}{\mathcal{E}}^0, 0\le m\le N_t.$$

The same holds for the L1–SOE memory (13) with $K_m$replaced by $(1+{\epsilon }_{SOE})K_m.$

Proof. Fix $n\ge 1$ and let $B\ge 0$ with ${\mathcal{E}}^j\le B$ for $j\le n-1$. By (18) and Lemma 1, $\parallel {\delta }_t{\mathbf{U}}^j{\parallel }_{\omega }^2\le 2{\mathcal{E}}^j$ and $a\left({\mathbf{U}}^n\right)\le 4{\mathcal{E}}^{n-1}/c^2$. Lemma 3 gives $\parallel H^n{\parallel }_{\omega }^2\le 2K_m^{2}B$, and the Lipschitz bound with Lemma 1(b) gives $\parallel f\left({\mathbf{U}}^n\right){\parallel }_{\omega }^2\le 4L_f^2{C}_P^{2}B/c^2$, so that $\parallel {\mathbf{F}}^n{\parallel }_{\omega }^2\le K_{1}B$ with $K_1=2K-1$. In Lemma 2 we use the Cauchy–Schwarz inequality, $\parallel {\mathbf{U}}^{n+1}-{\mathbf{U}}^{n-1}{\parallel }_{\omega }\le \tau (\parallel {\delta }_t{\mathbf{U}}^n{\parallel }_{\omega }+\parallel {\delta }_t{\mathbf{U}}^{n-1}{\parallel }_{\omega })$ and Young’s inequality to obtain ${\mathcal{E}}^n-{\mathcal{E}}^{n-1}\le {\displaystyle \frac{\tau }{2}}\left(K_{1}B+{\mathcal{E}}^n+B\right)$, hence $(1-\tau K){\mathcal{E}}^n\le (1+\tau K)B$ and, for $\tau K\le 1/2$, ${\mathcal{E}}^n\le (1+4\tau K)B$. Taking $B={\mathrm{m}\mathrm{a}\mathrm{x}}_{j<n}{\mathcal{E}}^j$ and arguing by induction gives ${\mathcal{E}}^m\le (1+4\tau K{)}^m{\mathcal{E}}^0$. For the SOE variant, Lemma 3 is used with $(1+{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}})K_m$. □

The time-step restriction $\tau K\le 1/2$ depends on $\gamma$, $\alpha$, $T$, $L_f$, $c$ and $L$, but not on $N$: in contrast to explicit schemes, there is no CFL condition, which matters because the largest eigenvalue of $-A$ grows like $N^4$. For $\gamma =0$ and $f=0$, Lemma 2 shows that ${\mathcal{E}}^n$ is exactly conserved for every $\tau$.

Corollary 1 (A priori bound for the cubic nonlinearity). Let $f\left(u\right)=-{\beta }^2{u}^3$, $M>0$, and let $K$be defined with $L_f={\beta }^2{M}^2$. If $(2\sqrt{L}/c)(e^{4KT}{\mathcal{E}}^0{)}^{1/2}\le M$and $\tau K\le 1/2$, then ${max}_j\left|U_j^n\right|\le M$for all $n\le N_t$, and the conclusion of Theorem 1 holds.

Proof. If ${\mathrm{m}\mathrm{a}\mathrm{x}}_j\left|U_j^n\right|\le M$, then $\left|f\right(U_j^n\left)\right|\le {\beta }^2{M}^2\left|U_j^n\right|$ and the Lipschitz assumption holds at level $n$. Conversely, by Lemma 1(c) and (18), ${\mathrm{m}\mathrm{a}\mathrm{x}}_j\left|U_j^n\right|\le \sqrt{L}\hspace{0.17em}a({\mathbf{U}}^n{)}^{1/2}\le (2\sqrt{L}/c\left)\right({\mathcal{E}}^{n-1}{)}^{1/2}$. Induction on $n$, using Theorem 1 up to level $n-1$, gives the claim. □

Corollary 1 is a small-data result; for the defocusing cubic term a modified energy that includes the potential ${\displaystyle \frac{{\beta }^2}{4}}\sum _j{\omega }_j(U_j^n{)}^4$ could remove the smallness assumption, but this is not pursued here.

Theorem 2 (Convergence). Let the solution of (1)–(2) satisfy $u\in C^4\left(\right[0,T];H^m(\Omega \left)\right)$with$m\ge 2$and$\parallel u{\parallel }_{L^{\infty }}\le M/2$, let${\Pi }_N$denote the$H_0^1$-orthogonal projection onto${\mathbb{P}}_N\cap H_0^1\left(\Omega \right)$, and let${\mathit{\theta }}^n={\Pi }_{N}u(\cdot ,t_n){|}_{x_j}-{\mathit{U}}^n$. Then there exist${\tau }_0>0$and$C$, independent of$\tau$and$N$, such that for$\tau \le {\tau }_0$and$N$sufficiently large

$${\mathrm{m}\mathrm{a}\mathrm{x}}_{0\le n\le N_t}\left(\parallel {\delta }_t{\mathsf{\theta }}^n{\parallel }_{\omega }+a({\mathsf{\theta }}^n{)}^{1/2}\right)\le C\left({\tau }^{2-\alpha }+N^{1-m}+{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}\right).$$

Proof. By Lemma 1(a), the collocation scheme is equivalent to a Galerkin method with numerical integration, since $⟨-A\mathbf{v},\mathbf{w}{⟩}_{\omega }={\int }_0^L{p}_{\mathbf{v}}^{\mathrm{\text{'}}}{p}_{\mathbf{w}}^{\mathrm{\text{'}}}\hspace{0.17em}dx$. Hence the error ${\mathsf{\theta }}^n$ satisfies (15) with ${\mathbf{F}}^n=\gamma H^n\left[\mathsf{\theta }\right]+{\mathbf{G}}^n+{\mathbf{R}}^n$, where ${\mathbf{G}}^n=f\left({\mathsf{\Pi }}_N{u}^n\right)-f\left({\mathbf{U}}^n\right)$ and the residual ${\mathbf{R}}^n$ is defined by $⟨{\mathbf{R}}^n,\mathbf{v}{⟩}_{\omega }=r^n\left(p_{\mathbf{v}}\right)$ for the consistency functional $r^n$. Because ${\mathsf{\Pi }}_N$ is the $H_0^1$-orthogonal projection, ${\int }_0^L({\mathsf{\Pi }}_{N}u{)}^{\mathrm{\text{'}}}{p}_{\mathbf{v}}^{\mathrm{\text{'}}}dx={\int }_0^L{u}^{\mathrm{\text{'}}}{p}_{\mathbf{v}}^{\mathrm{\text{'}}}dx$, so that the stiffness term does not contribute to $r^n$. The remaining contributions are: the truncation errors of the second difference and of the averaged stiffness term, $O\left({\tau }^2\right)$ for $u\in C^4$; the L1 truncation error, $O\left({\tau }^{2-\alpha }\right)$ [39]; the projection errors $(u-{\mathsf{\Pi }}_{N}u{)}_{tt}$ and the LGL quadrature errors for $u_{tt}{p}_{\mathbf{v}}$, $f\left(u\right)p_{\mathbf{v}}$ and $g{p}_{\mathbf{v}}$, $O\left(N^{1-m}\right)$ [16,37,38]; and the SOE error, $O\left({\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}\right)$. Since $\parallel p_{\mathbf{v}}{\parallel }_{L^2}\le \parallel \mathbf{v}{\parallel }_{\omega }$, it follows that $\parallel {\mathbf{R}}^n{\parallel }_{\omega }\le \rho :=C({\tau }^{2-\alpha }+N^{1-m}+{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}})$. As long as ${\mathrm{m}\mathrm{a}\mathrm{x}}_j\left|\right({\mathsf{\Pi }}_N{u}^n\left)\right(x_j\left)\right|\le M$ and ${\mathrm{m}\mathrm{a}\mathrm{x}}_j\left|U_j^n\right|\le M$, $\parallel {\mathbf{G}}^n{\parallel }_{\omega }\le 3{\beta }^2{M}^2\parallel {\mathsf{\theta }}^n{\parallel }_{\omega }$. The energy argument of Theorem 1 with the additional residual (formalised as l1_scheme_error_bound, Section 4.1) gives ${\mathcal{E}}_{\theta }^m\le e^{4K{t}_m}({\mathcal{E}}_{\theta }^0+3t_m{\rho }^2)$. The starting step (16) yields ${\mathcal{E}}_{\theta }^0\le C({\tau }^2+N^{1-m}{)}^2$. The bound on $\parallel {\mathbf{U}}^n{\parallel }_{\mathrm{\infty }}$ required for the Lipschitz estimate follows by induction from the error bound and Lemma 1(c) for $\tau$ small and $N$ large. □

Remark 2 (Regularity). The assumption$u\in C^4\left(\right[0,T];H^m)$, needed for the$O\left({\tau }^2\right)$consistency of the three-level scheme, is standard but restrictive for fractional problems [41,42].For (1) with$v_0\ne 0$the Caputo term behaves like$v_0{t}^{1-\alpha }/\Gamma (2-\alpha )$near$t=0$, so that$u_{tt}$is only Hölder continuous and$u_{ttt}\sim t^{-\alpha }$. The resulting initial layer is milder than for subdiffusion. In the test ofSection 5.2 (unforced problem with $v_0\ne 0$) no loss of the temporal rate at$t=T$is observed; for stronger singularities graded temporal meshes [41] can be used.

4.1. Machine-Checked Part of the Analysis

The algebraic core of Lemmas 2 and 3 and Theorems 1 and 2 was formalised in Lean 4 [33] using Mathlib [34]. The discrete space ${\mathbb{R}}^{N-1}$ with $⟨\cdot ,\cdot {⟩}_{\omega }$ is modelled by an arbitrary real inner-product space, the stiffness operator by a linear map, and the L1 operator (10) is defined explicitly. Table 1 lists the formalised statements. The development contains no sorry, and the command #print axioms reports only the standard axioms propext, Classical.choice and Quot.sound. It was checked with Lean v4.35.0-rc2 and a pinned Mathlib commit (Supplementary Materials).

The following ingredients are assumptions of the formal statements and are proved on paper only: the symmetry and semi-definiteness of $A$ and the Poincaré inequality (Lemma 1), the Lipschitz bound on $f$ (Corollary 1), the truncation-error estimates in the proof of Theorem 2 and the SOE perturbation bound in Lemma 3. The formalisation thus certifies the energy argument itself, which is the part of such proofs where sign and index errors are most likely.

5. Numerical Experiments

All computations were performed in Julia 1.12.4 with the packages FastGaussQuadrature 1.3.0, SpecialFunctions 2.9.0 and Plots 1.41.7, on a single core of an Intel Core i5-9400F processor (2.9 GHz, 16 GB RAM); the formal proofs were checked with Lean v4.35.0-rc2 and Mathlib (commit b0df0c4). The complete code and the Lean project are provided in the Supplementary Materials and are archived on Zenodo [43]. The generative AI tool Claude (Anthropic, Claude Opus 5) was used for translation and generation of the manuscript text, and to assist in implementing and testing the SOE solver, in developing the Lean 4 formalisation and in producing the figures; all text, code, proofs and results were checked by the authors (see Acknowledgments). To test accuracy we use (1) with $L=T=c=1$, $\gamma =-1$, $f\left(u\right)=-u^3$ and the manufactured solution $u(x,t)=\mathrm{s}\mathrm{i}\mathrm{n}\left(\pi x\right)(1+t^2)$, for which

$$g(x,t)=\left(2+c^2{\pi }^2(1+t^2)-{\displaystyle \frac{2\gamma \hspace{0.17em}{t}^{2-\alpha }}{\mathsf{\Gamma }(3-\alpha )}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\pi x\right)+{\mathrm{s}\mathrm{i}\mathrm{n}}^3\left(\pi x\right)(1+t^2{)}^3.$$

(19)

Errors are measured in the discrete norm $\parallel \cdot {\parallel }_{\omega }$ at $t=T$. Unless stated otherwise, the L1–SOE memory with ${\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}={10}^{-8}$ is used.

5.1. Temporal Convergence

Table 2 shows the errors for $N=48$, for which the spatial error is negligible. Two error components are present: an $O\left({\tau }^2\right)$ contribution of the three-level scheme and an $O\left({\tau }^{2-\alpha }\right)$ contribution of the L1 formula, whose constant is proportional to $\left|\gamma \right|$. For $\alpha =0.7$ the observed rate increases monotonically towards $2-\alpha =1.3$. For $\alpha =0.5$ the two contributions have opposite signs and cancel near $N_t\approx 700$, which produces the non-monotone behaviour; for larger $N_t$ the rate approaches $1.5$. For $\alpha =0.3$ the $O\left({\tau }^{2-\alpha }\right)$ contribution is smaller than the $O\left({\tau }^2\right)$ one over the whole range, and the observed rate is close to 2; the asymptotic rate $1.7$ is not reached before round-off errors, which are amplified like ${ϵ}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{h}}/{\tau }^2$ by the second difference, become visible at errors of about ${10}^{-9}$.

5.2. Unforced Problem with an Initial Layer

The manufactured solution above is smooth in time. To test the effect of the weak initial singularity discussed in Remark 2, we solve the unforced problem $u_{tt}=u_{xx}+\gamma \hspace{0.17em}{}^{C}D_t^{\alpha }u$, $u_0=\mathrm{s}\mathrm{i}\mathrm{n}\left(\pi x\right)$, with $v_0=0$ and $v_0=\mathrm{s}\mathrm{i}\mathrm{n}\left(\pi x\right)$ ($L=T=c=1$, $\gamma =-1$, $N=16$). Since the exact solution is not available in closed form, and a very fine reference solution would be contaminated by round-off (Section 5.1), we use the self-convergence differences $d\left(N_t\right)=\parallel {\mathbf{U}}_{N_t}\left(T\right)-{\mathbf{U}}_{2N_t}\left(T\right){\parallel }_{\omega }$, which behave like the error for a fixed asymptotic rate. Table 3 shows the results for $v_0=\mathrm{s}\mathrm{i}\mathrm{n}\left(\pi x\right)$. For $\alpha =0.5$ and $0.7$ the observed rates approach $2-\alpha$ as for the smooth solution; for $\alpha =0.3$ the $O\left({\tau }^2\right)$ component again dominates. The results for $v_0=0$ are similar (script exp_initial_layer.jl). Thus the initial layer does not reduce the rate at $t=T$ in this test.

5.3. Spatial Convergence

Figure 1 shows the error as a function of $N$ for $\alpha =0.5$ and $N_t=8192$. The error decreases exponentially until it reaches the temporal error level $3.5\times {10}^{-8}$ (cf. Table 2) at $N=10$. The exponential rate reflects the analyticity of the manufactured solution; for $u\in H^m$ Theorem 2 predicts an algebraic rate $N^{1-m}$.

5.4. Efficiency of the SOE Acceleration

Table 4 compares the CPU time spent on the memory term for the direct L1 formula and the L1–SOE formula ($N=48$, $\alpha =0.5$; best of three runs). The two solutions agree to within $4\times {10}^{-11}$ in the relative discrete norm, far below the discretisation error. The direct evaluation grows quadratically with $N_t$ (the growth beyond a factor of four between the last two rows is due to cache effects), whereas the SOE evaluation grows linearly. The SOE variant becomes faster than the direct one at $N_t\approx 1000$–$2000$ and is $19$ times faster for $N_t=16,384$ and $60$ times faster for $N_t=32,768$. In addition, the SOE variant stores $N_e\approx 130$ history vectors instead of $N_t$.

${\mathit{N}}_{\mathit{t}}{\mathit{N}}_{\mathit{e}}$

5.5. Energy Behaviour

Figure 2 shows the discrete energy for a Gaussian pulse $u_0=\mathrm{e}\mathrm{x}\mathrm{p}(-((x-0.5)/0.05{)}^2)$, $v_0=0$, $L=c=1$, $\gamma =-1$, $N=64$, $N_t=4000$, $T=4$, without and with the cubic term. In the nonlinear case we take ${\beta }^2=400$, for which the cubic term is comparable to the stiffness term, and add the discrete potential ${\displaystyle \frac{{\beta }^2}{8}}(\parallel {\mathbf{U}}^{n+1}{\parallel }_{L_{\omega }^4}^4+\parallel {\mathbf{U}}^n{\parallel }_{L_{\omega }^4}^4)$ to ${\mathcal{E}}^n$. At $T=4$ the energy is $0.911$, $0.584$ and $0.144$ of its initial value for $\alpha =0.2$, $0.5$ and $0.8$ in the linear case, and $0.909$, $0.582$ and $0.142$ with the cubic term; the nonlinearity changes the energy balance only slightly. For this dimensionless problem the dominant frequencies satisfy $\omega \gg 1$, so the stronger decay for larger $\alpha$ agrees with Section 2.1. The decay is not strictly monotone: on individual steps the energy may increase by up to $8\times {10}^{-5}{\mathcal{E}}^0$, which is consistent with Remark 1 and with Theorem 1, which only asserts a bound.

6. Application to Seismic Wave Attenuation

6.1. Set-Up and Verification

We consider a 10 km profile with $c=3000$ m/s and an initial Ricker wavelet $u_0\left(x\right)=(1-2a)e^{-a}$, $a=\left(\pi f_0\right(x-x_0)/c{)}^2$, with $f_0=10$ Hz, $x_0=5$ km and $v_0=0$. The amplitude is normalised to one. The wavelet splits into two pulses travelling in opposite directions, which reach the boundaries after $1.67$ s; with $T=1.5$ s the results are free of boundary reflections. The parameters are listed in Table 5. Receiver traces are obtained by barycentric interpolation of the LGL solution.

The discretisation parameters were chosen by comparison with the exact elastic solution $u={\displaystyle \frac{1}{2}}\left(u_0\right(x-ct)+u_0(x+ct\left)\right)$ ($\gamma =0$, $f=0$) at receivers at 6, 7 and 8 km (Table 6). With $N=128$ the pulse is severely under-resolved (relative error 77%) and the numerical solution exhibits spurious dispersive wakes; for $N\ge 384$ the spatial error is negligible and the remaining error is temporal and decreases at second order. The choice $N=768$, $N_t=6000$ gives a relative error of $0.5\mathrm{\%}$.

6.2. Attenuation in Three Media

We fix $\alpha =0.5$ and calibrate $\gamma$ by (6) to $Q_0=20$, $60$ and $300$ at $f_0=10$ Hz, illustrative values for strongly attenuating unconsolidated or hydrate-bearing sediments, sandstone and crystalline rock, respectively; this gives $\gamma =-33.54$, $-11.55$ and $-2.34{\mathrm{s}}^{-1.5}$. Figure 3 shows the wavefield at $t=0.5$ s and Figure 4 the peak amplitude as a function of the propagation distance, normalised by the elastic peak. The dashed lines show the plane-wave prediction $\mathrm{e}\mathrm{x}\mathrm{p}(-\pi f_{0}d/(c{Q}_0\left)\right)$ at the dominant frequency. The agreement is close for $Q_0=60$ and $300$; for $Q_0=20$ the deviations reflect the broadband nature of the pulse and the frequency dependence of $Q$, which shifts the dominant frequency during propagation. Figure 5 shows the space–time wavefield for $Q_0=60$.

6.3. Validation of the Attenuation Law

To test whether the discretisation reproduces the attenuation law of the continuous model, we calibrate $\gamma$ to $Q_0=30$ at 10 Hz for $\alpha =0.2$, $0.5$ and $0.8$ and record the right-going pulse at $x_1=6$ km and $x_2=8$ km. The attenuation coefficient is estimated from the spectral ratio, $\widehat{a}\left(f\right)=-\mathrm{l}\mathrm{n}\left(\left|\widehat{u}\right(x_2,f\left)\right|/\left|\widehat{u}\right(x_1,f\left)\right|\right)/(x_2-x_1)$, and converted into $\widehat{Q}\left(f\right)=\mathrm{R}\mathrm{e}\hspace{0.17em}k/\left(2\widehat{a}\right)$, where $k$ is the root of (4) with $\mathrm{I}\mathrm{m}\hspace{0.17em}k<0$. Figure 6 compares $\widehat{Q}$ with $Q_{\mathrm{e}\mathrm{f}\mathrm{f}}=\mathrm{R}\mathrm{e}\hspace{0.17em}k/(-2\hspace{0.17em}\mathrm{I}\mathrm{m}\hspace{0.17em}k)$. At 10 Hz the measured values are $29.98$, $30.00$ and $30.01$ (theory: $30.01$), and the maximal relative deviation between 5 and 20 Hz is $0.48\mathrm{\%}$, $0.09\mathrm{\%}$ and $0.02\mathrm{\%}$ for $\alpha =0.2$, $0.5$ and $0.8$. Thus the scheme reproduces the attenuation law, including its frequency dependence $Q\propto {\omega }^{2-\alpha }$, and the transient memory effects of the Caputo initial-value problem are negligible at these offsets.

6.4. Influence of the Fractional Order

Although the three media in Section 6.3 have the same $Q$ at 10 Hz, their seismograms differ considerably (Figure 7). For $\alpha =0.2$ the quality factor varies from about 4 at 4 Hz to about 130 at 22 Hz (Figure 6), so that low frequencies are strongly suppressed and the pulse becomes narrower-band and oscillatory, with a pronounced trailing wavelet and a phase shift due to velocity dispersion; for $\alpha =0.8$ the frequency dependence is weaker and the pulse shape is closer to the elastic one. Hence $\alpha$ is not determined by the value of $Q$ at a single frequency: it has to be constrained by the observed frequency dependence of attenuation.

6.5. Effect of the Cubic Term

Finally, Figure 8 shows seismograms at $x=8$ km for the sandstone model ($Q_0=60$, $\alpha =0.5$) with ${\beta }^2=\nu {\omega }_0^2$ and $\nu =0$, $0.1$ and $0.3$ (see (7) with $A_0=1$). The relative ${\mathcal{l}}^2$ difference from the linear seismogram is $0.31$ for $\nu =0.1$ and $0.83$ for $\nu =0.3$. The cubic term shifts the main phase to earlier times, reduces its period and, for $\nu =0.3$, generates a slowly decaying tail. As discussed in Section 2.2, such values of $\nu$ require a very strong on-site nonlinearity; the computation illustrates the robustness of the scheme in a moderately nonlinear regime rather than a quantitative prediction.

6.6. Comparison with a Finite-Difference Discretisation

To quantify the benefit of the spectral discretisation, we replace the LGL matrix $A$ by the standard second-order finite-difference matrix on a uniform grid with $M$ intervals and keep the time integrator unchanged. For the elastic test of Section 6.1 with $N_t=24,000$ (temporal error about $3\times {10}^{-4}$), Table 7 and Figure 9 compare the relative receiver errors and the CPU times. The LGL scheme reaches an error of ${10}^{-3}$ with $N=320$ unknowns in 1.3 s, whereas the finite-difference scheme needs $M=16,000$ intervals and 11.1 s for the same accuracy; with $N=384$ the LGL error reaches the temporal error level. The finite-difference matrix is sparse, so its cost per step is lower, but the number of unknowns required for a given accuracy is larger by a factor of about 50. The advantage of the spectral method grows with the required accuracy and with the number of wavelengths in the domain.

6.7. Extension to Caputo’s Constitutive Model

As discussed in Section 2.1, the damping model (1) yields $Q\propto {\omega }^{2-\alpha }$, whereas the attenuation of rocks is often nearly frequency independent. A classical model with this property is Caputo’s constitutive law [1,4], in which the fractional derivative acts on the strain:

$${\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}=c^2\left({\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+\eta \hspace{0.17em}{}^{C}D_t^{\alpha }{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}\right), \eta >0.$$

(20)

For a plane wave, the complex modulus is $1+\eta (\mathrm{i}\omega {)}^{\alpha }$ and

$$Q^{-1}\left(\omega \right)={\displaystyle \frac{\eta \hspace{0.17em}{\omega }^{\alpha }\mathrm{s}\mathrm{i}\mathrm{n}(\alpha \pi /2)}{1+\eta \hspace{0.17em}{\omega }^{\alpha }\mathrm{c}\mathrm{o}\mathrm{s}(\alpha \pi /2)}}, \eta ={\displaystyle \frac{{\omega }_0^{-\alpha }}{Q_0\mathrm{s}\mathrm{i}\mathrm{n}(\alpha \pi /2)-\mathrm{c}\mathrm{o}\mathrm{s}(\alpha \pi /2)}},$$

(21)

where the second relation calibrates $Q\left({\omega }_0\right)=Q_0$. Since $Q\propto {\omega }^{-\alpha }$ for weak attenuation, small $\alpha$ gives a nearly constant $Q$. The scheme of Section 3 carries over with one modification: because the constitutive term contains the unbounded operator $A$, it is averaged like the stiffness term, $c^2\eta A(H^{n+1}+H^{n-1})/2$, and the local L1 part of $H^{n+1}$ is treated implicitly. The system matrix becomes $I-{\displaystyle \frac{c^2{\tau }^2}{2}}\left(1+\eta {\tau }^{-\alpha }/\mathsf{\Gamma }(2-\alpha )\right)A$, which is again constant and symmetric positive definite with respect to $⟨\cdot ,\cdot {⟩}_{\omega }$; the history part is evaluated with the same SOE recursion. The stability analysis of Section 4 does not cover this variant, and its stability is observed numerically.

Figure 10 shows the quality factor measured from simulated seismograms (set-up of Section 6.3, $N=768$, $N_t=6000$) for $Q_0=30$ at 10 Hz. For $\alpha =0.1$, $Q$ varies only from 31.7 at 5 Hz to 28.4 at 20 Hz, and the measured values agree with (21) to within $0.02\mathrm{\%}$ between 5 and 20 Hz for both $\alpha =0.1$ and $0.3$; direct L1 and L1–SOE solutions agree to within $2\times {10}^{-9}$. The seismograms show two differences from the damping model with the same $Q$ at 10 Hz. First, the pulse retains a Ricker-like, broadened shape without the oscillatory trailing wavelet of the damping model. Second, it arrives earlier: the modulus $1+\eta (\mathrm{i}\omega {)}^{\alpha }$ exceeds one, so that the phase velocity at 10 Hz is about $1.13c$ for $\alpha =0.1$ and $1.03c$ for $\alpha =0.3$. In applications, $c$ should therefore be chosen such that the phase velocity at $f_0$ matches the observed velocity, in addition to the calibration of $\eta$.

7. Discussion

The analysis in Section 4 differs from arguments in which the discrete memory term is discarded as dissipative. Remark 1 shows that the corresponding quadratic form is indefinite for the centred pairing of the three-level scheme, and Figure 2 confirms that the discrete energy is not strictly monotone. Bounding the memory term instead (Lemma 3) leads to an exponential factor $e^{4KT}$ in the stability bound, and $K$ grows like ${\gamma }^2$. For the test problems of Section 5 the constants are moderate ($\tau K\le 0.02$). For the seismic runs, however, the calibrated coefficients are large: for $\alpha =0.5$ and $Q_0=300$, $60$, $30$ and $20$ one obtains $\tau K\approx 0.008$, $0.19$, $0.74$ and $1.6$, and for $\alpha =0.2$, $Q_0=30$ even $\tau K\approx 47$. The step-size condition of Theorem 1 is therefore violated in some of the simulations of Section 6, and where it holds the bound $e^{4KT}$ is far too large to be informative. The stability of these computations is thus an empirical observation, supported by the verification against the exact solution and by the agreement of the measured attenuation with the dispersion relation, and not a consequence of Theorem 1. A sharper analysis that exploits the dissipativity of the continuous operator, for instance by pairing the memory term with a backward difference in a modified scheme, is an interesting direction.

The formalisation covers the energy argument, i.e., the part of the proof in which algebraic errors are most likely to go unnoticed. The approximation-theoretic ingredients (LGL quadrature, projection and truncation error estimates, SOE approximation) are assumed. Formalising them would require substantial additional libraries.

The physical interpretation of the model requires care. The damping term $\gamma \hspace{0.17em}{}^{C}D_t^{\alpha }u$ produces a power-law attenuation with $Q\propto {\omega }^{2-\alpha }$ and not a constant $Q$; for fixed $\gamma$ the attenuation grows with $\alpha$, and only a calibration such as (6) makes the parameters comparable across media. For nearly constant-$Q$ modelling, the fractional derivative should be placed in the constitutive relation [6,7,8,35]. Section 6.7 shows that the LGL discretisation and the SOE compression carry over to Caputo’s constitutive model and reproduce its attenuation law; extending the energy analysis to this model, whose memory term involves the unbounded operator $A$, remains open. The resolution study in Section 6.1 shows that under-resolved spectral simulations produce dispersive wakes that can be mistaken for physical dispersion; convergence studies should therefore accompany applications of fractional wave models.

Section 6.6 shows that, for a propagation distance of several wavelengths, the spectral discretisation needs about 50 times fewer unknowns than a second-order finite-difference scheme and is about eight times faster for the same accuracy. The main limitations are the restriction to one space dimension, the small-data assumption in Corollary 1, the use of Dirichlet boundaries on a large domain instead of absorbing boundaries, and the reduced temporal order $2-\alpha$ of the L1 formula. Natural extensions are tensor-product LGL or spectral-element discretisations in two and three dimensions, perfectly matched layers or absorbing boundary conditions [44,45,46], graded meshes [41], and second-order formulas such as ${\mathrm{L2-1}}_{\sigma }$ [17].

8. Conclusions

We have presented an LGL spectral collocation method with a fast L1–SOE treatment of the Caputo damping term for a one-dimensional nonlinear fractional wave equation. The scheme is linearly implicit, requires one LU factorisation, and its memory cost grows only logarithmically with the number of time steps. We proved stability under a time-step restriction independent of the polynomial degree and an error bound of order ${\tau }^{2-\alpha }+N^{1-m}+{\epsilon }_{\mathrm{S}\mathrm{O}\mathrm{E}}$, without relying on a sign of the discrete memory form, and we machine-checked the energy argument in Lean 4. The numerical experiments confirm the temporal rates, spectral accuracy in space and speed-ups of the memory evaluation of up to 60 times. In the seismic application the damping coefficient was calibrated from a prescribed quality factor, and the quality factor recovered from simulated seismograms agreed with the dispersion relation to within 0.5%. Compared with second-order finite differences, the spectral scheme required about 50 times fewer unknowns and was about eight times faster at equal accuracy, and the method carries over to Caputo’s constitutive model, whose nearly constant quality factor it reproduces to within 0.02%.