Numerical calculation and accuracy analysis of the highly oscillatory integral in the complete geodetic excitation solution of Earth's polar motion equation

Variations in polar motion (PM) observed through space geodetic techniques imply comprehensive information on geophysical excitation. Traditional geodetic excitation series derived from PM observations typically excluded the calculation of the excitation derivative term d ψ / d t in the classical Liouville equation (LE) and relied on the low-accuracy numerical integration method. In this paper, we derive the complete LE (CLE) geodetic excitation solution, in which the integral term involves a highly oscillatory integral with a daily-frequency oscillatory kernel e i Ω τ . To address the limitations of traditional numerical integration methods when applied with daily PM observations as input, we apply two Filon-type methods—one based on Lagrange polynomials and the other on cubic spline interpolation to derive the CLE geodetic excitation. Using both simulated and observed datasets, we perform a recovery process for CLE solutions (input PM → excitation → recovered PM) to assess numerical calculation errors under the accuracy standard in the PM domain. The observed dataset demonstrates that the recovered PM errors, with standard deviations (STD) greater than 30 microarcseconds (μas) for the X and Y components, primarily vary within tens of μas, comparable to the observed PM errors during the same period. This also proves that the derived CLE geodetic excitation, with the geophysical completeness, satisfies the analysis research under the current PM observational accuracy, providing a valuable reference indicator for various geophysical processes. The differences between CLE and traditional geodetic excitations, with STDs of about 5–6 milliarcseconds (mas), arise from high-frequency biases associated with numerical calculations. Regarding future PM investigations in daily and sub-daily frequency bands, simulation tests indicate that a time sampling interval of less than 4 h for the input PM can achieve accuracies below the μas level, which demonstrates the capability of CLE numerical calculations to support future high-resolution PM applications.