Stability Analysis of Fractional‑Order Mathieu Equation with Forced Excitation

: The advantage of fractional‑order derivative has attracted extensive attention in the field of dynamics. In this paper, we investigated the stability of the fractional‑order Mathieu equation under forced excitation, which is based on a model of the pantograph–catenary system. First, we obtained the approximate analytical expressions and periodic solutions of the stability boundaries by the multi‑scale method and the perturbation method, and the correctness of these results were verified through numerical analysis by Matlab. In addition, by analyzing the stability of the k’T‑ periodic solutions in the system, we verified the existence of the unstable k’T‑resonance lines through numerical simulation, and visually investigated the effect of the system parameters. The results show that forced excitation with a finite period does not change the position of the stability boundaries, but it can affect the expressions of the periodic solutions. Moreover, by analyzing the properties of the resonant lines, we found that when the points with k’T‑periodic solutions were perturbed by the same frequency of forced excitation, these points became unstable due to resonance. Finally, we found that both the damping coefficient and the fractional‑order parameters in the system have important influences on the stability boundaries and the resonance lines.


Introduction
Since the concept of fractional-order derivative was first formulated by Hospital and Leibnitz in 1695 [1], much research on the definition, operation, and application of fractionalorder derivative has been conducted [2][3][4][5][6].Gradually, Fractional-order derivative was used to deal with numerous engineering problems, due to its indispensable role in the research of dynamic behavior, system optimization, and other engineering problems [7][8][9][10][11].
The applications of fractional-order derivative in engineering problems are mainly in the fields of dynamics and control [12][13][14][15][16].In the field of dynamics, fractional-order derivative is commonly used to model engineering materials with memory properties, such as viscoelastic components.By defining more accurate constitutive relations of materials, the vibration characteristics of nonlinear systems, such as air springs in the pantograph systems, can be more reasonably analyzed.For instance, Qi et al. [12] used fractional-order derivative to describe the viscoelastic damping characteristics of rubber airbags and optimized the aerodynamic model by analyzing the established fractional-order-modified model.In addition, fractional-order derivative can be very efficient in discrete systems.Lu et al. [13] constructed a new fractional-order discrete memristor model with prominent nonlinearity and simulated the dynamical behaviors of the Ruelkov neuron under electromagnetic radiation by introducing the proposed discrete memristor.Chen et al. [14] proposed a three-dimensional fractional-order discrete Hopfield neural network (FODHNN) within the left Caputo discrete delta's sense and studied the dynamic behavior and synchronization of FODHNN.
In the field of control, the fractional-order derivative is mainly used to influence the closed-loop control properties to improve the robustness of the system.Hou et al. [15] established the dynamic model of a gear transmission system based on speed feedback fractional-order proportional-integral-derivative (PID) control and discussed the robustness of the fractional-order PID controller in meshing the stiffness coefficient and excitation amplitude.Chang et al. [16] improved the traditional linear suspension model, studied the active control of the suspension with nonlinear stiffness and fractional-order damping, and verified the rationality of the model and the effectiveness of fractional-order PID control.The parametric excitation system has been widely used in engineering [17][18][19][20].For a parametrically excited system containing viscoelastic devices, such as the pantographcatenary system, the use of the Mathieu equation can adequately explain some resonance phenomena in the system, and the Mathieu equation with a fractional-order term can describe the constitutive relation of the system more accurately.In recent years, there has been a great deal of research on the fractional-order Mathieu equation.For example, Li et al. [21] constructed the fractional-order Mathieu equation for viscoelastic simply supported beams and analyzed the influence of the axial periodic excitation on the system dynamics by the averaging method.Wen et al. [22] studied the dynamic characteristics of the fractional-order Mathieu equation with damping via the multiple-scale method and discussed the influence of the fractional-order parameter on the stability boundaries.Guo et al. [23] introduced fractional-order derivative theory into the quasi-periodic Mathieu equation and studied the effects of fractional-order term parameters on the stability of the equation.In addition, several authors studied the influence of forced excitation factors on parametric excitation systems.For example, Aurora et al. [24] discussed the stable boundaries and unstable periodic solutions for the non-homogeneous Hill equation.
Note that most of the existing studies focused on the stability analysis of the homogeneous fractional Mathieu equation.However, due to the complexity of the actual engineering environment, the influence of external excitation must be considered.At the same time, the forced excitation will lead to some obvious changes in the stability of the fractional Mathieu equation, and some new unstable resonance phenomena will appear.Therefore, a more comprehensive stability analysis of the inhomogeneous fractional Matheiu equation is necessary.Motivated by the above considerations, the novelty and contribution of this paper are summarized as follows: (1) The multi-scale method combined with the perturbation method to obtain analytical solutions for the stable boundaries and periodic solutions was used; (2) Through theoretical analysis and numerical simulation by Matlab, we found that under certain circumstances, the forced excitation can lead to the divergence of the originally stable periodic solution in the stable region, which can constitute the resonance line.(3) We found that the parameters of fractional-order differential term affect the stability boundary and resonance line in the form of equivalent linear stiffness and damping.
In this paper, the system modeling, the stability determination theorem, and the analytical solutions of the stability boundary are discussed in Section 2. The accuracy of the analytical solutions result is verified by Matlab numerical analysis in Section 3. In Section 4, the existence conditions and properties of the unstable periodic solutions (resonance line) in the system are described, and the influence of fractional-order parameters is analyzed.Finally, some conclusions are provided in Section 5.

Stability Analysis of Fractional-Order Parametrically Excited System
The vibration problem is common in a high-speed vehicle system [25].Because the spring and damper in the vehicle have obvious fractional-order characteristics, it is necessary to use a fractional-order model instead of an integer-order model to describe the memory characteristics of such materials.In this section, we take the pantograph-catenary system model as an example to analyze the stability of the fractional-order parametrically excited system.

The Fractional-Order Model of the Pantograph-Catenary System
The fractional-order model is shown in Figure 1.It contains an air spring with fractional characteristics [26].
The vibration problem is common in a high-speed vehicle system [25].Because the spring and damper in the vehicle have obvious fractional-order characteristics, it is necessary to use a fractional-order model instead of an integer-order model to describe the memory characteristics of such materials.In this section, we take the pantograph-catenary system model as an example to analyze the stability of the fractional-order parametrically excited system.

The Fractional-Order Model of the Pantograph-Catenary System
The fractional-order model is shown in Figure 1.It contains an air spring with fractional characteristics [26].
The fractional-order model of the pantograph-catenary system.
In the simplification process of the pantograph-catenary system model, the pantograph is equivalent to the mass block, m , and the catenary force applied to the pantograph is equivalent to the variable stiffness spring, In Equation (1), ) (t F w stands for the excitation of the train on the pantograph, c is the damper, and is the air spring with fractional characteristics.According to this model, the fractional dynamic equation of the pantograph-catenary system can be constructed: Substituting , and Equation (3) can be transformed to It is obvious that the Equation ( 4) is a fractional-order Mathieu equation under forced excitation, which is the object of the following analysis.

Extended Floquet Theory
Linear homogeneous differential equations with T-periodic coefficients are given by In the simplification process of the pantograph-catenary system model, the pantograph is equivalent to the mass block, m, and the catenary force applied to the pantograph is equivalent to the variable stiffness spring, k(t), where In Equation ( 1), F w (t) stands for the excitation of the train on the pantograph, c is the damper, and KD p [x(t)] is the air spring with fractional characteristics.According to this model, the fractional dynamic equation of the pantograph-catenary system can be constructed: Substituting τ = πvt L into Equation (2), yields (3) , and Equation (3) can be transformed to It is obvious that the Equation ( 4) is a fractional-order Mathieu equation under forced excitation, which is the object of the following analysis.

Extended Floquet Theory
Linear homogeneous differential equations with T-periodic coefficients are given by where A(t) de f = A(t + T), ∀t ∈ R 1 is time-periodic and T is the time period, which can be evaluated by applying the Floquet-Lyapunov theorem [27].
According to the Floquet-Lyapunov theorem, any fundamental matrix X(t) with Tperiodic coefficients in Equation (5) can be written as where q(t) is a continuous T-periodic n × n matrix-function, and µ represents the Floquet characteristic exponents, which determine the stability of the trivial solution to Equation (5).
In addition, Equation (6) satisfies X(t + T) = λ X(t), where λ represents the Floquet multipliers, and the relationship between λ and µ can be expressed as In terms of the periodicity of X(t) and q(t), we can obtain the stability condition for system (6): 1.

3.
If Re(µ) = 0, |λ| = 1, then the solution is stable without the interference of forced excitation, and when there exists a positive integer m, also satisfying λ k ′ = 1, we can obtain the solutions with k'T-periodic.
If the linear system is inhomogeneous (with a forcing excitation), then Equation (5) should be changed to where g(x) = g(x + T) is a time-varying forcing excitation.If X(0) initial conditions are given and bounded, then the solutions to the inhomogeneous system of Equation ( 8) can be expressed as follows When we analyze the inhomogeneous linear system, the stability conditions 1 and 2 mentioned above do not change.However, condition 3 requires further refinement and adjustment: 3.1 Based on |λ| = 1, when λ 1 = λ 2 = ±1, the Floquet multipliers correspond to the πperiodic solutions and 2π-periodic solutions (corresponds to the stability boundary).
For an inhomogeneous linear system, the forced excitation with finite period does not change the convergence of the stability boundary; only when the period of the forced excitation is large enough will the stability boundary become unstable.
The solutions to the inhomogeneous system after nT time periods can be expressed as follows: In this case, the stability of the solutions depends on the convergence characteristics of the summation when it exists in the Cesaro [28] limit In the limit, the summation term becomes a diverging algebraic series increasing by G with each additional term.Therefore, the summation term causes the result to diverge.The solution is given by where Λ is the integral term and is a constant vector.
3.2 When the Floquet multipliers are λ 1 = λ 2 , |λ 1 | = |λ 2 | = 1, they correspond to the periodic solutions in the system.Under the forced excitation with the same period, such k'T-periodic solution will become unstable due to resonance [28].This will be described in more detail in Section 4.1.3.3 When there are multiple eigenvalues, not semisimple, the solutions of the system are unbounded in both homogeneous and inhomogeneous cases.

The Analytical Study of Stability Boundaries
In this section, the analytical solution of the stability boundary of Equation ( 4) is discussed.When the forced excitation f (t) is a harmonic excitation, such as F cos(wt), Equation (4) can be transformed into where ξ(0 ≤ ξ << 1) is the linear damping of the system, 2ε cos(2t) is the time varying stiffness of the system, KD p [x(t)] is the p derivative of x(t) with respect to t, and K(K > 0) is the coefficient of the fractional-order derivative.
There are several definitions for fractional-order derivative, such as Riemann-Liouville and Caputo definitions [1][2][3]29] and some newly proposed definitions, such as the generalized fractional derivative (GFD) definition [30,31].Under wide senses, they are equivalent for most mathematical functions.Without generality, Caputo's definition is adopted in Equation ( 13) with the following form: where Γ(y) is the Gamma function, satisfying the relation: into Equation (13), we obtain ..
According to the multi-scale method [22], introducing time variables of different scales, we can obtain the following expansions The approximate solutions of Equation ( 13) can be set as Combined with the perturbation method, δ is assumed to have the following form: If Equations ( 16) -( 18) are substituted into Equation ( 15) and the coefficients ε of the same power are compared, a system of equations can be obtained: The periodic solution of Equation ( 19) is where δ 0 = n 2 , n = 0, 1, 2, . .., and cc represent the complex conjugation of its preceding terms.The stability boundaries and periodic solutions for δ 0 = 0, 1, 4 are discussed below.
To avoid the secular term, aδ 1 must be 0, so Thus, the particular solution of Equation ( 21) can be written as follows: Use Equation D p t e iλt = (iλ) p e iλt and Euler formula, and the fractional-order differential term can be transformed into the following: Consequently, by substituting Equations ( 23), (25), and (26) into Equation ( 21), we can obtain the following: ..
By eliminating the secular term in Equation ( 28), one would arrive at δ 2 = − 1 2 and the particular solution that is Therefore, the stability boundary near δ 0 = 0 is given by There is a periodic response with a period of 2π in this boundary, and it is In order to find the relationship between fractional-order parameters and the stability boundary, we need to carry out higher-order approximate calculation, as follows: Similarly, by eliminating the secular term, we will have Hence, the stability boundary for this case can be written as Choose Substitute δ 0 = 1 into Equation ( 19) and it could yield By substituting δ 0 and x 0 into Equation (20), we arrive at where D 1 A is the partial derivative of A with respect to T 1 and A is the conjugate of A.
In order to obtain a periodic solution, we must let To solve Equation (37), set A = a 2 − b 2 i, separate the real and imaginary parts, and it can be obtained: If the system of equations has a non-zero solution, we can obtain the following: Hence, the particular solution of Equation ( 36) is Similarly, eliminating the secular terms of Equation ( 21), we obtain Using the above equations, the stability boundaries near δ 0 = 1 can be written as follows: where the equivalent linear damping is Using a similar procedure, we can obtain the stability boundaries near δ 0 = 4, which are where the equivalent linear damping is . Its corresponding periodic solution is presented: The accuracy of the above analytical results is verified through the analysis of special cases.First, when the forced excitation F = 0, the calculated results of stability boundaries and periodic solutions are the same as those of the homogeneous fractional-order Mathieu equation [22].Second, when the fractional -order derivative terms are removed from the results, the calculation results are consistent with those of the typical Mathieu equation.
According to the above results of δ 0 = n 2 , n = (0 , 1, 2), it can be found that both damping and fractional-order derivative in Equation ( 4) affect the stable boundaries and periodic solutions of the system, and their change affects the position of the stable boundaries and the form of the periodic solutions at the same time.Meanwhile, the magnitude and frequency of the forced excitation with a definite period do not affect the stability of the stable boundaries, but they will change the properties of the periodic solutions of the stable boundaries.

Comparative Analysis of Numerical Results and Analytical Results
To verify the accuracy of the analytical results, we numerically tested the influence of forced excitation on the stability boundaries of Equation (4) by the numerical method.First, the stability boundaries were drawn according to Equations (34), (43), and (45), as shown in Figure 2. The three shaded parts in Figure 2 represent the unstable region near δ 0 = 0 , 1, 4 respectively, and the blank region represents the stable region of the system.The boundaries of the intersection of the two regions are the stability boundaries of parametric excitation resonance.
change.Therefore, some typical points for numerical simulation were selected in the bottom area of the stability boundaries and near the boundary lines of  The responses of each typical point were compared and verified by the numerical method.The system parameters were selected as 0.005 , and 1 w = .Figures 3-6 show the time-domain response of points A, B, C, and D respectively, in which sub-graph (a) is the response without forced excitation and sub-graph (b) is the response with forced excitation.It can be seen from the comparison diagrams in Figures 3-6 that whether the point is in the stable region or the unstable region, the forced excitation did not change the stability of each point, but changed the speed of system convergence and divergence.This is consistent with the analytical results obtained in this paper, and the results confirmed the correctness of this paper's conclusions.
Figures 7 and 8 are the time-domain response of points E and F respectively; subgraph (a) is the response without forced excitation; and sub-graph (b) is the response with forced excitation with a frequency equal to the frequency of the corresponding periodic solutions.By comparison, it was found that when there was no forced excitation, the system responses corresponding to the E and F parameter points were periodic, but became divergent when the forced excitation was loaded.This phenomenon proved that when the forced excitation frequency was close to the periodic solutions frequency of the stability boundaries, those points on the stability boundaries became unstable.This result is consistent with the conclusion that the forced excitation frequency did not meet the definition domain in the analytical deduced process.4), marked with the positions of randomly selected points A-F.
When the system parameters change, the change of the stability boundaries near δ 0 = 1 and δ 0 = 4 are more obvious than the change near δ 0 = 0, and the stability of the parameter points near the bottom of the stability boundaries is more sensitive and easy to change.Therefore, some typical points for numerical simulation were selected in the bottom area of the stability boundaries and near the boundary lines of δ 0 = 1 or δ 0 = 4, The positions of these points are shown in Figure 2, and they are marked as points A, B, C, and D. In order to verify the relationship between the periodic solutions on the stability boundaries and the forced excitation, two points E and F on the boundary lines near δ 0 = 1 or δ 0 = 4 were selected to simulate randomly, as shown in Figure 2.
The responses of each typical point were compared and verified by the numerical method.The system parameters were selected as ξ = 0.005, K = 0.005, p = 0.5, and w = 1.Figures 3-6 show the time-domain response of points A, B, C, and D respectively, in which sub-graph (a) is the response without forced excitation and sub-graph (b) is the response with forced excitation.It can be seen from the comparison diagrams in Figures 3-6 that whether the point is in the stable region or the unstable region, the forced excitation did not change the stability of each point, but changed the speed of system convergence and divergence.This is consistent with the analytical results obtained in this paper, and the results confirmed the correctness of this paper's conclusions.Figures 7 and 8 are the time-domain response of points E and F respectively; subgraph (a) is the response without forced excitation; and sub-graph (b) is the response with forced excitation with a frequency equal to the frequency of the corresponding periodic solutions.By comparison, it was found that when there was no forced excitation, the system responses corresponding to the E and F parameter points were periodic, but became divergent when the forced excitation was loaded.This phenomenon proved that when the forced excitation frequency was close to the periodic solutions frequency of the stability boundaries, those points on the stability boundaries became unstable.This result is consistent with the conclusion that the forced excitation frequency did not meet the definition domain in the analytical deduced process.

Effects of Damping and Fractional-Order Parameters on the Stability Boundaries
In this section, the effect of system parameters on stability boundaries is analyzed.The basic parameters were selected as 0.005 K = , 0.5 p = , and 0 F = ;  was selected as 0, 0.03, and 0.05, respectively The stability boundaries are shown in Figure 9.It can be seen that the increase in the damping  made the stability boundaries move upward and the range of the stable region expanded gradually.

Effects of Damping and Fractional-Order Parameters on the Stability Boundaries
In this section, the effect of system parameters on stability boundaries is analyzed.The basic parameters were selected as K = 0.005, p = 0.5, and F = 0; ξ was selected as 0, 0.03, and 0.05, respectively The stability boundaries are shown in Figure 9.It can be seen that the increase in the damping ξ made the stability boundaries move upward and the range of the stable region expanded gradually.

Effects of Damping and Fractional-Order Parameters on the Stability Boundaries
In this section, the effect of system parameters on stability boundaries is analyzed.The basic parameters were selected as 0.005 K = , 0.5 p = , and 0 F = ;  was selected as 0, 0.03, and 0.05, respectively The stability boundaries are shown in Figure 9.It can be seen that the increase in the damping  made the stability boundaries move upward and the range of the stable region expanded gradually.
When the damping and the order of the fractional-order derivative were selected as  = 0.005, p= 0.5 and the coefficients of the fractional-order derivative K were selected as 0, 0.005, and 0.01, respectively, the stability boundaries near 0 4  = were as shown in Figure 10.As the coefficient K increased, the stability boundaries moved to the left and gradually moved upward.The main reason for this phenomenon was that the equivalent stiffness and the equivalent damping of the fractional-order derivative were both gradually increasing.In addition, when the fractional-order differential coefficient K was selected as 0.005 and the fractional-order differential order p was selected as 0, 0.5, and 0.8, respectively.The stability boundaries were as shown in Figure 11.It can also be seen that with the increase in p, due to the equivalent stiffness decreasing and the equivalent damping increasing, the stability boundaries moved upward and the position obviously shifted to the right.In short, we found that the fractional-order derivative not only affects the position of the stability boundaries but also changes the range of the stability region...
When the damping and the order of the fractional-order derivative were selected as ξ = 0.005, p = 0.5 and the coefficients of the fractional-order derivative K were selected as 0, 0.005, and 0.01, respectively, the stability boundaries near δ 0 = 4 were as shown in Figure 10.As the coefficient K increased, the stability boundaries moved to the left and gradually moved upward.The main reason for this phenomenon was that the equivalent stiffness and the equivalent damping of the fractional-order derivative were both gradually increasing.In addition, when the fractional-order differential coefficient K was selected as 0.005 and the fractional-order differential order p was selected as 0, 0.5, and 0.8, respectively.The stability boundaries were as shown in Figure 11.It can also be seen that with the increase in p, due to the equivalent stiffness decreasing and the equivalent damping increasing, the stability boundaries moved upward and the position obviously shifted to the right.In short, we found that the fractional-order derivative not only affects the position of the stability boundaries but also changes the range of the stability region.

Periodic Solutions in the Mathieu Equation
The Mathieu equation is a special case of Hill's equation, and Hill's equation can be expressed as

Periodic Solutions in the Mathieu Equation
The Mathieu equation is a special case of Hill's equation, and Hill's equation can be expressed as

Stability Analysis of Periodic Solutions 4.1. Periodic Solutions in the Mathieu Equation
The Mathieu equation is a special case of Hill's equation, and Hill's equation can be expressed as where α and β are two independent parameters, √ α is the natural frequency of the system, β is the parametric excitation coefficient, f (t) = f (t + T), and T is the least period of f (t).
According to the Floquet theory, the k'T-periodic solutions of Equation ( 47) can be obtained when λ k ′ = 1.To find the value of λ for the k'T-periodic solutions, we used its polar form, i.e.: x(t According to r k ′ = 1, e jk ′ θ = cos(k ′ θ) + j sin(k ′ θ) = 1, the angle conditions for k'Tperiodic solutions is as follows: From the angle conditions [24], we determined that the unit circle shows the eigenvalues' positions of k'T-periodic solutions, as shown in Figure 12.Consequently, the eigenvalues  corresponding to any k' can be obtained: From the unit circle diagram, the corresponding eigenvalues of different periods can be clearly defined.Therefore, when the value of (  ,  ) is fixed and the Floquet multiplier satisfies Equation (49), the k'T-periodic solution can be obtained.As k' is increased, the angle condition  gradually approached zero, and the k'T-periodic solutions were grad- ually closer to T -periodic solutions, thus affecting the stability of the stability boundary.Next, we considered the homogeneous Mathieu equation: In Equation ( 51), the least period T  = , and the ( ) cos2 f t t = satisfies the following equation: The stability-chart of Equation ( 51) is shown in Figure 13.There, exit  -period solu- tions or 2 -period solutions are shown on the stability boundaries.In addition, there Consequently, the eigenvalues λ corresponding to any k' can be obtained: From the unit circle diagram, the corresponding eigenvalues of different periods can be clearly defined.Therefore, when the value of (δ, ε) is fixed and the Floquet multiplier satisfies Equation (49), the k'T-periodic solution can be obtained.As k' is increased, the angle condition θ gradually approached zero, and the k'T-periodic solutions were gradually closer to T-periodic solutions, thus affecting the stability of the stability boundary.
Next, we considered the homogeneous Mathieu equation: ..
In Equation ( 51), the least period T = π, and the f (t) = cos 2t satisfies the following equation: The stability-chart of Equation ( 51) is shown in Figure 13.There, exit π-period solutions or 2π-period solutions are shown on the stability boundaries.In addition, there were stable periodic solutions of k'π-period distributed in the stable region (i.e., the blank region).Then, we analyzed the non-homogeneous Mathieu equation, as follows: where the angular frequency of the forced excitation is wi, and the least period of the forced excitation satisfies When there was no forced excitation, the k'T-periodic solutions remained stable in the stable region.However, when the forced excitation was loaded and the least period of the forced excitation was the same as the least period of the k'T-periodic solutions, those parameter points corresponding to the k'T-periodic solutions in δ -ε plane became di- vergent.Meanwhile, the angular frequency and the coefficient k' satisfies ' 2/ i k w = .

Numerical Simulation
In this section, the stability of k'T-periodic solutions is analyzed through numerical simulation to verify the correctness of the theoretical research.
When 3 δ = , 1 ε = , the system response of the Mathieu equation without the forced excitation is as shown in Figure 14.It is a steady periodic motion.The spectrum diagram of this response is shown in Figure 15.We found that the frequency of this periodic motion was 0.053 Hz, meaning the angular frequency was 6π .We selected , and this tim, the forced excitation frequency was 6π .The response of the non-homogeneous Mathieu equation is shown in Figure 16.It became divergent.From the results of these three figures, it was obvious that when the points with k'T-periodic solutions in the stable region were affected by the same frequency of forced excitation, their stability changed from stable to divergent, due to resonance.x + (δ + 2ε cos(2t))x = 0.
Then, we analyzed the non-homogeneous Mathieu equation, as follows: ..
where the angular frequency of the forced excitation is w i, and the least period of the forced excitation satisfies T i = 2π/w i .When there was no forced excitation, the k'T-periodic solutions remained stable in the stable region.However, when the forced excitation was loaded and the least period of the forced excitation was the same as the least period of the k'T-periodic solutions, those parameter points corresponding to the k'T-periodic solutions in δ-ε plane became divergent.Meanwhile, the angular frequency and the coefficient k' satisfies k ′ = 2/w i .

Numerical Simulation
In this section, the stability of k'T-periodic solutions is analyzed through numerical simulation to verify the correctness of the theoretical research.
When δ = 3, ε = 1, the system response of the Mathieu equation without the forced excitation is as shown in Figure 14.It is a steady periodic motion.The spectrum diagram of this response is shown in Figure 15.We found that the frequency of this periodic motion was 0.053 Hz, meaning the angular frequency was 6π.We selected w i = 1/3, and this tim, the forced excitation frequency was 6π.The response of the non-homogeneous Mathieu equation is shown in Figure 16.It became divergent.From the results of these three figures, it was obvious that when the points with k'T-periodic solutions in the stable region were affected by the same frequency of forced excitation, their stability changed from stable to divergent, due to resonance.
To analyze the distribution law of such unstable points more carefully, five kinds of forced excitation with different angular frequencies (w i = 1/n , n = (3, 4 . . .7)) were selected.The simulated results are shown in Figure 17.As can be seen from the figure, a large number of unstable points appeared in the stable region, due to the influence of forced excitation, and these points formed multiple linear regions with extremely thin widths (hereinafter referred to as the resonance line).The positions of these resonance lines were affected by forced-excitation frequency, but these resonance lines had a similar trend.In addition, we found that as the value of n increased, the position of the k'T-resonance lines gradually approached the stability boundaries.These results were consistent with the theoretical results in Section 4.1.
was 0.053 Hz, meaning the angular frequency was 6 .We selected , and this tim, the forced excitation frequency was 6 .The response of the non-homogeneous Mathieu equation is shown in Figure 16.It became divergent.From the results of these three figures, it was obvious that when the points with k'T-periodic solutions in the stable region were affected by the same frequency of forced excitation, their stability changed from stable to divergent, due to resonance.
To analyze the distribution law of such unstable points more carefully, five kinds of forced excitation with different angular frequencies ( ) were selected.The simulated results are shown in Figure 17.As can be seen from the figure, a large number of unstable points appeared in the stable region, due to the influence of forced excitation, and these points formed multiple linear regions with extremely thin widths (hereinafter referred to as the resonance line).The positions of these resonance lines were affected by forced-excitation frequency, but these resonance lines had a similar trend.In addition, we found that as the value of n increased, the position of the k'T-resonance lines gradually approached the stability boundaries.These results were consistent with the theoretical results in Section 4.1.) were selected.The simulated results are shown in Figure 17.As can be seen from the figure, a large number of unstable points appeared in the stable region, due to the influence of forced excitation, and these points formed multiple linear regions with extremely thin widths (hereinafter referred to as the resonance line).The positions of these resonance lines were affected by forced-excitation frequency, but these resonance lines had a similar trend.In addition, we found that as the value of n increased, the position of the k'T-resonance lines gradually approached the stability boundaries.These results were consistent with the theoretical results in Section 4.1.To analyze the distribution law of such unstable points more carefully, five kinds of forced excitation with different angular frequencies ( i w n n = =) were selected.The simulated results are shown in Figure 17.As can be seen from the figure, a large number of unstable points appeared in the stable region, due to the influence of forced excitation, and these points formed multiple linear regions with extremely thin widths (hereinafter referred to as the resonance line).The positions of these resonance lines were affected by forced-excitation frequency, but these resonance lines had a similar trend.In addition, we found that as the value of n increased, the position of the k'T-resonance lines gradually approached the stability boundaries.These results were consistent with the theoretical results in Section 4.1.

Effects of Damping and Fractional-Order Parameters on Resonance Lines
To analyze more intuitively the influence of the system parameters on the resonance line in Equation (4), we selected the parameter points of the 6T-resonance line in δ = 2.7~3.1 and ε = 0.5~0.8 for local amplification analysis.The system parameters were selected as K= 0.075, p = 0.5, ξ = 0, w i = 1/3, and F = 1; the step size of ε was 0.002, and the step size of δ was 0.001.The local chart of the 6T-resonance line of the Mathieu equation with forced excitation is shown in Figure 18.The red points represent unstable points, and the blue points represent stable points.It can be seen that the resonance line was in the shape of a long strip in the case of local magnification, and its width gradually increased as the value ε became larger.

Effects of Damping and Fractional-Order Parameters on Resonance Lines
To analyze more intuitively the influence of the system parameters on the resonance line in Equation (4), we selected the parameter points of the 6T-resonance line in  = 2.7~3.1 and  = 0.5~0.8 for local amplification analysis.The system parameters were se- lected as K= 0.075, p = 0.5, 0  = , wi = 1/3, and F = 1; the step size of  was 0.002, and the step size of  was 0.001.The local chart of the 6T-resonance line of the Mathieu equation with forced excitation is shown in Figure 18.The red points represent unstable points, and the blue points represent stable points.It can be seen that the resonance line was in the shape of a long strip in the case of local magnification, and its width gradually increased as the value  became larger.
Next, the effects of linear damping and fractional-order parameters on resonance lines were analyzed.The basic parameters of Equation (4) were selected as F = 1, wi = 1, K = 0.005, and p = 0.5, and the linear damping was selected as 0, 0.001, and 0.005, respectively.The 6T-resonance line is shown in Figure 19.As can be seen from Figure 19, the resonance line gradually shrank upward and narrowed in width as the damping increased.Therefore, the damping changed the position of the resonance line, and the increase in damping not only reduced the influence of forced excitation, but also enlarged the range of the system stability region...
Next, the effects of linear damping and fractional-order parameters on resonance lines were analyzed.The basic parameters of Equation (4) were selected as F = 1, w i = 1, K = 0.005, and p = 0.5, and the linear damping was selected as 0, 0.001, and 0.005, respectively.The 6T-resonance line is shown in Figure 19.As can be seen from Figure 19, the resonance line gradually shrank upward and narrowed in width as the damping increased.Therefore, the damping changed the position of the resonance line, and the increase in damping not only reduced the influence of forced excitation, but also enlarged the range of the system stability region.

Effects of Damping and Fractional-Order Parameters on Resonance Lines
To analyze more intuitively the influence of the system parameters on the resonance line in Equation (4), we selected the parameter points of the 6T-resonance line in  = 2.7~3.1 and  = 0.5~0.8 for local amplification analysis.The system parameters were se- lected as K= 0.075, p = 0.5, 0  = , wi = 1/3, and F = 1; the step size of  was 0.002, and the step size of  was 0.001.The local chart of the 6T-resonance line of the Mathieu equation with forced excitation is shown in Figure 18.The red points represent unstable points, and the blue points represent stable points.It can be seen that the resonance line was in the shape of a long strip in the case of local magnification, and its width gradually increased as the value  became larger.
Next, the effects of linear damping and fractional-order parameters on resonance lines were analyzed.The basic parameters of Equation (4) were selected as F = 1, wi = 1, K = 0.005, and p = 0.5, and the linear damping was selected as 0, 0.001, and 0.005, respectively.The 6T-resonance line is shown in Figure 19.As can be seen from Figure 19, the resonance line gradually shrank upward and narrowed in width as the damping increased.Therefore, the damping changed the position of the resonance line, and the increase in damping not only reduced the influence of forced excitation, but also enlarged the range of the system stability region.The effects of fractional parameters on resonance lines are discussed below.The basic parameters of Equation (4) were selected as ξ = 0, F = 1, w i = 1/3, and p = 0.5, and K was 0.01, 0.075, and 0.005, respectively.The 6T-resonance lines are shown in Figure 20.It can be seen that the change in the fractional-order derivative coefficient K changed the position and range of the resonance line.When K increased gradually, the position of the resonance line shifted to the left, due to the increase in the equivalent stiffness of the system caused by the fractional-order derivative term.In addition, the resonance line moved upwards and became narrower, due to the gradual increase in the equivalent damping.
resonance line shifted to the left, due to the increase in the equivalent stiffness of the sys-tem caused by the fractional-order derivative term.In addition, the resonance line moved upwards and became narrower, due to the gradual increase in the equivalent damping.
When K was 0.005 and p was 0, 0.5, and 1, the 6T-resonance lines were as shown in Figure 21.It also can be seen that as p increased, the resonance line gradually moved to the right.This was because the equivalent stiffness of the system decreased gradually, due to the fractional-derivative term.In addition, as p increased, the equivalent linear damping gradually increased, and the resonance line gradually moved upward.In sum, both the damping and the fractional-order parameters had important effects on the position and range of the resonance line, caused by forced excitation.

Conclusions
On the basis of the model of the pantograph-catenary system, the stability analysis of the non-homogeneous fractional-order Mathieu equation was concluded by the perturbation method and numerical simulation.We proved that forced excitation with a finite period does not change the convergence of the stability boundary.Only when the period of the forced excitation is large enough will the stability boundary become unstable.We verified that k'T-periodic solutions in the stable region of the non-homogeneous Mathieu equation are unstable under the influence of forced excitation.When the frequency of forced excitation is close to the frequency of the k'T-periodic solutions, unstable resonance lines appear in the stable region.Finally, by analyzing the influence of system parameters on the resonance line, we found that both the damping and fractional-order derivative term in the system affect the position and range of the resonance line region.These results When K was 0.005 and p was 0, 0.5, and 1, the 6T-resonance lines were as shown in Figure 21.It also can be seen that as p increased, the resonance line gradually moved to the right.This was because the equivalent stiffness of the system decreased gradually, due to the fractional-derivative term.In addition, as p increased, the equivalent linear damping gradually increased, and the resonance line gradually moved upward.In sum, both the damping and the fractional-order parameters had important effects on the position and range of the resonance line, caused by forced excitation.
Fractal Fract.2022, 6, x FOR PEER REVIEW 18 of 20 be seen that the change in the fractional-order derivative coefficient K changed the position and range of the resonance line.When K increased gradually, the position of the resonance line shifted to the left, due to the increase in the equivalent stiffness of the system caused by the fractional-order derivative term.In addition, the resonance line moved upwards and became narrower, due to the gradual increase in the equivalent damping.
When K was 0.005 and p was 0, 0.5, and 1, the 6T-resonance lines were as shown in Figure 21.It also can be seen that as p increased, the resonance line gradually moved to the right.This was because the equivalent stiffness of the system decreased gradually, due to the fractional-derivative term.In addition, as p increased, the equivalent linear damping gradually increased, and the resonance line gradually moved upward.In sum, both the damping and the fractional-order parameters had important effects on the position and range of the resonance line, caused by forced excitation.

Conclusions
On the basis of the model of the pantograph-catenary system, the stability analysis of the non-homogeneous fractional-order Mathieu equation was concluded by the perturbation method and numerical simulation.We proved that forced excitation with a finite period does not change the convergence of the stability boundary.Only when the period of the forced excitation is large enough will the stability boundary become unstable.We verified that k'T-periodic solutions in the stable region of the non-homogeneous Mathieu equation are unstable under the influence of forced excitation.When the frequency of forced excitation is close to the frequency of the k'T-periodic solutions, unstable resonance lines appear in the stable region.Finally, by analyzing the influence of system parameters on the resonance line, we found that both the damping and fractional-order derivative term in the system affect the position and range of the resonance line region.These results

Conclusions
On the basis of the model of the pantograph-catenary system, the stability analysis of the non-homogeneous fractional-order Mathieu equation was concluded by the perturbation method and numerical simulation.We proved that forced excitation with a finite period does not change the convergence of the stability boundary.Only when the period of the forced excitation is large enough will the stability boundary become unstable.We verified that k'T-periodic solutions in the stable region of the non-homogeneous Mathieu equation are unstable under the influence of forced excitation.When the frequency of forced excitation is close to the frequency of the k'T-periodic solutions, unstable resonance lines appear in the stable region.Finally, by analyzing the influence of system parameters on the resonance line, we found that both the damping and fractional-order derivative term in the system affect the position and range of the resonance line region.These results reveal the complex dynamic behavior of fractional-order parametric excitation systems under the influence of forced excitation, information that can be used to guide system design and optimization.

Figure 1 .
Figure 1.The fractional-order model of the pantograph-catenary system.

0 1 δ 1 δ
= or 0 4 δ = , The positions of these points are shown in Figure2, and they are marked as points A, B, C, and D. In order to verify the relationship between the periodic solutions on the stability boundaries and the forced excitation, two points E and F on the boundary lines near 0 = or 0 4 δ = were selected to simulate randomly, as shown in Figure2.

Figure 2 .
Figure 2. Stability chart of system (4), marked with the positions of randomly selected points A-F.

Figure 2 .
Figure 2. Stability chart of system (4), marked with the positions of randomly selected points A-F.

Figure 11 . 4 
Figure 11.The effects of p on stability boundaries ( 0

Figure 11 . 4 
Figure 11.The effects of p on stability boundaries ( 0

FractalFigure 12 .
Figure 12.The unit circle shows the eigenvalues' positions of k'T-periodic solutions.

Figure 12 .
Figure 12.The unit circle shows the eigenvalues' positions of k'T-periodic solutions.

Figure 16 .
Figure 16.System response of

Figure 18 .
Figure 18.Local chart of the 6T-resonance line:

Figure 19 .
Figure 19.The effects of  on the 6T-resonance line.

Figure 18 .
Figure 18.Local chart of the 6T-resonance line:

Figure 18 .
Figure 18.Local chart of the 6T-resonance line:

Figure 19 .
Figure 19.The effects of  on the 6T-resonance line.The effects of fractional parameters on resonance lines are discussed below.The basic parameters of Equation (4) were selected as 0  = , F = 1, wi = 1/3, and p = 0.5, and K was 0.01, 0.075, and 0.005, respectively.The 6T-resonance lines are shown in Figure20.It can

Figure 19 .
Figure 19.The effects of ξ on the 6T-resonance line.

Figure 20 .
Figure 20.The effects of K on the 6T-resonance line.

Figure 21 .
Figure 21.The effects of p on the 6T-resonance line.

Figure 20 .
Figure 20.The effects of K on the 6T-resonance line.

Figure 20 .
Figure 20.The effects of K on the 6T-resonance line.

Figure 21 .
Figure 21.The effects of p on the 6T-resonance line.

Figure 21 .
Figure 21.The effects of p on the 6T-resonance line.