Stability and responses of parametrically excited systems using higher-order piecewise polynomial approximations

The stability and responses of dynamic systems subjected to parametric excitation are frequently encountered in many fields of engineering and applied science, such as pile foundations under axial loadings. This paper presents a new numerical simulation method for constructing stability diagrams and predicting the responses of such dynamic systems governed by Mathieu-Hill equations. In the proposed method, the excitation period is divided into time intervals, and the excitation function within each interval is approximated using constant, linear, or quadratic piecewise polynomial representations. The corresponding local state transition matrices are obtained analytically using trigonometric, Bessel, and Weber functions, respectively, and are assembled over one excitation period. The method is applied to a pile foundation, modelled as a beam on a Winkler foundation, to develop stability diagrams and evaluate the dynamic responses. The results show that the constant approximation may fail to capture certain instability regions when a small number of time intervals is used, whereas the linear and quadratic approximations improve the accuracy of the stability prediction. In particular, the quadratic approximation identifies narrow instability regions using fewer intervals and predicts the critical dynamic loads of the first three instability regions with differences of less than 1% relative to the higher-order Hill-ID results. The time-history responses are also compared with the classical RKF45 solution, with the quadratic approximation giving the coefficient of determination of 0.9986 at the selected unstable point. The convergence results further confirm that higher-order approximations require fewer intervals to provide converged predictions. The proposed approach provides an accurate framework for the stability and response analysis of parametrically excited systems. The results also provide a calibration of the conventional Hill infinite determinant method in terms of harmonic balance.