This paper develops a structured analytical framework and a robust numerical methodology for the spectral stability of time-varying bilateral quaternion differential equations of the form q˙=A(t)q+qB(t). By systematically extending classical real matrix theory to non-commutative dynamical systems via exact isometric real representations, this study utilizes the Kronecker product of real adjoint matrices to rigorously elucidate the underlying tensor structure of the bilateral evolution operator. This tensor-based reformulation proves that the Floquet multipliers of the bilaterally coupled system can be strictly decoupled into the product of the spectra corresponding to the left and right unilateral subsystems. Second, a “Scalar-Vector Stability Separation Principle” based on logarithmic norms is proposed, demonstrating that the transient energy evolution of the system is governed exclusively by the Hermitian real parts of the coefficient matrices, remaining entirely independent of the anti-Hermitian imaginary parts (rotation terms). Furthermore, for constant-coefficient and slowly varying systems, the Riesz projection from holomorphic functional calculus is introduced to establish algebraic criteria for exponential dichotomies, thereby revealing a cubic scaling law that relates the robustness threshold to the spectral gap (ε0∝β3). Numerically, a Quaternion Chebyshev Spectral Collocation Method (Q-CSCM) is embedded within this exact vectorization framework to ensure that the algebraic symmetries of the bilateral system are strictly preserved through the isomorphic mapping. By explicitly constructing the fully discrete Kronecker product matrix via the exact real vectorization isomorphism, discrete energy estimates are utilized to rigorously prove that the numerical scheme successfully inherits the intrinsic spectral accuracy of the Chebyshev approximation. Comprehensive numerical experiments demonstrate that, within the low-dimensional regime, this methodology exhibits substantial temporal approximation efficiency advantages and superior numerical robustness compared to an alternative Legendre spectral baseline, as well as traditional explicit and state-of-the-art implicit symplectic Runge–Kutta methods, particularly when solving stiff and critically stable problems such as nonlinear Riccati oscillators.
Si et al. (Sat,) studied this question.