## 4 Results and Discussion  
### The HRVD Oscillator Model with a Periodic External Excitation  
The hybrid van der Pol-Duffing-Rayleigh (HRVD) oscillator accurately describes the lateral forces generated during walking, whose parameter variations can destabilize the system, leading to bursting oscillations and sway. Researching such systems aids in reducing economic and injury risks from engineering vibrations. We start by applying the HAVOK analysis framework for bursting oscillations to an HRVD oscillator [9] with an external periodic excitations model governed by  
\[\ddot{x}+\left(a_{3}\dot{x}^{2}+\alpha x^{4}+\beta x^{2}+\lambda\right)\dot{x}+ a_{1}x+a_{2}x^{3}=f_{1}\cos(\omega_{1}t). \tag{23}\]  
For Parameter Set \(\mathrm{Para1}=(a_{1},a_{2},a_{3},\alpha,\beta,\lambda,\omega,f_{1})=(-5,5,0.1, 0.1,-1,1.2,0.1,3)\), the uniform sampling time interval is set to \(dt=0.001\), with the time range spanning \([0,1000]\). The initial conditions are set as \(\left(x,\dot{x}\right)=(0.1,0.1)\). The state training data are obtained using the Runge-Kutta numerical solver in MATLAB R2022a. To mitigate the effects of initial chaos and focus on the system's steady-state or periodic behavior over a longer period, the first 5000 data points are discarded during the numerical simulation. Since the model (23) is a one-dimensional system, the embedding dimension should be at least 3. We initially set the embedding dimension \(m=3\) to observe the singular values (Sigs) and plot the variation of \(E(k)\) with \(k\), as shown in Figure 1.  
\[Sigs=\Big{[}1.8\times 10^{4},1.2,1.1\times 10^{-3}\Big{]}. \tag{24}\]  
According to model (24) and Figure 1, when \(r=2\), \(E(k)\) can exceed 99% and approach the convergence value. Choosing \(r=2\) implies that in this example, we can only plot the two-dimensional embedded flow field and the one-dimensional reconstructed flow field (since with \(r=2\), the output matrix of the forced linear model consists of only one column vector). By plotting the fit between the first column of the output matrix and the first column of the input matrix (Figure 2), we can observe the reconstruction performance of the forced linear model. Additionally, the mean squared error (MSE) of the reconstruction at this point is \(1.3751\times 10^{-17}\). Figure 1: Plot of cumulative energy variation with the number of singular values when the parameter set is taken as \(\mathrm{Para1}\), \(\mathrm{m}=3\), \(\mathrm{r}=2\). Alternatively model (23) can be rewritten as  
\[\left\{ \begin{array}{l}\dot{x}=y\\ \dot{y}=-a_{2}x^{3}-a_{1}x-a_{3}y^{3}-\alpha x^{4}y-\beta x^{2}y-\lambda y+f_{1} \text{cos}(\omega_{1}t),\end{array} \right. \tag{25}\]  
This treatment allows the researcher to plot the 2D phase diagram (Figure 3b). This allows us to comparatively examine the reconstruction of the 2D-embedded flow field in the 2D time-lag embedding space. Figure 3: Series of experimental plots when the parameter set is taken to be Para1, m = 3, r = 2. (**a**) 2D-embedded flow field; (**b**) Original 2D phase diagram; (**c**) Upper subplot for the eigentime series \(v_{1}\) and lower subplot for the forced eigentime series \(v_{r}\); (**d**) Localized enlargement of (**c**). Figure 2: The output matrix of the forced linear model when the parameter set is Para1, m = 3, r = 2 is plotted against a fit to the first column of the output matrix, where only uniformly spaced (100 \(\Delta\)t) discrete trajectory points are retained in the first column of the output matrix, labeled as blue points to avoid an overly dense and complex image.