## 2 Methods  
Figure 1 shows the vibro-impact model with mass, linear spring, and nonlinear elastic contact. The spring \(k\) and \(k_{ex}\) are linear spring coefficients, and the system can be under external kinematic excitation forcing. The spring \(k_{ex}\) at the left end has an amplitude \(f_{0}\) and an excitation frequency \(\omega_{ex}\), and is excited harmonically. The distance between the impact  
surface at the static equilibrium position of the system is \(r\). In addition, the nonlinear elastic model of the collision force was defined as [PERSON]'s contact model [12]. The dynamic behavior of the system can be made dimensionless for the generalization and efficiency of interpretation, as shown in Figure 1b. Consider a mass m attached to a spring stiffness coefficient \(k\), \(k_{e}\) and coefficient of viscous damping \(c\). If the impact force is 0, it is only a single degree-of-freedom equation of motion by excitation. The equation of motion is as follows:  
\[m\ddot{x}+c\dot{x}+(k+k_{ex})x=f_{0}k_{ex}\cos\omega_{ex}t. \tag{1}\]  
Using the dimensionless time, \(\tau=t\sqrt{(k+ke)/m}\) and the coordinate transformation \(x(t)=X(t)\), the dimensionless equation of motion can be written as:  
\[X^{\prime\prime}+X^{\prime}+X=\xi f_{0}\cos\eta\tau \tag{2}\]  
where prime is the differentiation with respect to \(\tau(\geq 0)\), and \(\omega_{n}=(\sqrt{(k+k_{e})/m})\) is the natural frequency of the undamped system. \(\eta=\omega_{ex}/\omega_{n}\) is the relation between the excitation and natural frequency. \(d=c/\sqrt{(k+k_{e})m}\) is the dimensionless damping ratio, \(\xi=k_{e}/(k+k_{e})\) is the relation between the spring coefficient, and \(u=X/f_{0}\) is the dimensionless displacement related to the magnitude of excitation forcing,  
\[u^{\prime\prime}+u^{\prime}+u=\xi\cos\eta\tau. \tag{3}\]  
The most popular [PERSON]'s model of contact is the nonlinear elastic model, and the dimensionless collision force is given by:  
\[F_{L}=0\text{ if }u<L \tag{4}\] \[F_{L}=k_{c}(u-L)^{3/2}\text{ if }u\geq L\]  
Figure 1: Dynamical impact system with external excitation, (**a**) 1-D impacting oscillator, (**b**) dimensionless form of (**a**). where \(k_{c}=k_{col}\sqrt{f_{0}}/(k+k_{e})\) is the dimensionless ratio of the stiffness of the surface, which is dependent on the elastic properties and geometry of the rigid body. The dimensionless equation of motion with impact force can be written as:  
\[u^{\prime\prime}+u^{\prime}+u=\xi\cos\eta\tau+F_{L}H(u-L)^{3/2} \tag{5}\]  
where \(H\) is the Heaviside function. For numerical integration, Equation (4) can be rewritten in the state space such that:  
\[ \begin{array}{l}\dot{z}_{1}=z_{2}\\ \dot{z}_{2}=-z_{2}-z_{1}+\xi\cos z_{3}+F_{L}H(u-L)^{3/2}\\ \dot{z}_{3}=\eta\end{array} . \tag{6}\]  
The dynamical equation of motion can be converted to the state-space vector form as follows:  
\[\dot{\mathbf{z}}\mathbf{=f}(\mathbf{z}),\mathbf{z}(t_{0})\mathbf{=z}_{0}. \tag{7}\]  
Here, \(\mathbf{y}(t_{0})\) is the initial condition of the perturbed solution, and Equation (7) is taken for the Taylor series at \(\mathbf{z}(t)\) such that:  
\[\delta\dot{\mathbf{z}}\mathbf{=\frac{\partial\mathbf{f}}{\partial\mathbf{z} }\bigg{|}_{\mathbf{z}=z_{i}}\cdot\delta\mathbf{z}+O_{(2)},\delta\mathbf{z}(t _{0})=\delta\mathbf{z}_{0}. \tag{8}\]  
By letting \(\delta\mathbf{z}=[\mathbf{\Phi}_{t}(\mathbf{z}_{0})]\delta\mathbf{z}_{0}\) substitute into perturbation Equation (8), the variation equation can be written as:  
\[\Big{[}\dot{\mathbf{\Phi}}_{t}(\mathbf{z}_{0})\Big{]}\delta\mathbf{z}_{0}=[D _{z}\mathbf{f}][\mathbf{\Phi}_{t}(\mathbf{z}_{0})]\dot{\mathbf{z}}\mathbf{z}_ {0},[\mathbf{\Phi}_{t_{0}}(\mathbf{z}_{0})]\dot{\mathbf{z}}\mathbf{z}_{0}=[ \mathbf{I}]\dot{\mathbf{z}}\mathbf{z}_{0} \tag{9}\]  
where \([D_{z}\mathbf{f}]\), \([\mathbf{I}]\), and \([\mathbf{\Phi}_{t}(\mathbf{z}_{0})]\) denotes the Jacobian matrix, identity matrix, and solution of the variational equation, respectively.