Abstract
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 modelling viscoelastic memory damping. The spatial discretisation employs Chebyshev–Gauss–Lobatto collocation, while the temporal integration uses a Newmark scheme (βNM=1/4) combined with an implicit–explicit linearisation in which the linear spatial operator is treated implicitly and the nonlinear terms are treated explicitly through a second-order extrapolation. This linearisation eliminates the need for Newton–Raphson iterations at each time step. To overcome the dense memory bottleneck arising from two distinct fractional orders α≠β, the convolution memory kernels are compressed by independent sum-of-exponentials approximations obtained from a double-exponential quadrature of the kernel’s integral representation, which significantly reduces the computational complexity of the history term. A rigorous stability estimate and a global convergence bound are established using a discrete Grönwall inequality. Numerical experiments confirm the theoretical temporal and spatial convergence rates and demonstrate the practical speed-up afforded by the sum-of-exponentials acceleration. A solitary wave collision scenario illustrates the method’s capability to capture asymmetric dispersive wakes generated by the fractional memory.