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

We develop a fast Chebyshev spectral collocation method for a coupled system of nonlinear Klein–Gordon equations augmented by Caputo-type fractional memory integrals. The governing equations retain the classical second-order time derivative as the leading operator and incorporate weakly singular convolution integrals with power-law kernels t−α, α∈(0,1), modelling viscoelastic memory damping rather than replacing the wave operator. The spatial discretisation employs Chebyshev–Gauss–Lobatto collocation, while the temporal integration uses a Newmark scheme (βNM=1/4) with the spatial operator treated implicitly and both the L1 memory sums and the cubic nonlinearities evaluated explicitly at the known time level; a linear extrapolation of the nonlinear terms eliminates the need for Newton–Raphson iterations. The disparate memory tails arising from two distinct fractional orders α≠β are compressed by independent Sum-of-Exponentials (SOE) approximations, reducing the per-step memory cost from O(Nt) to O(p+Nexp) and the total complexity from O(Nt2) to O(Nt(p+Nexp)). A rigorous stability estimate and a global convergence bound are established using a discrete Gronwall inequality. Numerical experiments confirm the temporal convergence rate O(Δtmin(2−α,2−β)), spectral spatial accuracy, and the practical speedup afforded by the SOE acceleration. A solitary wave collision scenario illustrates the method’s capability to capture asymmetric dispersive wakes generated by the fractional memory. The algorithmic architecture is dimension-independent by construction; a concrete extension pathway to multi-dimensional tensor-product Chebyshev grids, including Kronecker-product operators and Sylvester-based solvers, is presented.