Abstract
A dynamical analysis of a Mathieuvan der PolDuffing nonlinear system with fractionalorder derivative under combined parametric and forcing excitation is studied in this paper. The approximate analytical solution is researched for 1/2 sub harmonic resonance coupled with primary parametric resonance based on the improved averaging approach. The steadystate periodic solution including its stability condition is established. The equivalent linear stiffness coefficient (ELDC) and equivalent linear damping coefficient (ELSC) for this nonlinear fractionalorder oscillator are defined. Then, the numerical simulations are presented in three typical cases by iterative algorithms. The time history, phase portrait, FFT spectrum and Poincare maps are shown to explain the system response. Some different responses, such as quasiperiodic, multiperiodic and single periodic behaviors are observed and investigated. The results of comparisons between the numerical solutions and the approximate analytical solutions in three typical cases show the correctness of the analytical solutions. The influences of the fractionalorder parameters on the system dynamical response are researched based on the ELDC and ELSC. Through analysis, it could be found that the increase of the fractionalorder coefficient would result in the rightward and downward movements of the amplitudefrequency curves. The increase of the fractionalorder coefficient will also move the bifurcation point rightwards and will make the existing range of steadystate solution larger. It could also be found that the ELSC will become larger and ELDC smaller when the fractional order is closer to zero, so that the decrease of the fractional order would make the response amplitude larger. At last, the detailed conclusions are summarized, which is beneficial to design and control this kind of fractionalorder nonlinear system.
1. Introduction
With the development of modern computation technology in recent years, the fractional calculus has caused wide concerns in many fields, especially in mechanical, civil, biological and other engineering fields. The fundamental characteristics of fractional calculus are studied, and great progress has been made in the basic theory^{}[14]. For example, Petras [12] studied the modeling, analysis and simulation of the fractionalorder nonlinear systems, and they also focused on the modeling and control of the fractionalorder systems. Caputo et al. [3] proposed a new definition of the fractional derivative. Yang et al. [4] investigated the local fractional integral transforms and the applications problem, etc.
At present, the applications of the fractional calculus in engineering fields have been attracted much attention. Many engineering problems can be described by the fractional calculus equations to simulate some dynamical systems accurately. So, the fractional calculus on the dynamical system was important and interesting, which had been investigated by many researchers in some relevant fields [5, 6]. For example, Yang et al. [7] analyzed the response of a fractionalorder Duffing system by the harmonic balance method, and concluded that the fractionalorder damping could change the number of the steady stable states. Wang et al. [8] investigated the dynamical characteristics of a linear oscillator with fractionalorder derivative under forcing excitation, and found that the system solution had two parts in the fractionalorder linear oscillator based on the theory of the stability switches. Shen et al. [9, 10] investigated several fractionalorder oscillators by the averaging method, and proposed the fractionalorder derivative that had both the function of damping and the function of stiffness. Many mathematical theories of the fractional calculus were researched by Li et al. [11, 12]. They also compared the differences and connections between three typical definitions of fractionalorder derivatives, and presented some characteristics of the Caputo fractional calculus. Chen et al. [13] studied the dynamical response of the fractionalorder stochastic oscillator based on the stochastic averaging approach, and they mainly analyzed the influences of the fractionalorder derivative on the dynamics and the control of the stochastic system. Leung et al. [14] raised the residue harmonic balance approach to solve the fractionalorder van der Pol (VDP) system and the higherorder approximate solution could be deduced according to this method. You et al. [15] studied the optimal control of a fractionalorder vehicle suspension system, and the viscoelastic characteristics of the vehicle suspension system were investigated. Shen et al. [16] studied the fractionalorder PID controller on a nonsmooth oscillator and found the fractionalorder PID controller was better than the traditional PID controller. Wang et al. [17] investigated the chaos behavior of the fractionalorder Lorenz systems by the active feedback approach. Then the fractional order of 0.98 was implemented and the result demonstrated the validity of this method. Yang et al. [18, 19] studied several fractionalorder operators in the sense of Caputo definition, and discussed their characteristics based on the Laplace transform and Fourier transform. Those results were very efficient and helpful to investigate some complex phenomenon. Li et al. [20] investigated the stochastic bifurcations of DuffingVDP system under colored noise by the stochastic averaging method, and found that the change of fractional order, noise intensity and correlation time would lead to bifurcation of this system. Niu et al. [21] studied the Duffing oscillator with delayed fractionalorder PID control based on the KBM asymptotic method, and analyzed the influences of time delay. Giresse et al. [22] studied the synchronization of chaotic MathieuVDP and chaotic DuffingVDP systems with the fractionalorder derivative, and found fractionalorder derivative would lead to faster synchronization than the integerorder derivative.
There are many nonlinear dynamics phenomena in engineering fields, such as nonsmoothness, large deformation, backlash, and parametric excitation, etc., which will deteriorate the performance of engineering system [23, 24]. For example, Li et al. [2526] found the bursting phenomenon in a mechanical system and analyzed the slowfast effect. Belhaq et al. [27] analyzed the resonances of VDPMathieuDuffing oscillator in the frequencylocking area of 2:1 and 1:1, and found that the fast harmonic excitation could change the spring behavior. Based on the authors’ knowledge, the published papers were mainly focused on the dynamical responses of the fractionalorder system with forcing excitation and the integerorder nonlinear system. A few researchers studied the fractionalorder complex nonlinear system. For example, Wen et al. [28, 29] investigated the dynamical characteristics of MathieuDuffing oscillator with the fractionalorder delayed feedback, and found that the fractionalorder delayed feedback had the characteristics of both delayed velocity feedback and displacement feedback. Xu et al. [30] analyzed the dynamical responses of fractionalorder DuffingRayleigh oscillator with the Gauss white noise excitation. But the dynamical analysis of the complex fractionalorder nonlinear system with combined parametric and forcing excitation hardly has been investigated.
In this paper, the dynamical responses of Mathieuvan der PolDuffing nonlinear system with the fractionalorder derivative for 1/2 subharmonics resonance coupled with the primary parametric resonance are studied. The paper is arranged as follows: In Section 2, the approximate analytical solution of the fractionalorder Mathieuvan der PolDuffing nonlinear system is obtained. Then the ELDC and ELSC are also defined in this section. In Section 3, the amplitudefrequency equation and the stability conditions of the steadystate nonzero periodic solution are presented. Then in Section 4 the numerical solutions in three typical cases are obtained. The comparisons between the analytical results and the numerical solutions in three cases are presented to prove the correctness of the approximate analytical solution. The effects of the fractionalorder parameters are studied in Section 5. At last, the detailed conclusions are summarized and analyzed.
2. Firstorder approximate analytical solution
The considered fractionalorder nonlinear system with 1/2 subharmonics resonance coupled with primary parametric resonance is:
$+{\alpha}_{2}{x}^{3}\left(t\right)+{K}_{2}{D}^{p}\left[x\left(t\right)\right]=F\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega}t\right),$
where $m$, ${K}_{1}$, ${\alpha}_{1}$, ${\alpha}_{2}$, $F$, $\mathrm{\Omega}$ and $2B\mathrm{c}\mathrm{o}\mathrm{s}2\omega t$ are the system mass, linear stiffness coefficient, nonlinear damping coefficient, nonlinear stiffness coefficient, forcing excitation coefficient, forcing excitation frequency and the periodic timevarying stiffing coefficient respectively. ${D}^{p}\left[x\right(t\left)\right]$ is the $p$order derivative of $x\left(t\right)$ to $t$ with the fractional coefficient ${K}_{2}$ (${K}_{2}>0$) and fractional order $p$ ($0\le p\le 1$). Here the Caputo definition is used below:
where $\mathrm{\Gamma}\left(x\right)$ is the Gamma function which satisfies $\mathrm{\Gamma}\left(x+1\right)=x\mathrm{\Gamma}\left(x\right)$.
The following transformations are used:
Eq. (1) becomes:
$+\epsilon {\alpha}_{2}^{\text{'}}{x}^{3}\left(t\right)+\epsilon {k}_{2}{D}^{p}\left[x\left(t\right)\right]=\epsilon f\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega}t\right).$
Here, we focus on the 1/2 sub harmonic resonance coupled with the primary parametric resonance. It means the parametric frequency is close to the natural frequency and the forcing excitation frequency is close to the double of the parametric frequency, i.e. $\mathrm{\Omega}=2\omega $, $\omega \approx {\omega}_{0}$. The formula ${\omega}^{2}={\omega}_{0}^{2}+\epsilon \sigma $ is introduced to denote the approximate degree, where $\sigma $ is the detuning factor, so that Eq. (4) becomes:
Supposing the solution of Eq. (5) is:
where $\phi =\omega t+\theta $. Based on the improved averaging approach [9, 10], we could integrate the righthand part of Eq. (5) within the time interval [0, $T$] and then average it. The approximate solution of amplitude $a$ and phase $\theta $ are:
$a\dot{\theta}=\frac{1}{T\omega}{\int}_{0}^{T}\left[{P}_{1}\left(a,\theta \right)+{P}_{2}\left(a,\theta ,\tau \right)\right]\mathrm{c}\mathrm{o}\mathrm{s}\phi d\phi ,$
where:
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}2\beta \mathrm{c}\mathrm{o}\mathrm{s}\left[2\left(\phi \theta \right)\right]a\mathrm{c}\mathrm{o}\mathrm{s}\phi {\alpha}_{2}^{\text{'}}{a}^{3}\mathrm{c}\mathrm{o}\mathrm{s}\phi \},$
${P}_{2}\left(a,\theta \right)=\epsilon {k}_{2}{D}^{p}\left[a\mathrm{c}\mathrm{o}\mathrm{s}\phi \right].$
When Eq. (7) is expanded, it will yield:
The first part of Eq. (8a) and Eq. (8b) are the periodic function with $T=2\pi $. Integrating the first part of Eq. (8), one could get:
$a{\dot{\theta}}_{1}=\frac{1}{2\pi \omega}{\int}_{0}^{2\pi}{P}_{1}\left(a,\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\phi d\phi =\frac{\epsilon a\beta \mathrm{c}\mathrm{o}\mathrm{s}2\theta}{2\omega}\frac{\epsilon \sigma a}{2\omega}+\frac{3\epsilon {\alpha}_{2}^{\text{'}}{a}^{3}}{8\omega}.$
The second part of Eq. (8a) and Eq. (8b) are the aperiodic functions, and the time terminal $T$ is selected as $\infty $:
$a{\dot{\theta}}_{2}=\underset{T\to \infty}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{1}{T\omega}{\int}_{0}^{T}{P}_{2}\left(a,\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\phi d\phi =\underset{T\to \infty}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\epsilon {k}_{2}}{T\omega}{\int}_{0}^{T}{D}^{p}\left[a\mathrm{cos}\left(\phi \right)\right]\mathrm{c}\mathrm{o}\mathrm{s}\phi d\phi .$
Considering Eq. (2), one can obtain:
Supposing $s=tu$, $ds=du$, and when they are substituted into Eq. (11), it becomes:
$=\frac{\epsilon a{k}_{2}}{\mathrm{\Gamma}(1p)}\underset{T\to \infty}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{1}{T}{\int}_{0}^{T}\left\{\left[{\int}_{0}^{t}\frac{\mathrm{c}\mathrm{o}\mathrm{s}\omega s}{{s}^{p}}ds\right]\mathrm{s}\mathrm{i}\mathrm{n}(\omega t+\theta )\mathrm{s}\mathrm{i}\mathrm{n}(\omega t+\theta )\right\}dt$
$+\frac{\epsilon a{k}_{2}}{\mathrm{\Gamma}(1p)}\underset{T\to \infty}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{1}{T}{\int}_{0}^{T}\left\{\left[{\int}_{0}^{t}\frac{\mathrm{s}\mathrm{i}\mathrm{n}\omega s}{{s}^{p}}ds\right]\mathrm{c}\mathrm{o}\mathrm{s}(\omega t+\theta )\mathrm{s}\mathrm{i}\mathrm{n}(\omega t+\theta )\right\}dt.$
Two important formulae in Refs [910] are introduced:
$\underset{T\to \infty}{\mathrm{l}\mathrm{i}\mathrm{m}}{\int}_{0}^{T}\frac{\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)}{{t}^{p}}dt={\omega}^{p1}\mathrm{\Gamma}\left(1p\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{p\pi}{2}\right).$
In Eq. (12), defining the first part as ${\dot{a}}_{21}$ and the second part as ${\dot{a}}_{22}$, integrating it by parts respectively and considering Eq. (13), one can get:
$+\frac{\epsilon a{k}_{2}}{4\omega \mathrm{\Gamma}(1p)}\underset{T\to \infty}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{1}{T}{\int}_{0}^{T}\left\{[2\omega t\mathrm{s}\mathrm{i}\mathrm{n}(2\omega t+2\theta \left)\right]\frac{\mathrm{c}\mathrm{o}\mathrm{s}\omega t}{{t}^{p}}\right\}dt=\frac{\epsilon a{k}_{2}{\omega}^{p1}}{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{p\pi}{2}\right).$
Similarly, one can find the integral of the second part ${\dot{a}}_{22}$ in Eq. (12) approach 0 when $T\to \infty $, so that:
When the same method is used, it will yield:
Combining Eq. (9) and Eq. (15), one can obtain:
$a\dot{\theta}=\frac{\epsilon a\beta \mathrm{c}\mathrm{o}\mathrm{s}2\theta}{2\omega}\frac{\epsilon \sigma a}{2\omega}+\frac{3\epsilon {\alpha}_{2}^{\text{'}}{a}^{3}}{8\omega}+\frac{\epsilon a{k}_{1}{\omega}^{p1}}{2}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{p\pi}{2}\omega \tau \right).$
When the original parameters are substituted into Eq. (16), it will yield:
$a\dot{\theta}=\frac{Ba\mathrm{c}\mathrm{o}\mathrm{s}2\theta}{2m\omega}+\frac{3{\alpha}_{2}{a}^{3}}{8m\omega}\frac{\omega a}{2}+\frac{a}{2m\omega}{K}_{e}\left(p\right),$
where:
${K}_{e}\left(p\right)={K}_{1}+{K}_{2}{\omega}^{p}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{p\pi}{2}\right),$
are defined as the equivalent linear damping coefficient (ELDC) and the equivalent linear stiffness coefficient (ELSC) respectively.
Based on Eq. (17) and Eq. (18), it could be found that the fractional coefficient ${K}_{2}$ and the fractional order $p$ are all significant in system dynamics. The fractional coefficient ${K}_{2}$ is linear with the ELDC and ELSC in this system. The fractional order $p$ affects this system dynamical response with the ELDC and ELSC in the form of trigonometric function. When $p\ne 0$ and ${K}_{2}\ne 0$, the fractionalorder derivative will possess both the function of stiffness and the function of damping simultaneously.
3. Steadystate periodic solutions and stability conditions
Letting $\dot{a}=0$ and $\dot{\theta}=0$, Eq. (17) becomes:
$\frac{B\stackrel{}{a}\mathrm{c}\mathrm{o}\mathrm{s}2\stackrel{}{\theta}}{2m\omega}+\frac{3{\alpha}_{2}{\stackrel{}{a}}^{3}}{8m\omega}\frac{\omega \stackrel{}{a}}{2}+\frac{\stackrel{}{a}}{2m\omega}{K}_{e}\left(p\right)=0.$
When $\stackrel{}{\theta}$ is eliminated from Eq. (19), the amplitudefrequency equation is:
That is:
Letting $u=a\mathrm{c}\mathrm{o}\mathrm{s}\theta $, $v=a\mathrm{s}\mathrm{i}\mathrm{n}\theta $, and when they are substituted into Eq. (17), it becomes:
The Jacobian matrix for Eq. (21) is:
The Jacobian matrix of the nonzero periodic solution is:
Letting:
and using Eq. (19), the characteristic equation can be get as follows:
From the analysis of Eq. (22), it can be found that the stability conditions for the nonzero periodic solution are:
$4({C}_{1}^{2}+{C}_{3}^{2}{C}_{2}^{2})+4({C}_{1}{C}_{4}+{C}_{3}{C}_{5}){a}^{2}<0.$
According to Eq. (20b), one could get the following results.
When $\left[{K}_{e}\right(p)m{\omega}^{2}{]}^{2}+[{C}_{e}\left(p\right)\omega {]}^{2}<{B}^{2}$, Eq. (20b) has the zero solution and one nonzero periodic solution as follows:
${a}_{2}={\left\{\frac{\begin{array}{c}4\left\{{\omega}^{2}{\alpha}_{1}{C}_{e}^{}\left(p\right)+3{\alpha}_{2}\left[{K}_{e}\left(p\right)m{\omega}^{2}\right]\right\}\\ +{\left[[{{\alpha}_{1}}^{2}{\omega}^{2}+9{{\alpha}_{2}}^{2}]{B}^{2}\left\{{\alpha}_{1}\omega \right[{K}_{e}\left(p\right)m{\omega}^{2}]3{\alpha}_{2}\omega {C}_{e}(p){\}}^{2}\right]}^{1/2}\end{array}}{{{\alpha}_{1}}^{2}{\omega}^{2}+9{{\alpha}_{2}}^{2}}\right\}}^{1/2}.$
When:
and:
Eq. (20b) has the zero solution and two nonzero periodic solutions as follows:
According to the stability conditions of Eq. (23), it can be found that the nonzero solution given by Eq. (24) is stable and there is only one nonzero solution for Eq. (1) at this moment. When Eq. (1) has two nonzero solutions, one nonzero solution given in Eq. (25a) is stable and the other one given in Eq. (25b) is unstable.
4. Numerical simulation
Numerical simulation solutions of Eq. (1) are presented in the following part. The numerical iteration methods introduced in [13] are adopted and the numerical formula is:
where ${C}_{j}^{p}$ is the fractional binomial coefficient which satisfies the iterative relationship of Eq. (27), ${t}_{l}=lh$ is the time sample points and $h$ is the sample step.
Based on Eq. (26) and Eq. (27), the numerical iterative algorithm of Eq. (1) could be obtained:
${Z}_{2}\left({t}_{l}\right)=\left\{F\mathrm{c}\mathrm{o}\mathrm{s}\right(2\omega {t}_{l}){\alpha}_{1}({Z}_{1}({t}_{l1}{)}^{2}1){Z}_{2}\left({t}_{l1}\right){\alpha}_{2}{Z}_{1}({t}_{l1}{)}^{3}$
$[k+2b\mathrm{c}\mathrm{o}\mathrm{s}(2\omega {t}_{l}\left)\right]{Z}_{1}\left({t}_{l1}\right){K}_{1}{Z}_{3}\left({t}_{l1}\right)\}h{\sum}_{j=1}^{l}{C}_{j}^{1}{Z}_{2}\left({t}_{lj}\right),$
${Z}_{3}\left({t}_{l}\right)={Z}_{2}\left({t}_{l1}\right){h}^{1p}{\sum}_{j=1}^{l}{C}_{j}^{1p}{Z}_{3}\left({t}_{lj}\right),$
where ${Z}_{1}$ is the displacement, ${Z}_{2}$ is the velocity, and ${Z}_{3}$ is the fractionalorder derivative. Here the sample step $h$ is selected as 0.001. The temporary response is omitted, and the maximum steady amplitude of the posterior steadystate response is adopted. Numerical results and subsequent analysis of Eq. (1) are presented in three typical topology cases.
The system parameters of the first illustrative example are selected as: $m=\text{1}$, ${K}_{1}=\text{1}$, ${\alpha}_{1}=\text{0.1}$, ${\alpha}_{2}=\text{0.1}$, $B=\text{0.4}$, ${K}_{2}=\text{0.1}$, $p=\text{0.5}$ and $F=\text{0.5}$. The numerical results are shown in Fig. 1. As the frequency is swept forward, the amplitudefrequency curve based on the numerical iterative algorithm is shown in Fig. 1 by small triangles and circles. From the observation of Fig. 1, it can be found that in the amplitudefrequency curve there are three sudden change points and four different response regions. The following parts will simulate the dynamical motions in different regions. In order to analyze the system motion, the time history figure only shows a partial steady state process in Figs. 27.
Fig. 1Amplitudefrequency curves: (m=1, K1=1, α1=0.1, α2=0.1, B=0.4, K2=0.1, p=0.5, F=0.5)
(1) When $\omega =$ 0.55, the time history, phase portrait, FFT spectrum and Poincare maps are as shown in Fig. 2 respectively. It can be found that there is only one main frequency in the FFT spectrum of Fig. 2(c), and there is only one point in Poincare maps of Fig. 2(d), so that the system response is simple singleperiodic motion in this case.
(2) When $\omega =$0.75, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 3 respectively. It can be found that there are three irreducible frequency components in the FFT spectrum of Fig. 3(c), and there is a closed loop in Poincare maps of Fig. 3(d), so that the system response is quasiperiodic motion in this case.
(3) When $\omega =$0.8, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 4 respectively. There are several frequency components in the FFT spectrum of Fig. 4(c) and several points in Poincare maps of Fig. 4(d). It can be found that the system response is multiperiodic motion in this case which is in the midst of quasiperiodic motion and singleperiodic motion.
(4) When $\omega =$1.0, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 5 respectively. The motion characteristics are similar to those shown in Fig. 2, so that the system response also is singleperiodic motion in this case.
(5) When $\omega =$1.5, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 6 respectively. When $\omega =$1.55, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 7 respectively. It can be found that the motion characteristics are similar to those shown in Fig. 4, so that the system responses are all multiperiodic motion in those cases.
Fig. 2System response for ω= 0.55
a) Time history
b) Phase portrait
c) FFT spectrum
d) Poincare maps
Fig. 3System response for ω= 0.75
a) Time history
b) Phase portrait
c) FFT spectrum
d) Poincare maps
Fig. 4System response for ω= 0.8
a) Time history
b) Phase portrait
c) FFT spectrum
d) Poincare maps
Fig. 5System response for ω= 1.0
a) Time history
b) Phase portrait
c) FFT spectrum
d) Poincare maps
We can find that there are quasiperiodic behaviors and multiperiodic behaviors (QMP) in Region II and IV denoted by small triangles, and singleperiodic behaviors (SP) in Region I and III denoted by small circles in Fig. 1. With the frequency increase, the SP, QMP, SP and then QMP responses will take place in this system gradually. The SP response in Region I is the super harmonic resonance, and the SP response in Region III is the primary resonance.
The numerical results of the other two cases are shown in Fig. 8 and Fig. 9. The amplitudefrequency curves of the approximate analytical solution of Eq. (20) in three typical cases are also shown in Fig. 1, Fig. 8 and Fig. 9 respectively, where the solid line shows the stable solution and the dash line shows the unstable one. From the observation of Fig. 1, Fig. 8 and Fig. 9, one can find that the analytical solutions are all in a good agreement with the numerical ones at the regions $\omega \approx $1 in three topology cases. Here 1/2 sub harmonic resonance coupled with the primary parametric resonance is focused, so that the correctness and accuracy of the analytical solutions can be verified. Furthermore, we also find that the approximate analytical solution does not agree very well with the numerical one near the bifurcation points on the both sides in Region III of Fig. 1. In order to analyze the error reason, two points $\omega =$0.87 and $\omega =$0.875 near the bifurcation points are selected to simulate time histories in Fig. 10(b) and Fig. 10(a) respectively. From the comparison of Fig. 10, it can be found that the system response approaches the steady state after 100s when $\omega =$0.875, but the system response approaches the steady state after 500 s when $\omega =$0.87. Accordingly, it could be concluded that the closer $\omega $ is to the bifurcation point, the longer the system needs the convergence time. Because of the calculation and time troubles, the fittings on both sides of the amplitudefrequency curve are not realized.
Fig. 6System response for ω= 1.5
a) Time history
b) Phase portrait
c) FFT spectrum
d) Poincare maps
Fig. 7System response for ω= 1.55
a) Time history
b) Phase portrait
c) FFT spectrum
d) Poincare maps
5. Influences of fractionalorder derivative parameters
In order to optimize the system design and to control the system response through choosing the fractionalorder derivative parameters, the influences of the fractionalorder derivative parameters on the system dynamical responses are investigated. When ${K}_{2}$ is changed, the system amplitudefrequency curves of ${\alpha}_{1}=\text{0.4}$ and ${\alpha}_{1}=\text{0.8}$ are given in Fig. 11(a, b) respectively. The system basic parameters are selected as $m=\text{1}$, ${K}_{1}=\text{1}$, ${\alpha}_{2}=\text{0.1}$, $B=\text{0.4}$, $F=\text{0.5}$ and $p=\text{0.5}$.
Fig. 8Comparison between analytical results and numerical solutions: (m=1, K1=1, α1=0.42, α2=0.01, B=0.4, K2=0.1, p=0.5, F=0.5)
Fig. 9Comparison between analytical results and numerical solutions: (m=1, K1=1, α1=0.8, α2=0.01, B=0.4, K2=0.1, p=0.5, F=0.5)
Fig. 10System response for ω= 0.875 and 0.87: Time history
a)$\omega =$ 0.875
b)$\omega =$ 0.87
From the observation of Fig. 11, it could be found that the increase of ${K}_{2}$ would result into the rightward and downward movements of the amplitudefrequency curves. The main reason is that the ELSC will become larger with the increase of ${K}_{2}$. Besides the above characteristics of the dynamical system responses, the topological structures of the amplitudefrequency curves vary significantly too. The increase of fractionalorder coefficient ${K}_{2}$ would make the ovoid curve disappear, the open curve gradually turn up, and push the bifurcation point rightward. The proportion of steadystate solution in the amplitudefrequency curve is getting bigger and bigger with the increase of ${K}_{2}$.
Fig. 11Influences of K2
a)${\alpha}_{1}=\text{0.4}$
b)${\alpha}_{1}=\text{0.8}$
When $p$ changes from 1 to 0, the amplitudefrequency curves of the fractional coefficient ${K}_{2}=\text{0.1}$ and ${K}_{2}=\text{0.5}$ are given in Fig. 12(a, b) respectively. The basic system parameters are selected as $m=\text{1}$, ${K}_{1}=\text{1}$, ${\alpha}_{2}=\text{0.1}$, $B=\text{0.4}$, $F=\text{0.5}$ and ${\alpha}_{1}=\text{0.4}$. Through the analysis of Fig. 12, one can find that the decrease of $p$ will result into the rightward and upward movements of the amplitudefrequency curves. The open curve disappears, the ovoid curve gradually turns up, and the bifurcation point will move rightwards with the decrease of $p$. Through the analysis of Eq. (17), one can find that the reason is that the ELSC will become larger, and the ELDC will become smaller when the fractional order $p\to \text{0}$. The ELSC will become smaller, and the ELDC will become larger when the fractional order $p\to \text{1}$. Accordingly, the decrease of $p$ will make the maximum response amplitude larger.
Fig. 12Influences of the fractional order p
a)${K}_{2}=$ 0.1
b)${K}_{2}=\mathrm{}$0.5
6. Conclusions
The dynamical characteristics of a Mathieuvan der PolDuffing fractionalorder nonlinear system under 1/2 subharmonic resonance coupled with primary parametric resonance are investigated. The approximate analytical system solution is presented by the improved averaging approach. The numerical solutions are also studied by an iterative approach in three typical cases. The time history, phase portrait, FFT spectrum and Poincare maps are presented to identify the system response. There are Quasiperiodic, multiperiodic or single periodic behaviors. As the frequency is swept forward, the QMP response, SP response and QMP response will takes place in this system gradually, which is observed and investigated in three cases. Then the numerical results are compared with the analytical results. From three comparisons between the numerical and analytical results, one can find that the stable periodic solutions are in a good agreement with the numerical solutions, which verify the correctness and high accuracy of the approximate analytical solution. Finally, the influences of the fractionalorder derivative parameters on the system dynamical response are analyzed. Through the analysis of the fractionalorder derivative parameters, it can be found that the ELSC and ELDC both will become larger with the increase of fractionalorder coefficient ${K}_{2}$. The ELSC will become larger, and the ELDC will become smaller when the fractional order $p\to \text{0}$. The ELSC will become smaller, and the ELDC will become larger if the fractional order $p\to \text{1}$. A change of the fractionalorder derivative parameters would lead to the motion of the amplitudefrequency curve and to a change of the topological structure of the amplitudefrequency curve. Those conclusions are very useful for the system design and control of the system dynamical response by choosing suitable fractionalorder derivative parameters.
References

Petras I. Fractionalorder Nonlinear Systems: Modeling, Analysis and Simulation. Higher Education Press, Beijing, 2011.

Caponetto R., Dongola G., Fortuna L., Petras I. Fractional Order Systems: Modeling and Control Applications. World Scientific, New Jersey, 2010.

Caputo M., Fabrizio M. A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation and Applications, Vol. 1, Issue 2, 2015, p. 7385.

Yang X. J., Baleanu D., Srivastava H. M. Local Fractional Integral Transforms and Their Applications. Academic Press, London, 2015.

Agnieszka B., Malinowska, Delfim Torres F. M. Fractional calculus of variations for a combined Caputo derivative. Fractional Calculus and Applied Analysis, Vol. 14, Issue 4, 2011, p. 523537.

Machado J. T., Kiryakova V., Mainardi F. Recent history of fractional calculus. Communications in Nonlinear Science and Numerical Simulation, Vol. 16, Issue 3, 2011, p. 11401153.

Yang J. H., Zhu H. The response property of one kind of fractionalorder linear system excited by different periodical signals. Acta Physica Sinica, Vol. 62, Issue 2, 2013, p. 24501, (in Chinese).

Wang Z. H., Hu H. Y. Stability of a linear oscillator with damping force of fractionalorder derivative. Science in China G: Physics, Mechanics and Astronomy, Vol. 39, Issue 10, 2009, p. 14951502, (in Chinese).

Shen Y. J., Wei P., Yang S. P. Primary resonance of fractionalorder van der Pol oscillator. Nonlinear Dynamics, Vol. 77, Issue 4, 2014, p. 16291642.

Shen Y. J., Yang S. P., Sui C. Y. Analysis on limit cycle of fractionalorder van der Pol oscillator. Chaos, Solitons and Fractals, Vol. 67, 2014, p. 94102.

Li C. P., Deng W. H. Remarks on fractional derivatives. Applied Mathematics and Computation, Vol. 187, Issue 2, 2007, p. 777784.

Deng W. H., Li C. P. The evolution of chaotic dynamics for fractional unified system. Physics Letters A, Vol. 372, Issue 4, 2008, p. 401407.

Chen L. C., Hu F., Zhu W. Q. Stochastic dynamics and fractional optimal control of quasi integrable Hamiltonian systems with fractional derivative damping. Fractional Calculus and Applied Analysis, Vol. 16, Issue 1, 2013, p. 189225.

Leung A. Y. T., Guo Z. J., Yang H. X. The residue harmonic balance for fractional order van der Pol like oscillators. Journal of Sound and Vibration, Vol. 331, Issue 5, 2012, p. 11151126.

You H., Shen Y. J., Yang S. P. Optimal design for fractionalorder active isolation system. Advances in Mechanical Engineering, Vol. 7, Issue 12, 2015, p. 111.

Shen Y. J., Niu J. C., Yang S. P., Li S. J. Primary resonance of dryfriction oscillator with fractionalorder PID Controller of velocity feedback. ASME Journal of Computational and Nonlinear Dynamics, Vol. 11, Issue 5, 2016, p. 51027.

Wang X. Y., Song J. M. Synchronization of the fractional order hyperchaos Lorenz systems with activation feedback control. Communications in Nonlinear Science and Numerical Simulation, Vol. 14, Issue 8, 2009, p. 33513357.

Yang X. J. Fractional derivatives of constant and variable orders applied to anomalous relaxation models in heattransfer problems. Thermal Science, Vol. 21, Issue 3, 2017, p. 11611171.

Yang X. J., Machado J. A. T. A new fractional operator of variable order: application in the description of anomalous diffusion. Physica A: Statistical Mechanics and its Applications, Vol. 481, 2017, p. 276283.

Li W., Zhang M. T., Zhao J. F. Stochastic bifurcations of generalized Duffingvan der Pol system with fractional derivative under colored noise. Chinese Physics B, Vol. 26, 2017, p. 90501.

Niu J. C., Shen Y. J., Yang S. P. Analysis of duffing oscillator with timedelayed fractionalorder PID controller. International Journal of NonLinear Mechanics, Vol. 92, 2017, p. 6675.

Giresse T. A., Crépin K. T. Chaos generalized synchronization of coupled MathieuVan der Pol and coupled DuffingVan der Pol systems using fractional orderderivative. Chaos, Solitons and Fractals, Vol. 98, 2017, p. 88100.

Nayfeh A. H. Nonlinear Oscillations. Wiley, New York, 1973.

Chen Y. S. Nonlinear Vibrations. Higher Education Press, Beijing, 2002.

Li X. H., Hou J. Y., Shen Y. J. Slowfast effect and generation mechanism of Brusselator based on coordinate transformation. Open Physics, Vol. 14, Issue 1, 2016, p. 261268.

Li X. H., Hou J. Y. Bursting phenomenon in a piecewise mechanical system with parameter perturbation in stiffness. International Journal of NonLinear Mechanics, Vol. 81, 2016, p. 65176.

Belhaq M., Fahsi A. 2:1 and 1:1 frequencylocking in fast excited van der PolMathieuDuffing oscillator. Nonlinear Dynamics, Vol. 53, Issues 12, 2008, p. 139152.

Wen S. F., Shen Y. J., Yang S. P., Wang J. Dynamical response of MathieuDuffing oscillator with fractionalorder delayed feedback. Chaos, Solitons and Fractals, Vol. 94, 2017, p. 5462.

Wen S. F., Shen Y. J., Li X. H., Yang S. P. Dynamical analysis of Mathieu equation with two kinds of van der Pol fractionalorder terms. International Journal of NonLinear Mechanics, Vol. 84, 2016, p. 130138.

Xu Y., Li Y. G., Liu D., Jia W. T., Huang H. Responses of Duffing oscillator with fractional damping and random phase. Nonlinear Dynamics, Vol. 74, Issue 3, 2013, p. 745753.
About this article
The authors are grateful to the support by the National Natural Science Foundation of China (No. 11602152), and the department project in Hebei Province (14212202D).