A high-order Rothe–Jacobi framework for multi-term fractional telegraphic models: error estimation, stability, and precision performance

Abstract This paper presents a high-fidelity numerical framework, the Rothe–GSJP algorithm, designed to resolve the multi-term fractional telegraph equation (FTE) subject to inhomogeneous initial-boundary conditions (IBCs). Unlike traditional Rothe-type methods, which typically address single-rate evolutions, the proposed approach integrates a specialized hybrid temporal semi-discretization utilizing L1 and L2 stencils to simultaneously handle multi-scale temporal memory effects, with a spatial spectral projection based on generalized shifted Jacobi polynomials (GSJP). By leveraging the operational matrices (OMs) of fractional derivatives (FDs), the governing partial differential equation (PDE) is transformed into a sequence of steady-state spatial fractional systems subject to inhomogeneous BCs. A central contribution of this work is the rigorous stability analysis conducted via an operator-theoretic approach. We establish the unconditional stability of the scheme in the $L_{\infty }$ L ∞ norm by proving that the system operator is strictly accretive and that its inverse is bounded by the reciprocal of the minimum temporal damping coefficient, ρ . This approach successfully circumvents the mathematical limitations of the traditional maximum principle, which is often inapplicable to non-local fractional spatial operators. Furthermore, we provide a unified convergence analysis demonstrating that the global error is bounded by the summation of the temporal truncation error— $\mathcal{O}((\Delta \tau )^{\min (2-\alpha , 3-\beta )})$ O ( ( Δ τ ) min ( 2 − α , 3 − β ) ) —and the spatial spectral decay rate, $N^{-s}$ N − s . By utilizing GSJP basis functions that incorporate BCs directly at the functional level through an inhomogeneous boundary lifting technique, the framework avoids the approximation saturation common in standard spatial methods. Numerical simulations demonstrate that the Rothe–GSJP algorithm consistently suppresses numerical dispersion and achieves a maximum absolute error (MAE) in the range of 10 −14 to 10 −15 . This near-machine precision floor represents a significant improvement over existing methods reported in the literature, which often saturate at significantly higher error margins. The results confirm that the framework is robust across the entire fractional spectrum and effectively suppresses numerical dispersion, making it a powerful tool for the precision modeling of fractional wave-diffusion phenomena.