A Fast Chebyshev Spectral Collocation Method for a Coupled System of Nonlinear Klein–Gordon Equations with Caputo Fractional Memory

1. Introduction

Nonlinear time-fractional partial differential equations (FPDEs) have become a cornerstone in mathematical physics, where fractional integrodifferential operators model anomalous diffusion, viscoelastic damping, and dispersive wave propagation in heterogeneous media [1,2,15]. Replacing classical wave models with fractional hyperbolic models allows one to capture “memory effects” — where the current state of the system depends on its entire historical trajectory [3].

The nonlinear Klein–Gordon equation is fundamental in relativistic quantum mechanics and in condensed matter physics for modeling dispersive solitary waves (solitons). Its extension to the time-fractional domain — incorporating Caputo-type memory integrals into the wave equation — introduces non-trivial long-range temporal interactions. Due to the non-local nature of fractional operators, exact analytical solutions are rare and generally limited to linear, uncoupled equations [4]. Developing stable and efficient numerical methods for such systems is therefore an active area of computational mathematics.

Significant progress has been made for single-component fractional wave equations. Sun and Wu [5] established the $L1$-type finite difference scheme achieving $\mathcal{O}(\Delta t^{2-\alpha })$ accuracy for the Caputo derivative. Berdyshev et al. [18] investigated non-local boundary problems for parabolic–hyperbolic equations with Riemann–Liouville fractional operators. Baigereyev et al. [19] developed convergent numerical methods for fractional models of flow in porous media. To improve spatial accuracy, Bhrawy et al. [6] and Zaky [7] introduced spectral collocation methods with Chebyshev and Legendre polynomials, achieving exponential spatial convergence $\mathcal{O}\left(e^{-cN}\right)$. Baishemirov et al. [20] proposed numerical methods for convection–diffusion equations with Caputo fractional derivatives, incorporating memory effects.

However, extending to coupled nonlinear fractional systems introduces three significant computational and analytical challenges:

• Nonlinear coupling: Coupled systems with cubic self-interactions $u^3$ and inter-component coupling $\eta u{v}^2$ require fully implicit treatments at each time step, typically via Newton–Raphson iterations with Jacobian assembly [8,9].
• Dense memory convolution: Evaluating the Caputo memory integral at time $t_n$ requires $\mathcal{O}\left(n\right)$ operations, leading to $\mathcal{O}\left(N_t^2\right)$ total complexity and $\mathcal{O}(N_t·N)$ memory storage, which becomes prohibitive for long-time simulations [10].
• Coupled stability and convergence analysis: Establishing rigorous a priori stability for systems with two distinct fractional memory orders $\alpha \ne \beta$ and inter-component cubic coupling requires controlling the cross-energy exchange between components through a unified discrete energy functional — a non-trivial extension of scalar fractional stability theory that, to our knowledge, has not been addressed in the existing literature.

To address the dense memory bottleneck, early approaches used the “short-memory principle” by truncating the historical tail, at the cost of convergence accuracy [4]. More recently, fast memory algorithms based on Sum-of-Exponentials (SOE) approximations of the fractional kernel, using Gauss–Laguerre quadratures or Talbot contours, have been developed by Jiang et al. [11] and Baffet & Hesthaven [10], reducing the temporal complexity to $\mathcal{O}\left(N_t{N}_{\exp }\right)$. However, implementing dual-SOE arrays for coupled nonlinear systems with disparate fractional orders remains largely unexplored.

Comparison with existing approaches. Table 1 summarises the key differences between the present work and the most closely related studies.

Specifically: (i) Sun and Wu [5] and Liao et al. [13] developed $L1$-type finite difference schemes for scalar subdiffusion and wave equations, achieving $\mathcal{O}(\Delta t^{2-\alpha })$ temporal accuracy but without spectral spatial discretization or nonlinear coupling. (ii) Bhrawy et al. [6] and Zaky [7] applied Chebyshev/Legendre spectral collocation to scalar nonlinear fractional Klein–Gordon equations, but did not consider coupled systems, disparate orders $\alpha \ne \beta$, or the $\mathcal{O}\left(N_t^2\right)$ memory bottleneck. (iii) Jiang et al. [11] introduced the SOE acceleration for linear fractional diffusion equations, and recently Shiyapov et al. [12] developed high-order structure-preserving spectral schemes with fast memory for nonlocal nonlinear diffusion equations; we extend this to a dual-SOE architecture handling two independent fractional orders within a coupled nonlinear Klein–Gordon system. (iv) Among coupled fractional systems, Li, Huang and Zhao [21] proposed a fast conservative finite difference scheme for the coupled space-fractional Klein–Gordon–Schrödinger equations, and Hu and Cai [22] developed structure-preserving spectral Galerkin methods for the same system. However, both works address space-fractional Riesz derivatives (not time-fractional Caputo memory), couple KG with Schrödinger (not KG with KG), use a single fractional order, and do not employ SOE acceleration. (v) Berdyshev et al. [18] and Baigereyev et al. [19] investigated fractional parabolic–hyperbolic and filtration models, establishing existence and convergence results for related problems, but without spectral discretization or SOE memory compression. (vi) To our knowledge, no existing work combines all of the following elements simultaneously: a coupled nonlinear Klein–Gordon system with two distinct time-fractional Caputo orders, Chebyshev spectral collocation, an IMEX Newmark temporal scheme with second-order extrapolation, dual-SOE memory acceleration, and a rigorous coupled stability–convergence analysis via a discrete Gronwall inequality.

Contributions. In this work, we propose a fast linearly extrapolated implicit Chebyshev spectral collocation method for a coupled nonlinear fractional Klein–Gordon system with Caputo-type memory integrals. The method is developed in one spatial dimension as a rigorous methodological foundation; the algorithmic architecture — temporal integration, IMEX linearisation, and dual-SOE memory compression — is dimension-independent by construction, and we provide a concrete extension pathway to $d\ge 2$ in Section 5.1. Our contributions are:

• An implicit Newmark temporal averaging (${\beta }_{\mathrm{NM}}=1/4$) of the spatial operator [14] that ensures stability for any $\Delta t>0$ at fixed spatial resolution N.
• A second-order linear extrapolation $U^{\ast }=2{\mathbf{U}}^n-{\mathbf{U}}^{n-1}$ that eliminates the need for Newton–Raphson iterations without degrading the $\mathcal{O}(\Delta t^{min(2-\alpha ,2-\beta )})$ temporal accuracy dictated by the L1 memory approximation.
• A dual-SOE acceleration accommodating disparate fractional orders ($\alpha \ne \beta$), reducing the $\mathcal{O}\left(N_t^2\right)$ memory bottleneck.
• A rigorous coupled stability–convergence analysis (Theorems 1–2) via a unified discrete energy method and a discrete Gronwall inequality, treating both components and their inter-component coupling simultaneously. The proof architecture is dimension-independent: only the Chebyshev inverse-inequality constant changes with the spatial dimension d.

1.1. Physical Derivation of the Coupled System

The governing equations arise from Newton’s second law applied to two interacting elastic waveguides whose internal damping exhibits long-range temporal correlations. We summarize the five-step derivation and then give a concrete physical example.

Step 1. Classical elastic wave equations. Consider two one-dimensional elastic waveguides with displacement fields $u(x,t)$ and $v(x,t)$ propagating in $\Omega =(0,L)$. The i-th waveguide has mass density ${\rho }_i$, elastic stiffness $E_i$, and wave speed $c_i=\sqrt{E_i/{\rho }_i}$. From Hooke’s law and the balance of linear momentum, the undamped wave equations are:

$$u_{tt}-c_1^2{u}_{xx}=f_1(x,t),v_{tt}-c_2^2{v}_{xx}=f_2(x,t).$$

(1)

Step 2. Fractional memory damping. In real materials — porous rock, biological tissue, polymers — each material point additionally experiences a hereditary drag body force whose intensity depends on the entire history of the velocity field through a power-law kernel [2]. Let $\alpha \in (0,1)$ be the memory order for the first waveguide and $\beta \in (0,1)$ for the second. Adding the fractional body damping to each equation in (1):

$$\begin{array}{cc}\hfill u_{tt}-c_1^2{u}_{xx}+{\sigma }_1{\int }_0^t{\kappa }_{\alpha }(t-\tau )u_{\tau }(x,\tau )d\tau & =f_1,\hfill \\ \hfill v_{tt}-c_2^2{v}_{xx}+{\sigma }_2{\int }_0^t{\kappa }_{\beta }(t-\tau )v_{\tau }(x,\tau )d\tau & =f_2,\hfill \end{array}$$

(2)

where ${\sigma }_1,{\sigma }_2>0$ are the damping strengths and ${\kappa }_{\gamma }\left(t\right)=t^{-\gamma }/\Gamma (1-\gamma )$ is the weakly singular memory kernel. Each memory integral is precisely the Caputo fractional derivative ^C$D_t^{\alpha }u$ (respectively ^C$D_t^{\beta }v$), scaled by ${\sigma }_i$.

Remark. This model is distinct from the fractional Kelvin–Voigt viscoelastic model, in which leads, via $\rho u_{tt}={\sigma }_x$ and $\epsilon =u_x$, to a memory integral acting on $u_{xx\tau }$. In Equation (2) the damping is a body force proportional to velocity, modelling energy loss at the material point level (e.g., friction with a surrounding porous matrix).

Step 3. Nonlinear Klein–Gordon restoring forces. When the displacement amplitude is not small, the linear restoring terms $m_{i}u$ and $m_{i}v$ are augmented by cubic self-interactions arising from quartic (${\varphi }^4$) potentials:

$$W_1\left(u\right)=\frac{m_1}{2}{u}^2+\frac{k_1}{4}{u}^4,W_2\left(v\right)=\frac{m_2}{2}{v}^2+\frac{k_2}{4}{v}^4.$$

(3)

The restoring forces $W_1^{\prime }\left(u\right)=m_{1}u+k_1{u}^3$ and $W_2^{\prime }\left(v\right)=m_{2}v+k_2{v}^3$ are added to the left-hand sides of (2).

Step 4. Inter-component coupling. If the two waveguides interact through a symmetric quartic coupling potential

$$W_{\mathrm{int}}(u,v)=\frac{\eta }{2}{u}^2{v}^2,$$

(4)

the partial derivatives $\partial W_{\mathrm{int}}/\partial u=\eta u{v}^2$ and $\partial W_{\mathrm{int}}/\partial v=\eta v{u}^2$ provide the coupling terms in the respective equations.

Step 5. Assembly of the coupled system. Combining Steps 1–4 for both components and appending the initial and boundary conditions, we arrive at the complete coupled system:

(5)

which is precisely the system (6) stated formally below. Both equations retain the classical wave operator ${\partial }_{tt}$, preserving the hyperbolic character (finite propagation speed, two initial conditions per component). The key distinction from models that replace $u_{tt}$ by ^C$D_t^{\alpha }u$ is that Equation (5) keep the integer-order inertia term; the fractional memory acts as hereditary damping, not as a replacement for ${\partial }_{tt}$.

Figure 1. (a) Comparison of fractional memory kernels ${\kappa }_{\alpha }\left(t\right)=t^{-\alpha }/\Gamma (1-\alpha )$ for several values of $\alpha \in (0,1)$ with the classical exponential decay $e^{-t}$. The power-law tails decay much more slowly, encoding long-range memory. (b) Schematic of two coupled viscoelastic waveguides. Each component $u,v$ propagates in $\Omega =(0,L)$ with its own fractional memory kernel (${\kappa }_{\alpha },{\kappa }_{\beta }$), and the dashed arrows represent the inter-component nonlinear coupling $\eta u{v}^2$

Figure 1. (a) Comparison of fractional memory kernels ${\kappa }_{\alpha }\left(t\right)=t^{-\alpha }/\Gamma (1-\alpha )$ for several values of $\alpha \in (0,1)$ with the classical exponential decay $e^{-t}$. The power-law tails decay much more slowly, encoding long-range memory. (b) Schematic of two coupled viscoelastic waveguides. Each component $u,v$ propagates in $\Omega =(0,L)$ with its own fractional memory kernel (${\kappa }_{\alpha },{\kappa }_{\beta }$), and the dashed arrows represent the inter-component nonlinear coupling $\eta u{v}^2$

Example (seismic waveguide). Consider two adjacent geological layers with distinct viscoelastic properties. The horizontal displacement in the i-th layer satisfies Equation (2) with parameters $(c_i,{\alpha }_i,{\sigma }_i)$ fitted to laboratory creep data. For a sandstone layer one typically finds $\alpha \approx 0.3$–$0.5$[2]; for a clay layer the memory is stronger, $\alpha \approx 0.6$–$0.8$. The inter-layer friction at the contact surface is modelled by the coupling term $\eta u{v}^2$ with $\eta$ calibrated from displacement correlation measurements. With $c_1=1500$ m/s, $c_2=800$ m/s, $\alpha =0.4$, $\beta =0.7$, ${\sigma }_1={\sigma }_2=1$, $m_1=m_2=1$, $k_1=k_2=1$, and $\eta =1$, the resulting system is precisely of the form (6). The numerical experiments in Section 4.5 simulate a related scenario with counter-propagating solitary pulses.

We consider the following coupled system. Let $\Omega =(0,L)$ be a bounded one-dimensional spatial domain with boundary $\partial \Omega$.

Given data:

• The fractional memory orders $\alpha ,\beta \in (0,1)$ characterizing the anomalous Caputo memory damping in each component.
• The memory strength coefficients ${\sigma }_1,{\sigma }_2>0$.
• The characteristic wave propagation velocities $c_1,c_2>0$.
• The linear restoring mass/stiffness coefficients $m_1,m_2\ge 0$.
• The cubic self-interaction strength parameters $k_1,k_2>0$ and the inter-component nonlinear coupling parameter $\eta >0$.
• The continuous external source forcing functions $f_1(x,t)$ and $f_2(x,t)$ defined on $\Omega \times (0,T]$.
• The initial displacement profiles ${\varphi }_1\left(x\right),{\varphi }_2\left(x\right)$ and the corresponding initial velocity fields ${\psi }_1\left(x\right),{\psi }_2\left(x\right)$.
• The prescribed inhomogeneous Dirichlet boundary data $g_1\left(t\right),g_2\left(t\right),h_1\left(t\right),h_2\left(t\right)$ defined on $\partial \Omega \times [0,T]$.

Problem (Find): Determine the real-valued scalar displacement fields $u(x,t)$ and $v(x,t)$ such that they satisfy the following coupled initial-boundary value system for all $x\in \Omega$ and $t\in (0,T]$:

$$\left\{\begin{array}{c}{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}-c_1^2\Delta u+{\sigma }_1{\int }_0^t{\kappa }_{\alpha }(t-\tau )\frac{\partial u}{\partial \tau }(x,\tau )d\tau }\hfill \\ {\displaystyle +m_{1}u+k_1{u}^3+\eta u{v}^2=f_1,}\hfill & x\in \Omega ,t\in (0,T],\hfill \\ {\displaystyle \frac{{\partial }^{2}v}{\partial t^2}-c_2^2\Delta v+{\sigma }_2{\int }_0^t{\kappa }_{\beta }(t-\tau )\frac{\partial v}{\partial \tau }(x,\tau )d\tau }\hfill \\ {\displaystyle +m_{2}v+k_2{v}^3+\eta v{u}^2=f_2,}\hfill & x\in \Omega ,t\in (0,T],\hfill \\ u(x,0)={\varphi }_1\left(x\right),{\partial }_{t}u(x,0)={\psi }_1\left(x\right),\hfill & x\in \Omega ,\hfill \\ v(x,0)={\varphi }_2\left(x\right),{\partial }_{t}v(x,0)={\psi }_2\left(x\right),\hfill & x\in \Omega ,\hfill \\ u(0,t)=g_1\left(t\right),u(L,t)=h_1\left(t\right),\hfill & t\in [0,T],\hfill \\ v(0,t)=g_2\left(t\right),v(L,t)=h_2\left(t\right),\hfill & t\in [0,T],\hfill \end{array}\right.$$

where $\Delta ={\displaystyle \frac{{\partial }^2}{\partial x^2}}$ is the one-dimensional spatial Laplacian. The convolution integrals represent Caputo-type fractional memory effects with the weakly singular power-law kernels

$${\kappa }_{\gamma }\left(t\right)={\displaystyle \frac{t^{-\gamma }}{\Gamma (1-\gamma )}},\gamma \in \{\alpha ,\beta \}\subset (0,1),$$

(6)

where ${\sigma }_1,{\sigma }_2>0$ are the memory strength coefficients and $\Gamma (·)$ denotes the Euler Gamma function. The memory integral thus reduces to the Caputo fractional derivative of u (respectively v) of order $\gamma \in (0,1)$, scaled by ${\sigma }_i$. Crucially, the leading temporal operator remains the classical second-order derivative${\partial }_{tt}$, ensuring the hyperbolic (wave) character of the system. The fractional memory terms act as viscoelastic damping, not as replacements for ${\partial }_{tt}$. Two initial conditions per component (${\varphi }_i$ and ${\psi }_i$) are therefore required, consistent with the second-order nature of the equations.

The primary obstacle in simulating Equation (6) is the strongly singular memory tail of the Caputo convolution integral (6), which requires continuously accumulating and evaluating the convolution over all previous time steps.

2. Numerical Methodology

We partition the time domain $t\in [0,T]$ into $N_t$ uniform steps $\Delta t=T/N_t$, with $t_n=n\Delta t$. The continuous physical domain $[0,L]$ is discretized using $N+1$ shifted Chebyshev–Gauss–Lobatto (CGL) collocation points $x_m$. The continuous problem (6) is translated into the following fully discrete scheme advancing from $t_n$ to $t_{n+1}$:

$$\left\{\begin{array}{c}{\displaystyle \frac{{\mathbf{U}}^{n+1}-2{\mathbf{U}}^n+{\mathbf{U}}^{n-1}}{\Delta t^2}-c_1^2{\mathbf{D}}^2{\overline{\mathbf{U}}}^n+{\sigma }_1{a}_0^{\left(\alpha \right)}{\mathbf{H}}_u^n+m_1{\mathbf{U}}^n+{\mathbf{N}}_u={\mathbf{f}}_1^n,}\hfill \\ {\displaystyle \frac{{\mathbf{V}}^{n+1}-2{\mathbf{V}}^n+{\mathbf{V}}^{n-1}}{\Delta t^2}-c_2^2{\mathbf{D}}^2{\overline{\mathbf{V}}}^n+{\sigma }_2{a}_0^{\left(\beta \right)}{\mathbf{H}}_v^n+m_2{\mathbf{V}}^n+{\mathbf{N}}_v={\mathbf{f}}_2^n,}\hfill \end{array}\right.$$

(7)

where we use the Newmark averaging ${\overline{\mathbf{U}}}^n:={\displaystyle \frac{1}{4}}({\mathbf{U}}^{n+1}+2{\mathbf{U}}^n+{\mathbf{U}}^{n-1})$, and the L1 Caputo memory is evaluated explicitly at time $t_n$:

$${\mathbf{H}}_u^n=\sum _{j=0}^{n-1}{w}_j^{\left(\alpha \right)}({\mathbf{U}}^{n-j}-{\mathbf{U}}^{n-j-1}),w_k^{\left(\gamma \right)}={(k+1)}^{1-\gamma }-k^{1-\gamma }.$$

(8)

Here $a_0^{\left(\alpha \right)}=\Delta t^{-\alpha }/\Gamma (2-\alpha )$ and ${\mathbf{D}}^2$ is the second-order Chebyshev spectral differentiation matrix. Thus ${\sigma }_i{a}_0^{\left(\gamma \right)}{\mathbf{H}}_u^n$ approximates with $\mathcal{O}(\Delta t^{2-\gamma })$ accuracy ([23]). The vector ${\mathbf{U}}^n={[u(x_0,t_n),\dots ,u(x_N,t_n)]}^T$ represents the spatial solution vector. The scheme is initialized by a Taylor expansion preserving $\mathcal{O}(\Delta t^2)$ accuracy:

$${\mathbf{U}}^1={\varphi }_1+\Delta t{\psi }_1+{\displaystyle \frac{\Delta t^2}{2}}\left(c_1^2{\mathbf{D}}^2{\varphi }_1-m_1{\varphi }_1-{\mathbf{N}}_u({\varphi }_1,{\varphi }_2)+{\mathbf{f}}_1^0\right),$$

(9)

and analogously for ${\mathbf{V}}^1$. Note that no memory term appears in the initialisation since ^C$D_t^{\gamma }u(·,0)=0$ for the smooth initial data used here (the Caputo derivative of a constant is zero). For $n\ge 1$ the nonlinear terms $\mathbf{N}$ are evaluated at the linearly extrapolated values (detailed in Section 2.4); in the initialisation step, they are evaluated at the exact initial data.

2.1. Spatial Chebyshev Collocation

The continuous physical domain $[0,L]$ is mapped to the standard Chebyshev domain $[-1,1]$ and discretized using the shifted Chebyshev-Gauss-Lobatto (CGL) collocation points:

$$x_j={\displaystyle \frac{L}{2}}\left(1-cos\left({\displaystyle \frac{j\pi }{N}}\right)\right),j=0,1,\dots ,N.$$

(10)

The continuous spatial derivatives are approximated exactly up to polynomial degree N by the spectral differentiation matrix D. For a discrete function vector $U={[u_0,u_1,\dots ,u_N]}^T$, the first derivative vector is given analytically by $U^{\prime }=DU$, where the matrix entries $D_{j,k}$ are defined on the standard Chebyshev nodes ${\xi }_j=cos(j\pi /N)$ as:

$$D_{j,k}=\left\{\begin{array}{cc}{\displaystyle \frac{c_j}{c_k}}{\displaystyle \frac{{(-1)}^{j+k}}{{\xi }_j-{\xi }_k}},\hfill & j\ne k,\hfill \\ -{\displaystyle \frac{{\xi }_j}{2(1-{\xi }_j^2)}},\hfill & 1\le j=k\le N-1,\hfill \\ {\displaystyle \frac{2N^2+1}{6}},\hfill & j=k=0,\hfill \\ -{\displaystyle \frac{2N^2+1}{6}},\hfill & j=k=N,\hfill \end{array}\right.$$

(11)

with the boundary weights fixed at $c_0=c_N=2$, and $c_j=1$ for internal nodes $1\le j\le N-1$. To account for the physical domain length L, the second-order spectral differentiation matrix on $[0,L]$ is assembled as ${\mathbf{D}}_{[0,L]}^2={(2/L)}^2{D}^2$, where D is the first-order Chebyshev matrix on $[-1,1]$[17]. This spectral formulation provides exponential spatial convergence $\mathcal{O}\left(e^{-cN}\right)$ for smooth solutions [16].

2.2. Dual Fast Memory Acceleration

Evaluating the Caputo memory convolution in (6) at time $t_n$ requires summing over all previous time steps, leading to $\mathcal{O}\left(n\right)$ operations per step and $\mathcal{O}\left(N_t^2\right)$ total temporal complexity. To reduce this cost, we employ a Sum-of-Exponentials (SOE) approximation following Jiang et al. [11]. The key observation is that the discrete $L1$ weights admit the integral representation

$$w_j={(j+1)}^{1-\gamma }-j^{1-\gamma }={\displaystyle \frac{1-\gamma }{\Gamma \left(\gamma \right)}}{\int }_0^{\infty }{s}^{\gamma -2}(1-e^{-s})e^{-sj}ds,$$

(12)

where the factor $(1-e^{-s})/s$ corrects the continuous kernel $j^{-\gamma }$ to the exact discrete difference. Applying a double-exponential (DE) quadrature ([24]) to this integral, the tail weights ($j\ge p$) are approximated as:

$$w_j\approx \sum _{l=1}^{N_{\exp }}{c}_l{q}_l^j,q_l=e^{-s_l},c_l={\displaystyle \frac{(1-\gamma )}{\Gamma \left(\gamma \right)}}{h}_{\mathrm{DE}}{s}_l^{\gamma -1}(1-e^{-s_l}){\phi }^{\prime }\left(x_l\right),$$

(13)

where $s_l=exp\left(\phi \left(x_l\right)\right)$ are the DE-transformed nodes with $\phi \left(x\right)={\displaystyle \frac{\pi }{2}}sinh\left(x\right)$, ${\phi }^{\prime }\left(x_l\right)={\displaystyle \frac{\pi }{2}}cosh\left(x_l\right)$, and $h_{\mathrm{DE}}$ is the DE step size. The exponential structure $q_l^j$ enables a one-step recursion for the history accumulators:

$$Z_l^n=q_l{Z}_l^{n-1}+({\mathbf{U}}^{n-p}-{\mathbf{U}}^{n-p-1}),Z_l^0=\mathbf{0}.$$

(14)

The full memory sum is then assembled from the exact near part (last p increments) and the SOE-accelerated tail:

$${\mathbf{H}}^n=\sum _{j=0}^{p-1}{w}_j({\mathbf{U}}^{n-j}-{\mathbf{U}}^{n-j-1})+\sum _{l=1}^{N_{\exp }}{c}_l{Z}_l^n.$$

(15)

Since updating $Z_l^n$ requires only the previous accumulator $Z_l^{n-1}$ and the local increment $({\mathbf{U}}^{n-p},{\mathbf{U}}^{n-p-1})$, earlier solution vectors can be discarded from memory. This reduces the per-step cost from $\mathcal{O}\left(n\right)$ to $\mathcal{O}(p+N_{\exp })$, and the total temporal complexity from $\mathcal{O}\left(N_t^2\right)$ to $\mathcal{O}\left(N_t(p+N_{\exp })\right)$.

2.3. Linearization of Nonlinear Terms via Explicit Extrapolation

The nonlinear terms in (6) — the cubic self-interactions $k_1{u}^3$, $k_2{v}^3$ and the inter-component coupling $\eta u{v}^2$, $\eta v{u}^2$ — depend on the unknown solution at time $t_{n+1}$. A fully implicit treatment would require solving a nonlinear algebraic system at every time step via Newton–Raphson iterations. To avoid this cost while retaining second-order temporal accuracy, we adopt an implicit–explicit (IMEX) strategy: the linear spatial operator (inertia and Laplacian via Newmark averaging) is treated implicitly, while the memory sums, nonlinear terms, and mass terms are evaluated explicitly at known data.

At each time step, the nonlinear arguments are approximated by the linearly extrapolated values:

$$U^{\ast }=2{\mathbf{U}}^n-{\mathbf{U}}^{n-1},V^{\ast }=2{\mathbf{V}}^n-{\mathbf{V}}^{n-1}.$$

(16)

By Taylor expansion, $U^{\ast }={\mathbf{U}}^{n+1}-\Delta t^2{u}_{tt}\left(t_n\right)+\mathcal{O}(\Delta t^3)$, so the predictors approximate ${\mathbf{U}}^{n+1}$ with $\mathcal{O}(\Delta t^2)$ accuracy. The nonlinear forcing vectors in (7) are then evaluated entirely from known data:

$${\mathbf{N}}_u({\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast })=k_1{\left(U^{\ast }\right)}^3+\eta U^{\ast }{\left(V^{\ast }\right)}^2,{\mathbf{N}}_v({\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast })=k_2{\left(V^{\ast }\right)}^3+\eta V^{\ast }{\left(U^{\ast }\right)}^2,$$

(17)

where all operations are understood componentwise on the collocation grid.

This explicit treatment offers two critical algorithmic advantages:

• Constant coefficient matrix. Since the nonlinear terms ${\mathbf{N}}_u,{\mathbf{N}}_v$ are evaluated at known data and the memory sum ${\mathbf{H}}_u^n$ consists entirely of previously computed solution vectors, these terms are moved to the right-hand side. The left-hand side linear operator

  $${\mathbf{A}}_u={\displaystyle \frac{1}{\Delta t^2}}\mathbf{I}-{\displaystyle \frac{c_1^2}{4}}{\mathbf{D}}^2$$

  (18)

  remains identical at every time step (and analogously ${\mathbf{A}}_v$ for the v-component). A single LU factorization of ${\mathbf{A}}_u$ and ${\mathbf{A}}_v$ computed once prior to the time loop suffices for the entire simulation. Each time step reduces to two back-substitutions.
• No CFL restriction. The stiff eigenvalues that typically necessitate a CFL condition originate from the spectral Laplacian ${\mathbf{D}}^2$, whose eigenvalues grow as $\mathcal{O}\left(N^4\right)$. Because ${\mathbf{D}}^2$ is treated implicitly via the Newmark averaging (${\beta }_{\mathrm{NM}}=1/4$), these large eigenvalues are unconditionally damped by the matrix ${\mathbf{A}}_u$. The explicitly evaluated nonlinear terms ${\mathbf{N}}_u({\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast })$ and the memory sum ${\mathbf{H}}_u^n$ are bounded functions of previously computed data and do not introduce additional spectral radius constraints on $\Delta t$. Hence, for fixed spatial resolution N and bounded solution norms, the scheme is stable for any $\Delta t>0$ (proved rigorously in Theorem 1).

3. Mathematical Analysis

We introduce the following discrete inner product and norm on the Chebyshev collocation grid ${\left\{x_j\right\}}_{j=0}^N$ equipped with the Clenshaw–Curtis quadrature weights $\left\{w_j\right\}$:

$${\langle \mathbf{U},\mathbf{V}\rangle }_N=\sum _{j=0}^N{w}_j{U}_j{V}_j,{\parallel \mathbf{U}\parallel }_N=\sqrt{{\langle \mathbf{U},\mathbf{U}\rangle }_N}.$$

(19)

Since the boundary conditions are inhomogeneous, we employ the standard boundary lifting decomposition. Let ${\mathbf{G}}_u^n$ denote the discrete boundary lift satisfying the prescribed Dirichlet values $g_1\left(t_n\right),h_1\left(t_n\right)$ at the endpoints, and define the homogeneous residual ${\tilde{\mathbf{U}}}^n={\mathbf{U}}^n-{\mathbf{G}}_u^n$. The lifted variable ${\tilde{\mathbf{U}}}^n$ vanishes at $\partial \Omega$. The Chebyshev second-derivative spectral matrix ${\mathbf{D}}^2$, when restricted to functions satisfying homogeneous Dirichlet conditions, possesses the fundamental semi-definiteness property:

$$-{\langle {\mathbf{D}}^2\tilde{\mathbf{U}},\tilde{\mathbf{U}}\rangle }_N={|\tilde{\mathbf{U}}|}_{H_N^1}^2\ge 0,$$

(20)

where ${|·|}_{H_N^1}$ is the discrete $H^1$ seminorm. We shall use property (20) as the core mechanism for deriving data-dependent bounds. The boundary lift contributions (derivatives of ${\mathbf{G}}_u^n,{\mathbf{G}}_v^n$) enter the right-hand side as additional known source terms.

3.1. Theorem 1: A Priori Stability Estimate (Data-Dependent Bound)

Theorem 1.

Let ${\{{\mathbf{U}}^n,{\mathbf{V}}^n\}}_{n=0}^{N_t}$ be the solution of the fully discrete coupled scheme (7), in which the cubic nonlinearities ${\mathbf{N}}_u,{\mathbf{N}}_v$ are evaluated at the linear extrapolants ${\mathbf{U}}^{\ast }=2{\mathbf{U}}^n-{\mathbf{U}}^{n-1}$ and the L1 memory sums ${\mathbf{H}}_u^n,{\mathbf{H}}_v^n$ are evaluated explicitly at $t_n$ from the exact discrete weights (without SOE compression). Then there exists a constant $C>0$, depending on the physical parameters $c_1,c_2,m_1,m_2,k_1,k_2,\eta$, the final time T, and the number of collocation nodes N (through the inverse inequality constant $C_N=\mathcal{O}\left(\sqrt{N}\right)$ used to bound the pointwise nonlinear terms), such that for all $1\le n\le N_t$:

$$\begin{array}{c}\parallel {\mathbf{U}}^n{\parallel }_N^2+\parallel {\mathbf{V}}^n{\parallel }_N^2\le C(\parallel {\varphi }_1{\parallel }_N^2+\parallel {\varphi }_2{\parallel }_N^2+\parallel {\psi }_1{\parallel }_N^2+{\parallel {\psi }_2\parallel }_N^2\hfill \\ \hfill +\underset{0\le k\le n}{max}{\parallel {\mathbf{G}}^k\parallel }_N^2+\Delta t\sum _{k=0}^n\left(\parallel {\mathbf{f}}_1^k{\parallel }_N^2+{\parallel {\mathbf{f}}_2^k\parallel }_N^2\right)),\end{array}$$

(21)

where ${\varphi }_i,{\psi }_i$ are the discrete initial data vectors, ${\mathbf{f}}_i^k$ are the discrete source vectors, and ${\mathbf{G}}^k$ collects the boundary lift contributions from the prescribed Dirichlet data $g_i\left(t_k\right),h_i\left(t_k\right)$.

Proof. 

We derive the estimates for both components $\mathbf{U}$ and $\mathbf{V}$ simultaneously, since the inter-component coupling terms $\eta U^{\ast }{\left(V^{\ast }\right)}^2$ and $\eta V^{\ast }{\left(U^{\ast }\right)}^2$ prevent treating each equation independently. First, we test the first equation in (7) with the symmetric discrete velocity ${\widehat{\mathbf{U}}}^n:={\displaystyle \frac{{\mathbf{U}}^{n+1}-{\mathbf{U}}^{n-1}}{2\Delta t}}={\displaystyle \frac{{\mathbf{d}}_u^n+{\mathbf{d}}_u^{n-1}}{2}}$ and the second equation with ${\widehat{\mathbf{V}}}^n:={\displaystyle \frac{{\mathbf{V}}^{n+1}-{\mathbf{V}}^{n-1}}{2\Delta t}}$, obtaining the pair:

$$\begin{array}{cc}\hfill & {〈{\displaystyle \frac{{\mathbf{U}}^{n+1}-2{\mathbf{U}}^n+{\mathbf{U}}^{n-1}}{\Delta t^2}},{\widehat{\mathbf{U}}}^n〉}_N-c_1^2{\langle {\mathbf{D}}^2{\overline{\mathbf{U}}}^n,{\widehat{\mathbf{U}}}^n\rangle }_N+a_0^{\left(\alpha \right)}{\sigma }_1{\langle {\mathbf{H}}_u^n,{\widehat{\mathbf{U}}}^n\rangle }_N\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +m_1{\langle {\mathbf{U}}^n,{\widehat{\mathbf{U}}}^n\rangle }_N+{\langle {\mathbf{N}}_u,{\widehat{\mathbf{U}}}^n\rangle }_N={\langle {\mathbf{f}}_1^n,{\widehat{\mathbf{U}}}^n\rangle }_N,\hfill \\ \hfill & {〈{\displaystyle \frac{{\mathbf{V}}^{n+1}-2{\mathbf{V}}^n+{\mathbf{V}}^{n-1}}{\Delta t^2}},{\widehat{\mathbf{V}}}^n〉}_N-c_2^2{\langle {\mathbf{D}}^2{\overline{\mathbf{V}}}^n,{\widehat{\mathbf{V}}}^n\rangle }_N+a_0^{\left(\beta \right)}{\sigma }_2{\langle {\mathbf{H}}_v^n,{\widehat{\mathbf{V}}}^n\rangle }_N\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & +m_2{\langle {\mathbf{V}}^n,{\widehat{\mathbf{V}}}^n\rangle }_N+{\langle {\mathbf{N}}_v,{\widehat{\mathbf{V}}}^n\rangle }_N={\langle {\mathbf{f}}_2^n,{\widehat{\mathbf{V}}}^n\rangle }_N.\hfill \end{array}$$

(23)

Subsequently, applying the identity $({\mathbf{U}}^{n+1}-2{\mathbf{U}}^n+{\mathbf{U}}^{n-1})=\Delta t({\mathbf{d}}_u^n-{\mathbf{d}}_u^{n-1})$, the inertial term yields a perfect discrete derivative:

$${〈{\displaystyle \frac{{\mathbf{U}}^{n+1}-2{\mathbf{U}}^n+{\mathbf{U}}^{n-1}}{\Delta t^2}},{\displaystyle \frac{{\mathbf{U}}^{n+1}-{\mathbf{U}}^{n-1}}{2\Delta t}}〉}_N={\displaystyle \frac{\parallel {\mathbf{d}}_u^n{\parallel }_N^2-{\parallel {\mathbf{d}}_u^{n-1}\parallel }_N^2}{2\Delta t}}.$$

(24)

For the Laplacian, using ${\overline{\mathbf{U}}}^n={\displaystyle \frac{1}{4}}({\mathbf{U}}^{n+1}+{\mathbf{U}}^n)+{\displaystyle \frac{1}{4}}({\mathbf{U}}^n+{\mathbf{U}}^{n-1})$, decomposing ${\mathbf{U}}^n={\tilde{\mathbf{U}}}^n+{\mathbf{G}}_u^n$, and incorporating the semi-definiteness (20), we extract the telescoping spatial strain energy:

$$-c_1^2{\langle {\mathbf{D}}^2{\overline{\tilde{\mathbf{U}}}}^n,{\widehat{\tilde{\mathbf{U}}}}^n\rangle }_N={\displaystyle \frac{1}{2\Delta t}}\left({\displaystyle \frac{c_1^2}{4}}|{\tilde{\mathbf{U}}}^{n+1}+{\tilde{\mathbf{U}}}^n{|}_{H_N^1}^2-{\displaystyle \frac{c_1^2}{4}}{|{\tilde{\mathbf{U}}}^n+{\tilde{\mathbf{U}}}^{n-1}|}_{H_N^1}^2\right).$$

(25)

The boundary lifting cross-terms are bounded by the Cauchy–Schwarz inequality and shifted to the right-hand side. Analogous relations hold for the $\mathbf{V}$-component.

To bound the explicit memory and mass terms, we invoke the Cauchy–Schwarz and Young inequalities, observing that $\parallel {\widehat{\mathbf{U}}}^n{\parallel }_N^2\le {\displaystyle \frac{1}{2}}(\parallel {\mathbf{d}}_u^n{\parallel }_N^2+\parallel {\mathbf{d}}_u^{n-1}{\parallel }_N^2)$:

$$\begin{array}{cc}\hfill |a_0^{\left(\alpha \right)}{\sigma }_1{\langle {\mathbf{H}}_u^n,{\widehat{\mathbf{U}}}^n\rangle }_N|& \le {\displaystyle \frac{1}{2}}{\left(a_0^{\left(\alpha \right)}{\sigma }_1\right)}^2\parallel {\mathbf{H}}_u^n{\parallel }_N^2+{\displaystyle \frac{1}{2}}{\parallel {\widehat{\mathbf{U}}}^n\parallel }_N^2,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill |m_1{\langle {\mathbf{U}}^n,{\widehat{\mathbf{U}}}^n\rangle }_N|& \le {\displaystyle \frac{m_1^2}{2}}\parallel {\mathbf{U}}^n{\parallel }_N^2+{\displaystyle \frac{1}{2}}{\parallel {\widehat{\mathbf{U}}}^n\parallel }_N^2.\hfill \end{array}$$

(27)

To bound the memory, we use the summation property ${\sum }_{j=0}^{n-1}{w}_j^{\left(\alpha \right)}=n^{1-\alpha }$ and ${\mathbf{U}}^{n-j}-{\mathbf{U}}^{n-j-1}=\Delta t{\mathbf{d}}_u^{n-j-1}$:

$$\parallel {\mathbf{H}}_u^n{\parallel }_N\le \sum _{j=0}^{n-1}{w}_j^{\left(\alpha \right)}\Delta t\parallel {\mathbf{d}}_u^{n-j-1}{\parallel }_N\le n^{1-\alpha }\Delta t\underset{0\le k\le n}{max}{\parallel {\mathbf{d}}_u^k\parallel }_N.$$

(28)

Since $a_0^{\left(\alpha \right)}=\Delta t^{-\alpha }/\Gamma (2-\alpha )$ and $(n\Delta t)\le T$, the scaled memory is bounded uniformly and explicitly:

$$\begin{array}{cc}\hfill {\left(a_0^{\left(\alpha \right)}\right)}^2{\parallel {\mathbf{H}}_u^n\parallel }_N^2& \le {\left({\displaystyle \frac{\Delta t^{-\alpha }}{\Gamma (2-\alpha )}}\right)}^2{\left(n^{1-\alpha }\Delta t\right)}^2\underset{0\le k\le n}{max}{\parallel {\mathbf{d}}_u^k\parallel }_N^2\hfill \\ \hfill & \le {\displaystyle \frac{T^{2-2\alpha }}{\Gamma {(2-\alpha )}^2}}\underset{0\le k\le n}{max}{\parallel {\mathbf{d}}_u^k\parallel }_N^2.\hfill \end{array}$$

(29)

Control over the nonlinear coupling structure is established by evaluating the terms ${\mathbf{N}}_u,{\mathbf{N}}_v$ at the explicit extrapolants ${\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast }$. By Cauchy–Schwarz:

$$|{\langle {\mathbf{N}}_u,{\widehat{\mathbf{U}}}^n\rangle }_N|\le {\displaystyle \frac{1}{2}}\parallel {\mathbf{N}}_u{\parallel }_N^2+{\displaystyle \frac{1}{2}}\parallel {\widehat{\mathbf{U}}}^n{\parallel }_N^2,|{\langle {\mathbf{N}}_v,{\widehat{\mathbf{V}}}^n\rangle }_N|\le {\displaystyle \frac{1}{2}}\parallel {\mathbf{N}}_v{\parallel }_N^2+{\displaystyle \frac{1}{2}}{\parallel {\widehat{\mathbf{V}}}^n\parallel }_N^2.$$

(30)

Bounding $\parallel {\mathbf{N}}_u{\parallel }_N$ requires controlling both components via the inter-component coupling. Using the Chebyshev inverse inequality ${\parallel W\parallel }_{L^{\infty }}\le C_N{\parallel W\parallel }_N$ (with $C_N=\mathcal{O}\left(\sqrt{N}\right)$) twice:

$$\parallel \eta U^{\ast }{\left(V^{\ast }\right)}^2{\parallel }_N\le \eta \parallel U^{\ast }{\parallel }_{L^{\infty }}\parallel V^{\ast }{\parallel }_{L^{\infty }}\parallel V^{\ast }{\parallel }_N\le \eta C_N^2\parallel {\mathbf{U}}^{\ast }{\parallel }_N{\parallel {\mathbf{V}}^{\ast }\parallel }_N^2.$$

(31)

The self-interaction $\parallel k_1{\left(U^{\ast }\right)}^3{\parallel }_N\le k_1{C}_N^2{\parallel {\mathbf{U}}^{\ast }\parallel }_N^3$ is similarly bounded. Summing both components:

$$\parallel {\mathbf{N}}_u{\parallel }_N^2+{\parallel {\mathbf{N}}_v\parallel }_N^2\le C(k_1,k_2,\eta ,N){\left(\parallel {\mathbf{U}}^n{\parallel }_N+\parallel {\mathbf{U}}^{n-1}{\parallel }_N+\parallel {\mathbf{V}}^n{\parallel }_N+{\parallel {\mathbf{V}}^{n-1}\parallel }_N\right)}^6.$$

(32)

The localized source terms are similarly bounded by standard inequalities:

$$|{\langle {\mathbf{f}}_1^n,{\widehat{\mathbf{U}}}^n\rangle }_N|\le {\displaystyle \frac{1}{2}}\parallel {\mathbf{f}}_1^n{\parallel }_N^2+{\displaystyle \frac{1}{2}}\parallel {\widehat{\mathbf{U}}}^n{\parallel }_N^2,|{\langle {\mathbf{f}}_2^n,{\widehat{\mathbf{V}}}^n\rangle }_N|\le {\displaystyle \frac{1}{2}}\parallel {\mathbf{f}}_2^n{\parallel }_N^2+{\displaystyle \frac{1}{2}}{\parallel {\widehat{\mathbf{V}}}^n\parallel }_N^2.$$

(33)

Defining the effective discrete kinetic and strain energy for the system:

$$\begin{array}{cc}\hfill {\mathcal{E}}^n& :=\parallel {\mathbf{d}}_u^n{\parallel }_N^2+\parallel {\mathbf{d}}_v^n{\parallel }_N^2+{\displaystyle \frac{c_1^2}{4}}{|{\tilde{\mathbf{U}}}^{n+1}+{\tilde{\mathbf{U}}}^n|}_{H_N^1}^2\hfill \\ \hfill & +{\displaystyle \frac{c_2^2}{4}}|{\tilde{\mathbf{V}}}^{n+1}+{\tilde{\mathbf{V}}}^n{|}_{H_N^1}^2+\parallel {\mathbf{U}}^{n+1}{\parallel }_N^2+{\parallel {\mathbf{V}}^{n+1}\parallel }_N^2.\hfill \end{array}$$

(34)

Multiplying the inner products (22)–(23) by $2\Delta t$, substituting Steps 2–5, absorbing the finite memory limit (29), and leveraging the fact that $\parallel {\mathbf{U}}^{n+1}{\parallel }_N^2$ is bounded by the history $max{\mathcal{E}}^k$ via ${\mathbf{U}}^{n+1}={\mathbf{U}}^0+\Delta t\sum {\mathbf{d}}_u^k$, we obtain a telescoping recursive bound:

$${\mathcal{E}}^n-{\mathcal{E}}^{n-1}\le C_1\Delta t\underset{0\le k\le n}{max}{\mathcal{E}}^k+C_2\Delta t\left(\parallel {\mathbf{f}}_1^n{\parallel }_N^2+\parallel {\mathbf{f}}_2^n{\parallel }_N^2+{\parallel {\mathbf{G}}^n\parallel }_N^2\right),$$

(35)

where $C_1$ depends on the physical constants, T, and N. Taking the maximum over n on the left-hand side transforms this into ${\mathcal{E}}_{max}^n\le {\mathcal{E}}_{max}^{n-1}+C_1\Delta t{\mathcal{E}}_{max}^n+\dots$, to which the classical integer-order discrete Gronwall lemma [25] readily applies. The stability bound (21) directly follows on any finite interval $[0,T]$ where the polynomial bound (32) is dynamically linearisable. □

Remark 1.

Estimate (21) is a data-dependent bound, meaning the discrete solution at any time level is controlled by the prescribed initial conditions, boundary values, and source magnitudes. The stability constant C inherently grows with N due to the use of the Chebyshev inverse inequality ${\parallel U\parallel }_{L^{\infty }}\le C_N{\parallel U\parallel }_N$ (with $C_N=\mathcal{O}\left(\sqrt{N}\right)$) required to bound the pointwise cubic nonlinearities. Consequently, for a fixed spatial resolution N and bounded solution norms, the scheme remains unconditionally stable for any $\Delta t>0$.

When the SOE fast evaluation is explicitly activated, the exact discrete memory sum ${\mathbf{H}}_u^n$ is replaced by its exponential approximation ${\widehat{\mathbf{H}}}_u^n$, introducing a truncation perturbation ${\delta }_{\mathrm{SOE}}^n={\widehat{\mathbf{H}}}_u^n-{\mathbf{H}}_u^n$. Because the memory integral is evaluated explicitly and shifted to the right-hand side of the linear system, this ${\delta }_{\mathrm{SOE}}^n$ effectively acts as an additive discrete source. The overarching stability bound (21) maintains its validity under the augmented source profile ${\mathbf{f}}_i^k+{\delta }_{\mathrm{SOE}}^k$. In practical implementations, the magnitude $\parallel {\delta }_{\mathrm{SOE}}^n{\parallel }_N$ is strictly governed by the Gauss–Legendre and Double Exponential quadrature fidelity, which remains negligible at double precision for $N_{\exp }\ge 40$ (as verified numerically in Section 4).

3.2. Theorem 2: Global Convergence Order

Theorem 2.

Let ${\mathbf{e}}_u^n={\mathbf{U}}^n-{\mathbf{u}}^n$ and ${\mathbf{e}}_v^n={\mathbf{V}}^n-{\mathbf{v}}^n$ denote the global discretization errors at the collocation nodes, where ${\mathbf{u}}^n,{\mathbf{v}}^n$ are the exact solution vectors. Assume that the exact solutions $u,v\in C^4([0,T];H^s(\Omega ))$ for some $s\ge 1$, and that the a priori bound of Theorem 1 ensures that both the numerical and exact solutions remain in a bounded set where the nonlinear map $(u,v)\mapsto (k_1{u}^3+\eta u{v}^2,k_2{v}^3+\eta v{u}^2)$ is locally Lipschitz continuous. Then there exists a constant $C>0$, independent of $\Delta t$ but dependent on N through the Lipschitz constant and the inverse inequality (see Theorem 1), such that:

$$\parallel {\mathbf{e}}_u^n{\parallel }_N^2+{\parallel {\mathbf{e}}_v^n\parallel }_N^2\le C\left(\Delta t^{2(2-\alpha )}+\Delta t^{2(2-\beta )}+N^{-2s}+{\epsilon }_{\mathrm{SOE}}^2\right)\mathrm{for}\mathrm{all}1\le n\le N_t,$$

(36)

where s is the Sobolev regularity exponent of the exact solution, and ${\epsilon }_{\mathrm{SOE}}$ denotes the maximum SOE quadrature error in approximating the memory tail (set ${\epsilon }_{\mathrm{SOE}}=0$ when no SOE compression is used).

Proof. 

Subtracting the exact discrete equations from the numerical scheme (7) produces the error system:

$$\left\{\begin{array}{c}{\displaystyle \frac{{\mathbf{e}}_u^{n+1}-2{\mathbf{e}}_u^n+{\mathbf{e}}_u^{n-1}}{\Delta t^2}-c_1^2{\mathbf{D}}^2{\overline{\mathbf{e}}}_u^n+a_0^{\left(\alpha \right)}{\sigma }_1{\mathbf{H}}_{e_u}^n+m_1{\mathbf{e}}_u^n}\hfill \\ +{\mathbf{N}}_u({\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast })-{\mathbf{N}}_u({\mathbf{u}}^{\ast },{\mathbf{v}}^{\ast })={\tau }_u^n+{\delta }_{\mathrm{SOE},u}^n,\hfill \\ {\displaystyle \frac{{\mathbf{e}}_v^{n+1}-2{\mathbf{e}}_v^n+{\mathbf{e}}_v^{n-1}}{\Delta t^2}-c_2^2{\mathbf{D}}^2{\overline{\mathbf{e}}}_v^n+a_0^{\left(\beta \right)}{\sigma }_2{\mathbf{H}}_{e_v}^n+m_2{\mathbf{e}}_v^n}\hfill \\ +{\mathbf{N}}_v({\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast })-{\mathbf{N}}_v({\mathbf{u}}^{\ast },{\mathbf{v}}^{\ast })={\tau }_v^n+{\delta }_{\mathrm{SOE},v}^n,\hfill \end{array}\right.$$

(37)

where ${\mathbf{H}}_{e_u}^n={\sum }_{j=0}^{n-1}{w}_j^{\left(\alpha \right)}({\mathbf{e}}_u^{n-j}-{\mathbf{e}}_u^{n-j-1})$ is the L1 memory of the error (evaluated explicitly at $t_n$), ${\tau }_u^n,{\tau }_v^n$ are the local truncation error vectors, and ${\delta }_{\mathrm{SOE},u}^n,{\delta }_{\mathrm{SOE},v}^n$ are the SOE compression errors (zero when full L1 weights are used).

The system’s truncation errors stem from three primary analytical components. First, the temporal approximation of the Caputo fractional derivative using the $L1$ scheme introduces a local error of $\mathcal{O}(\Delta t^{2-\alpha })$ for the u-component [23], alongside a similarly derived $\mathcal{O}(\Delta t^{2-\beta })$ error for the v-component. Second, the implicit spatial evaluation structured through Newmark averaging provides an $\mathcal{O}(\Delta t^2)$ geometric convergence (since symmetric Taylor expansion yields ${\displaystyle \frac{1}{4}}{u}^{n+1}+{\displaystyle \frac{1}{2}}{u}^n+{\displaystyle \frac{1}{4}}{u}^{n-1}=u^n+{\displaystyle \frac{\Delta t^2}{4}}{u}_{tt}\left(t_n\right)+\mathcal{O}(\Delta t^4)$), while the nonlinear explicit extrapolation contributes an identically ordered $\mathcal{O}(\Delta t^2)$ error. Finally, spatial collocation using Chebyshev clustering bounds the interpolation discrepancy by $\mathcal{O}\left(N^{-s}\right)$ for solutions maintaining Sobolev regularity $H^s(\Omega )$ [16]. Aggregating these sources leaves the temporal fractional derivative as the dominant source of error:

$$\parallel {\tau }_u^n{\parallel }_N\le C_{\tau }\left(\Delta t^{2-\alpha }+N^{-s}\right),{\parallel {\tau }_v^n\parallel }_N\le C_{\tau }\left(\Delta t^{2-\beta }+N^{-s}\right).$$

(38)

Furthermore, due to the local Lipschitz continuity of the cubic polynomial nonlinearities on bounded domains (dynamically guaranteed by the a priori stability established in Theorem 1), the nonlinear residual satisfies:

$$\parallel {\mathbf{N}}_u({\mathbf{U}}^{\ast },{\mathbf{V}}^{\ast })-{\mathbf{N}}_u({\mathbf{u}}^{\ast },{\mathbf{v}}^{\ast }){\parallel }_N\le L\left(\parallel {\mathbf{e}}_u^n{\parallel }_N+\parallel {\mathbf{e}}_u^{n-1}{\parallel }_N+\parallel {\mathbf{e}}_v^n{\parallel }_N+{\parallel {\mathbf{e}}_v^{n-1}\parallel }_N\right),$$

(39)

where the Lipschitz constant L absorbs $k_1,k_2,\eta$ and the a priori bounded solution magnitudes.

Testing the error equations (37) with the discrete error velocities ${\widehat{\mathbf{e}}}_u^n:={\displaystyle \frac{{\mathbf{e}}_u^{n+1}-{\mathbf{e}}_u^{n-1}}{2\Delta t}}$ and ${\widehat{\mathbf{e}}}_v^n$ exactly mirrors the energy technique used in Theorem 1. The explicit memory ${\mathbf{H}}_{e_u}^n$ is uniformly bounded by the discrete velocity history, the Lipschitz-bounded nonlinear residual collapses into the right-hand side, and the geometric truncation errors formulate the driving source mechanism. Defining the total error energy ${\mathcal{E}}_{err}^n$ identically to the physical kinetic-strain energy, we recover the telescoping inequality:

$$\begin{array}{cc}\hfill {\mathcal{E}}_{err}^n-{\mathcal{E}}_{err}^{n-1}& \le C_1\Delta t\underset{0\le k\le n}{max}{\mathcal{E}}_{err}^k\hfill \\ \hfill & +C_2\Delta t\left(\parallel {\tau }_u^n{\parallel }_N^2+\parallel {\tau }_v^n{\parallel }_N^2+\parallel {\delta }_{\mathrm{SOE},u}^n{\parallel }_N^2+{\parallel {\delta }_{\mathrm{SOE},v}^n\parallel }_N^2\right).\hfill \end{array}$$

(40)

Since the initial errors vanish (${\mathbf{e}}_u^0={\mathbf{e}}_v^0=\mathbf{0}$) and taking the maximum over n on the left-hand side, the standard discrete Gronwall inequality ([25]) yields:

$$\begin{array}{cc}\hfill {\mathcal{E}}_{err}^n& \le e^{C_{1}T}\Delta t\sum _{k=1}^n\left(\parallel {\tau }_u^k{\parallel }_N^2+\parallel {\tau }_v^k{\parallel }_N^2+\parallel {\delta }_{\mathrm{SOE},u}^k{\parallel }_N^2+{\parallel {\delta }_{\mathrm{SOE},v}^k\parallel }_N^2\right)\hfill \\ \hfill & \le C\left(\Delta t^{2(2-\alpha )}+\Delta t^{2(2-\beta )}+N^{-2s}+{\epsilon }_{\mathrm{SOE}}^2\right),\hfill \end{array}$$

(41)

where ${\epsilon }_{\mathrm{SOE}}={max}_{1\le k\le n}max(\parallel {\delta }_{\mathrm{SOE},u}^k{\parallel }_N,\parallel {\delta }_{\mathrm{SOE},v}^k{\parallel }_N)$. This establishes the convergence estimate (36). □

Remark 2.

The appearance of the fractional exponent $2(2-\alpha )$ in the established bound (36) analytically arises from targeting the squared spectral $L^2$-norm $\parallel {\mathbf{e}}_u^n{\parallel }_N^2$. Upon taking the square root, the global convergence rate gracefully scales as $\mathcal{O}(\Delta t^{2-\alpha }+N^{-s})$, fully consistent with standard $L1$ formulation characteristics. The global numerical error trajectory is strictly dictated by the slower fractional memory relaxation scale, equating to a temporal asymptote of $\mathcal{O}(\Delta t^{min(2-\alpha ,2-\beta )})$ which is extensively validated in the numerical experiments.

4. Numerical Results

4.1. Unified Experimental Setup

All numerical experiments are conducted on the domain $\Omega =(0,1)$ with final time $T=1$ using the Chebyshev spectral collocation scheme described in Section 2.1–2.4. To verify the theoretical convergence estimates of Theorem 2, we employ the method of manufactured solutions. The exact solution pair is prescribed as:

$$u(x,t)=t^{3}sin\left(\pi x\right),v(x,t)=t^{3}sin\left(2\pi x\right).$$

(42)

These functions satisfy homogeneous Dirichlet boundary conditions $u(0,t)=u(1,t)=0$, $v(0,t)=v(1,t)=0$ (i.e., $g_1=h_1=g_2=h_2=0$), and the initial conditions $u(x,0)=v(x,0)=0$, ${\partial }_{t}u(x,0)={\partial }_{t}v(x,0)=0$. The source terms $f_1(x,t)$ and $f_2(x,t)$ are computed analytically by substituting (42) into the governing system (6), using the Caputo derivative identity ^C$D_t^{\gamma }{t}^3={\displaystyle \frac{6}{\Gamma (4-\gamma )}}{t}^{3-\gamma }$[1]. The physical parameters are fixed at $c_1=c_2=1$, $m_1=m_2=1$, $k_1=k_2=1$, ${\sigma }_1={\sigma }_2=1$, and $\eta =1.0$ unless stated otherwise.

Remark 3.

The underlying manufactured profiles $t^{3}sin\left(\pi x\right)$ were strategically selected to guarantee infinite smoothness near the temporal origin, aggressively isolating the convergence rates of the spatial and temporal schemes without boundary-layer interference. While this isolates analytical accuracy, typical fractional models prominently exhibit weak solution singularities near $t=0$ (e.g., fractional power laws). Adapting the current framework over graded meshes to directly resolve these initial temporal layers represents a straightforward methodological progression.

Unified error norm (Clenshaw–Curtis quadrature).

Since the exact structural solution operates over CGL collocation nodes, analytical verification relies on spectrally converging Clenshaw–Curtis numerical quadrature, which dictates exact integration up to polynomials of degree $\le N$ on the Chebyshev grid [6]:

$$\parallel e_u{\parallel }_{L^2}={\left(\sum _{j=0}^N{w}_j^{\mathrm{CC}}{\left(U_j^{N_t}-u(x_j,T)\right)}^2\right)}^{1/2},$$

(43)

where $w_j^{\mathrm{CC}}$ are the predefined Clenshaw–Curtis weights across $[0,L]$. Validation sets uniformly assess the temporal scaling parameter via

$$\mathrm{Rate}={\displaystyle \frac{log(\parallel e{\parallel }_{\mathrm{prev}}/\parallel e{\parallel }_{\mathrm{curr}})}{log(\Delta t^{\mathrm{prev}}/\Delta t^{\mathrm{curr}})}},$$

(44)

and similarly track spectral exponentiation under variations of N.

4.2. Temporal Convergence

To isolate the temporal discretisation error, the spatial polynomial degree is fixed at $N=24$ (rendering the spectral spatial error negligible), while $N_t$ is successively doubled: $N_t\in \{64,128,256,512,1024\}$. The standard L1 evaluation (without SOE compression) is used to ensure the measured rates reflect only the temporal discretisation. The expected convergence rate is $\mathcal{O}(\Delta t^{2-\gamma })$ for $\gamma \in \{\alpha ,\beta \}$. Three representative fractional order pairs are tested.

Table 2. Temporal convergence for the coupled system at $N=24$, $T=1$.

Nt	Error u	Rate	Error v	Rate
$\alpha =0.1,\beta =0.3$, Expected rates $\approx 1.9$ and $1.7$
64	4.200E-05	—	3.725E-05	—
128	1.019E-05	2.04	9.491E-06	1.97
256	2.535E-06	2.01	2.493E-06	1.93
512	6.387E-07	1.99	6.687E-07	1.90
1024	1.621E-07	1.98	1.821E-07	1.88
$\alpha =0.4,\beta =0.7$, Expected rates $\approx 1.6$ and $1.3$
64	9.274E-05	—	1.281E-04	—
128	2.769E-05	1.74	4.748E-05	1.43
256	8.486E-06	1.71	1.825E-05	1.38
512	2.642E-06	1.68	7.165E-06	1.35
1024	8.305E-07	1.67	2.850E-06	1.33
$\alpha =0.8,\beta =0.9$, Expected rates $\approx 1.2$ and $1.1$
64	6.473E-04	—	3.374E-04	—
128	2.772E-04	1.22	1.518E-04	1.15
256	1.196E-04	1.21	6.958E-05	1.13
512	5.180E-05	1.21	3.217E-05	1.11
1024	2.247E-05	1.20	1.495E-05	1.11

Table 2. Temporal convergence for the coupled system at $N=24$, $T=1$.

Nt	Error u	Rate	Error v	Rate
$\alpha =0.1,\beta =0.3$, Expected rates $\approx 1.9$ and $1.7$
64	4.200E-05	—	3.725E-05	—
128	1.019E-05	2.04	9.491E-06	1.97
256	2.535E-06	2.01	2.493E-06	1.93
512	6.387E-07	1.99	6.687E-07	1.90
1024	1.621E-07	1.98	1.821E-07	1.88
$\alpha =0.4,\beta =0.7$, Expected rates $\approx 1.6$ and $1.3$
64	9.274E-05	—	1.281E-04	—
128	2.769E-05	1.74	4.748E-05	1.43
256	8.486E-06	1.71	1.825E-05	1.38
512	2.642E-06	1.68	7.165E-06	1.35
1024	8.305E-07	1.67	2.850E-06	1.33
$\alpha =0.8,\beta =0.9$, Expected rates $\approx 1.2$ and $1.1$
64	6.473E-04	—	3.374E-04	—
128	2.772E-04	1.22	1.518E-04	1.15
256	1.196E-04	1.21	6.958E-05	1.13
512	5.180E-05	1.21	3.217E-05	1.11
1024	2.247E-05	1.20	1.495E-05	1.11

The results confirm the theoretical $\mathcal{O}(\Delta t^{2-\gamma })$ estimate from Theorem 2: for the u-component, the empirical rate converges to $2-\alpha$, and for the v-component, to $2-\beta$. As expected, the global accuracy is bounded by the slower kernel, $\mathcal{O}(\Delta t^{min(2-\alpha ,2-\beta )})$. For the strongly anomalous regime ($\alpha =0.8$, $\beta =0.9$), the reduced convergence rate ($\approx 1.1$–$1.2$) is clearly observable even at modest $N_t$. For the weakly anomalous case ($\alpha =0.1$, $\beta =0.3$), the observed rates ($\approx 2.0$ and $\approx 1.9$) slightly exceed the theoretical prediction $2-\alpha =1.9$; this is consistent with the fact that the smooth manufactured solution $u=t^{3}sin\left(\pi x\right)$ has no initial-layer singularity, which is the scenario where the $L1$ scheme can exhibit mildly superconvergent behavior for small $\alpha$.

4.3. Spatial Spectral Convergence

To demonstrate the exponential spatial accuracy of the Chebyshev spectral method, we fix the temporal resolution at $N_t=16384$ (reducing the temporal error to $\approx 3.5\times {10}^{-6}$) and vary the polynomial degree $N\in \{4,8,12,16,24,32\}$. The fractional orders are $\alpha =0.4$, $\beta =0.7$.

Table 3. Spatial convergence ($N_t=16384$ fixed, $\alpha =0.4$, $\beta =0.7$, $T=1$). The spectral method saturates at the temporal error floor (†).

N	Error u	Rate	Error v	Rate
4	9.808E-04	—	5.986E-02	—
8	3.287E-06	8.22	4.349E-05	10.43
12	3.481E-06†	−0.14	3.307E-06†	6.35
16	3.481E-06†	0.00	3.302E-06†	0.00
24	3.481E-06†	0.00	3.303E-06†	0.00
32	3.480E-06†	0.00	3.302E-06†	0.00

Table 3. Spatial convergence ($N_t=16384$ fixed, $\alpha =0.4$, $\beta =0.7$, $T=1$). The spectral method saturates at the temporal error floor (†).

N	Error u	Rate	Error v	Rate
4	9.808E-04	—	5.986E-02	—
8	3.287E-06	8.22	4.349E-05	10.43
12	3.481E-06†	−0.14	3.307E-06†	6.35
16	3.481E-06†	0.00	3.302E-06†	0.00
24	3.481E-06†	0.00	3.303E-06†	0.00
32	3.480E-06†	0.00	3.302E-06†	0.00

The spectral method exhibits dramatic exponential convergence: from $N=4$ to $N=8$, the error drops by 2–3 orders of magnitude for the u-component and by 3 orders for the v-component. Already at $N=12$ the spatial error has reached the temporal error floor ($\approx 3.5\times {10}^{-6}$ for $N_t=16384$), and further spatial refinement yields no improvement. The near-zero convergence rates at $N\ge 12$ are not a degradation of the method — they reflect the fact that the total error is entirely dominated by the fixed temporal discretisation. This extreme spatial efficiency is a hallmark of spectral collocation: for smooth solutions, a polynomial degree of $N=8$–12 suffices to saturate any practical temporal resolution.

4.4. Fast SOE Memory Acceleration

To validate the computational efficiency of the SOE fast memory algorithm, we compare the wall-clock execution times of the standard direct $L1$ convolution (cost $\mathcal{O}\left(N_t^2\right)$) against the SOE-accelerated evaluation (cost $\mathcal{O}\left(N_t{N}_e\right)$) with $N_e=40$ exponential terms and a local window of $p=20$. Spatial resolution is fixed at $N=24$.

Table 4. Standard $L1$ vs. fast SOE: wall-clock time and $L^2$ error for the u-component ($N=24$, $\alpha =0.4$, $\beta =0.7$).

Nt	Std. time (s)	Fast time (s)	Speedup	Error ratio
64	0.01	0.01	0.9×	1.00
512	0.13	0.11	1.1×	1.01
1024	0.34	0.22	1.5×	1.01
2048	1.00	0.48	2.1×	1.69
4096	2.38	1.02	2.3×	3.41

Table 4. Standard $L1$ vs. fast SOE: wall-clock time and $L^2$ error for the u-component ($N=24$, $\alpha =0.4$, $\beta =0.7$).

Nt	Std. time (s)	Fast time (s)	Speedup	Error ratio
64	0.01	0.01	0.9×	1.00
512	0.13	0.11	1.1×	1.01
1024	0.34	0.22	1.5×	1.01
2048	1.00	0.48	2.1×	1.69
4096	2.38	1.02	2.3×	3.41

The wall-clock speedup grows from $1.0\times$ at moderate $N_t$ to $2.3\times$ at $N_t=4096$, reflecting the transition from $\mathcal{O}\left(N_t^2\right)$ to $\mathcal{O}\left(N_t{N}_e\right)$ scaling. For $N_t\le 1024$, the error ratio remains within $1\%$ of unity, confirming negligible accuracy loss from the SOE compression. At very large temporal resolutions ($N_t\ge 2048$), the error ratio is observed to increase. This is not an instability, but a strict mathematical consequence of the fixed SOE quadrature precision ($N_e=40$): as $\Delta t\to 0$, the $L1$ temporal truncation error $\mathcal{O}(\Delta t^{2-\alpha })$ drops below the absolute precision floor of the Double Exponential integral approximation. To maintain optimal accuracy for production-scale simulations ($N_t\gg {10}^4$) where the asymptotic linear speedup $\mathcal{O}\left(N_t{N}_e\right)$ becomes massively dominant, the number of exponential nodes $N_e$ must be logarithmically scaled upwards (e.g., $N_e=80\sim 100$).

4.5. Solitary Wave Collision Dynamics

To demonstrate the physical capabilities of the solver on a problem with non-trivial nonlinear dynamics, we remove the manufactured source terms and introduce localised Gaussian initial conditions acting as counter-propagating solitary pulses:

$$u(x,0)=exp\left(-{\displaystyle \frac{{(x-0.3)}^2}{0.005}}\right),v(x,0)=exp\left(-{\displaystyle \frac{{(x-0.7)}^2}{0.005}}\right),$$

(45)

with zero initial velocities and homogeneous boundary conditions. The coupling parameter is set to $\eta =2.0$ to produce a strong inter-component energy exchange during collision. The simulation uses $N=64$ CGL nodes, $N_t=2048$ time steps, and fractional orders $\alpha =0.6$, $\beta =0.4$.

As seen in Figure 2, the Caputo fractional derivatives manifest physically as pronounced memory “tails” trailing behind the propagating wavefronts. The amplitude asymmetry between the u and v wakes is a direct consequence of the different fractional orders ($\alpha \ne \beta$): the u-component, governed by the stronger memory ($\alpha =0.6$), exhibits a more slowly decaying wake due to the stronger hereditary damping, while the v-component ($\beta =0.4$, weaker memory) retains more of its initial amplitude. During the collision event, the inter-component coupling $\eta u{v}^2$ mediates a nonlinear energy exchange that is accurately captured by the implicit scheme without spurious oscillations.

5. Discussion

The proposed method addresses three key computational and analytical challenges of coupled fractional Klein–Gordon systems: the nonlinear iteration cost (resolved by IMEX linearisation), the dense memory bottleneck (resolved by dual-SOE compression), and the rigorous coupled stability–convergence analysis for multi-component systems with disparate fractional orders. We now present a concrete roadmap for extending the method to higher spatial dimensions and discuss the current limitations.

5.1. Extension to Multi-Dimensional Domains

The algorithmic architecture developed in Section 2–3 cleanly separates spatial discretisation from temporal integration and memory management. This modularity enables a systematic extension to higher spatial dimensions. We outline the concrete pathway for $d=2$ on a rectangular domain $\Omega =(0,L_x)\times (0,L_y)$; the extension to $d=3$ or to mapped curvilinear geometries follows analogously via tensor-product or spectral-element constructions [16].

Tensor-product Chebyshev grid. The spatial domain is discretised using $(N_x+1)\times (N_y+1)$ CGL collocation points $(x_i,y_j)$. The two-dimensional Laplacian is discretised via the Kronecker-product structure:

$${\mathbf{D}}_{2\mathrm{D}}^2={\mathbf{D}}_x^2\otimes {\mathbf{I}}_y+{\mathbf{I}}_x\otimes {\mathbf{D}}_y^2,$$

(46)

where ${\mathbf{D}}_x^2\in {\mathbb{R}}^{(N_x+1)\times (N_x+1)}$ and ${\mathbf{D}}_y^2\in {\mathbb{R}}^{(N_y+1)\times (N_y+1)}$ are the one-dimensional Chebyshev second-derivative matrices scaled to $[0,L_x]$ and $[0,L_y]$ respectively, and ${\mathbf{I}}_x,{\mathbf{I}}_y$ are identity matrices of appropriate dimensions.

Dimension-independent temporal components. The Newmark-$\beta$ time-stepping scheme (7), the L1 memory evaluation (8), the IMEX linearisation (16), and the dual-SOE acceleration all operate component-wise on the solution vectors and are therefore entirely dimension-independent. The SOE accumulators $Z_l^n$ increase in length from $N+1$ to $(N_x+1)(N_y+1)$ without any structural modification. Consequently, the dual-SOE compression ratio and the temporal convergence order $\mathcal{O}(\Delta t^{min(2-\alpha ,2-\beta )})$ are preserved identically in any spatial dimension.

Efficient linear algebra via Kronecker structure. The implicit Newmark operator in 2D,

$${\mathbf{A}}_{2\mathrm{D}}={\displaystyle \frac{1}{\Delta t^2}}\mathbf{I}-{\displaystyle \frac{c^2}{4}}\left({\mathbf{D}}_x^2\otimes {\mathbf{I}}_y+{\mathbf{I}}_x\otimes {\mathbf{D}}_y^2\right),$$

(47)

is time-independent, so a single LU factorisation suffices for the entire simulation. For moderate resolutions ($N_x,N_y\le 64$), direct factorisation of the $(N_x+1)(N_y+1)\times (N_x+1)(N_y+1)$ matrix is entirely feasible. For larger problems, the Kronecker-sum structure of (46) is exploited via a generalised Sylvester equation solver: after a one-time eigendecomposition ${\mathbf{D}}_x^2={\mathbf{P}}_x{\Lambda }_x{\mathbf{P}}_x^{-1}$, each time step reduces to solving $N_y+1$ decoupled systems of size $N_x+1$, yielding a per-step cost of $\mathcal{O}(N_x^2{N}_y+N_x{N}_y^2)$ [16].

Stability and convergence in higher dimensions. The a priori estimates of Theorems 1 and 2 extend to the multi-dimensional setting with two modifications: (i) the discrete inner product (19) is replaced by its tensor-product analogue with weights $w_{ij}=w_i^x{w}_j^y$, and (ii) the Chebyshev inverse inequality constant satisfies $C_N=\mathcal{O}\left(N^{d/2}\right)$ in dimension d, affecting the stability constant but not the convergence order. The full convergence estimate $\mathcal{O}(\Delta t^{min(2-\alpha ,2-\beta )}+N^{-s})$ (per coordinate direction) remains valid.

The present one-dimensional study thus provides the rigorous methodological foundation and validation of all algorithmic components — temporal integration, IMEX linearisation, dual-SOE compression, and coupled stability theory — for the multi-dimensional solver currently under development.

5.2. Limitations and Future Directions

While the spectral formulation exhibits exponential spatial convergence for smooth profiles, its direct application to problems with sharp gradients or discontinuities would trigger Gibbs oscillations. In such regimes, spectral filtering or hybridisation with finite element spatial discretisations on the CGL grid, akin to approaches studied in [5], would provide improved monotonicity properties.

From a theoretical standpoint, the a priori stability analysis requires the Chebyshev inverse inequality, embedding a growth factor of $\mathcal{O}\left(N^{d/2}\right)$ in the stability constant. While this is standard for spectral methods applied to nonlinear problems and is not observed to be restrictive in our numerical experiments, developing N-independent stability estimates — for instance via discrete maximum principles or energy-decreasing nonlinear schemes — remains an important open problem.

The dual-SOE compression introduces a quadrature floor ${\epsilon }_{\mathrm{SOE}}$ that becomes visible only when the temporal truncation error drops below $\sim {10}^{-12}$ at $N_{\exp }=40$. For production-scale simulations with $N_t\gg {10}^4$, adaptively scaling $N_{\exp }$ logarithmically ensures the SOE error remains sub-dominant. Finally, the current validation employs smooth manufactured solutions with $u\sim t^3$; resolving the weak initial-layer singularities ($u\sim t^{1+\gamma }$, $\gamma \in (0,1)$) inherent in realistic fractional models requires graded temporal meshes near $t=0$, which represent a natural and straightforward extension of the present framework.

6. Conclusions

We have developed a fast, fully discrete IMEX Chebyshev spectral collocation method for coupled nonlinear Klein–Gordon systems with multi-order Caputo fractional memory. The method combines four key algorithmic ingredients: (i) implicit Newmark averaging of the spatial operator for unconditional stability, (ii) explicit linear extrapolation of the nonlinear terms to eliminate Newton–Raphson iterations while preserving $\mathcal{O}(\Delta t^2)$ accuracy, (iii) explicit L1 evaluation of the fractional memory sums, and (iv) dual-SOE compression of the disparate memory tails to reduce the temporal complexity from $\mathcal{O}\left(N_t^2\right)$ to $\mathcal{O}\left(N_t(p+N_{\exp })\right)$. The resulting scheme requires only a single LU factorisation and two back-substitutions per time step, with no CFL restriction on $\Delta t$ at fixed spatial resolution.

Rigorous a priori stability estimates and a global convergence bound of $\mathcal{O}(\Delta t^{min(2-\alpha ,2-\beta )}+N^{-s}+{\epsilon }_{\mathrm{SOE}})$ have been established using a unified discrete energy technique and the Gronwall inequality, treating both coupled components and their inter-component nonlinear interactions simultaneously. The theoretical rates are confirmed numerically across three distinct fractional order pairs, and the spectral spatial discretisation achieves exponential convergence with as few as $N=8$–12 modes. A solitary wave collision scenario further validates the method’s ability to capture the asymmetric dispersive wakes characteristic of multi-order fractional memory.

Crucially, the algorithmic architecture — temporal integration, IMEX linearisation, dual-SOE memory compression, and the coupled stability analysis — is dimension-independent by construction. We have provided a concrete extension pathway to multi-dimensional tensor-product Chebyshev grids (Section 5.1), demonstrating that only the spatial discretisation matrix changes (via Kronecker products and Sylvester-based solvers), while all temporal and memory components carry over identically. The full multi-dimensional implementation and validation on physically motivated two-dimensional coupled viscoelastic wave problems constitute our immediate ongoing work, building directly on the methodological foundation established here. Additional extensions to graded temporal meshes for resolving initial-layer singularities and to spectral-element formulations for complex geometries are natural progressions of the present framework.