Abstract <p>This study employs high-precision numerical integration methods to analyze the phenomenon of coexistence of weak chaos and strong stability in the Sun–Earth–Moon system, taking into account factors such as tidal perturbations, planetary rotation, and relativistic corrections. Key chaos indicators were calculated, including Lyapunov exponents, Poincaré sections, and fractal dimensions, to quantify the system’s dynamic behavior. For the Sun–Earth system, the estimated Lyapunov exponent is on the order of 10^–6 per year, while the Earth–Moon system exhibits stronger chaotic characteristics with a Lyapunov exponent of approximately 10^–4 per year, corresponding to an orbital divergence time scale shorter than million years. Poincaré section analysis and the box-counting method reveal that the fractal dimension of the Earth–Moon subsystem is 1.768, indicating the existence of a mixed phase space where KAM tori and chaotic resonance layers coexist. Comparison with an idealized two-body model validates the numerical accuracy and highlights the destabilizing effect of solar perturbations. Furthermore, the influence of other planets is quantitatively proven to be negligible. Long-term stability analysis (over 10⁴ years) shows that variations in eccentricity and the Earth–Moon distance remain within limited ranges, with tidal dissipation acting as a stabilizing mechanism. The drift in energy and angular momentum is minimal, verifying the reliability of the numerical results. Multi-scale entropy analysis further supports the conclusion of weak chaos, showing low overall complexity but dynamic structures that vary over time.</p>