Primary resonance and feedback control of the fractional Duffing-van der Pol oscillator with quintic nonlinear-restoring force

1.   Introduction

Fractional calculus is an ancient mathematical topic from the 17th century. Although Fractional calculus is a mathematical subject with a history of more than 300 years, its application in physics and engineering did not attract wide attention until recent decades. In the past twenty years, significant progress has been made in fractional calculus and fractional order dynamical systems ^[1,2,3,4]. Previous studies have shown that fractional calculus is a mathematical tool, which is applicable to anomalous social and physical systems with non-local, frequency- and history-dependent properties, as well as intermediate states such as soft materials that are neither ideal solids nor ideal fluid ^[5,6,7]. Moreover, applications of fractional calculus have been reported in many areas such as signal processing, image processing, automatic control and robotics ^[8,9,10,11]. Research has even applied fractional derivatives to specific fields of research: psychological and life sciences ^[12,13]. In a word, virtually no area of classical analysis has been left untouched by fractional calculus. These examples and many other similar samples perfectly clarify the importance of consideration and analysis of dynamical systems with fractional-order models.

For the reason that nonlinear vibration systems exist in a large number of engineering technology and scientific practice, the research of nonlinear vibration systems has attracted the attention of many engineering and technical personnel and researchers. As a combination of two typical nonlinear systems, the nonlinear part of the Duffing-van der Pol system contains not only the nonlinear damping term for maintaining self-excited vibration of the van der Pol system, but also the cubic nonlinear restoring force term of the Duffing system, which has particularly rich dynamic characteristics. The Duffing-van der Pol oscillator is also a typical system to simulate some physical phenomena, and its simple nonlinear structure has led to an in-depth study of its dynamic behavior, which can be used as a research model in many fields, such as engineering, electronics, biology and neurology ^{[14,15,16,17]}. Furthermore, the fractional Duffing-van der Pol oscillator is considered to be one of the most commonly used equations in various system designs, including dynamics, biology, seismology, etc. ^[18,19]. Therefore, it is of practical significance to further study the dynamic behavior of such systems.

Many scientists have done a lot of research on the nonlinear dynamic behavior of the Duffing-van der Pol equation, including integer and fractional order cases ^{[20,21,22,23,24,25]}. Specially, in ^[22], a three-dimensional autonomous Van der Pol-Duffing type oscillator is proposed. Using the Routh-Hurwitz stability criterion and the linear stability analysis of equilibrium points, it was found that there is a Hopf bifurcation in the three-dimensional autonomous Van der Pol-Duffing oscillator. In ^[23], the damping characteristics of two Duffing-van der Pol oscillators with damping terms described by fractional derivative and time delay are studied by the residue harmonic balance method. A Duffing-van der Pol oscillator having fractional derivatives and time delays is investigated by the residue harmonic method in ^[24]. The dynamical properties of fractional order Duffing-van der Pol oscillators are studied in ^[25], and the amplitude-frequency response equation of the primary resonance is obtained by the harmonic balance method.

According to current research results, many studies focus on only Duffing-van der Pol equations with lower power terms under different excitations, while relatively few studies focus on the system with high power nonlinear terms. The Duffing-van der Pol systems with high power nonlinear terms are more worthy of study, because in the process of simplifying the actual model to the Duffing equation or van der Pol equation, in order to reduce the error between the system and the actual process, the nonlinear terms need to be approximately expressed as higher powers by power series to study the richer dynamic characteristics of the system. In addition, for the reason that the nonlinear dynamic behavior of the Duffing-van der Pol oscillator is very rich in bifurcation, chaos, primary resonance, superharmonic or subharmonic resonance, stochastic resonance and other vibration characteristics, the dynamic characteristics of the fractional Duffing-van der Pol oscillator are worthy of further exploration and research. It should be noted here that the primary resonance of a fractional Duffing-van der Pol oscillator can be studied by perturbation methods such as harmonic balance method or averaging method. Although a multiple scale method has been commonly used in the study of integer-order models, they are rarely involved in the study of fractional-order related problems. Therefore, it is also very necessary to use the multiple scale method to study the resonance characteristics of fractional Duffing-van der Pol oscillators, for example, to understand how the system parameters affect the periodic solutions resulted from resonance of the fractional Duffing-van der Pol oscillator.

The paper is organized as follows: In Section 2, the approximately analytical solution is obtained by the multiple scale method. The amplitude frequency equation is obtained and the stability conditions for the steady-state solution are given, while numerical simulations of the time history and phase diagrams are also provided. In Section 3, the bifurcations of primary resonance are investigated and the influence of parameters on fractional-order models is investigated. In Section 4, the comparisons of the approximate analytical solution and numerical solution are fulfilled. Finally, concluding remarks and implications in Section 5 close the paper.

2.   Mathematical analysis

2.1.   Primary resonance

The fractional Duffing-van der Pol oscillator with quintic nonlinear-restoring force is given by a second-order non-autonomous differential equation as follows

where $x$ is the position coordinate and is a function of time $t$, $D_t^{q}x$ is the $q$-order derivative of $x$ with respect to time, $\mu$ is a positive fractional damping coefficient, $\omega_0$ is the natural frequency, $\alpha$, $\lambda$ and $\beta$ are coefficients of nonlinear restoring forces, the external excitation have amplitude and frequency parameters which are $f$ and $\omega$, respectively. There are many definitions available for the fractional-order derivative, and in this study, $D_t^{q}x$ with $0 < q\leq1$ is the Caputo's fractional derivative of $x(t)$ described by

in which $\Gamma(z)$ is Gamma function satisfying $\Gamma(z+1) = z\Gamma(z)$. Under certain conditions, primary resonance with limited amplitude occurs in the fractional Duffing-van der Pol oscillator (2.1) when $\omega = \omega_0\overset{{\text{def}}}{ = }\sqrt{\alpha}$ or $\omega\approx\omega_0$. That is, the primary resonance means the excitation frequency is close to the natural one.

In order to find an approximate resonant solution of Eq (2.1), the multiple scale method will be used, for which a small parameter $\varepsilon$ is required. Let $T_n = \varepsilon^n t, (n = 0, 1)$ be two independent time variables, an approximate solution of Eq (2.1) with small amplitudes in three time scales can be represented by

With such a small parameter $\varepsilon$, in order to balance the effect of the nonlinearity, damping and excitation, the following scaling is used

Then, Eq (2.1) becomes

The derivatives with respect to $t$ can be expressed in terms of the new scaled times $T_n$ using chain rule as follows:

in which $D_0^q = \frac{\partial^q}{\partial T_0^q}$, $D_n = \frac{\partial}{\partial T_n}$, $D_n^2 = \frac{\partial^2}{\partial T_n^2}$ and $n = 0, 1$.

For the primary resonance with $\omega = \omega_0: = \sqrt{\alpha}$ or $\omega\approx\omega_0$, a detuning parameter $\sigma$ describing the nearness of $\omega$ to $\omega_0$ is introduced by

then $\omega t = \omega_0T_0+\sigma T_1$. Substituting (2.2), (2.5b) and (2.5c) into (2.4) leads to the following equation

Equating the coefficients of the same power of $\varepsilon$, a set of linear differential equations are obtained

from which $x_0$ and $x_1$ can solved one-by-one, respectively. In this way, the resonant solution $x$ is dominated by $x_0$ and collected by $\varepsilon x_1$.

In order to obtain the primary resonance solutions, it is necessary to eliminate the secular term when solving the Eqs (2.7) and (2.8) one by one, so as to find the amplitude-frequency response equation, from which the influence of various parameters on the resonance solution can be studied.

The general solution of Eq (2.7) is of the form,

where $A(T_1)$ and $\overline{A}(T_1)$ are unknown functions, $\overline{A}(T_1)$ denotes the complex conjugate of $A(T_1)$.

To solve Eq (2.8), for sufficient large $t\gg1$, the $q$th-order $(0\leq q\leq1)$ derivative of ${\rm e}^{{\rm i}\Omega t}$ is approximated with

according to the "Short Memory Principle" ^[13,26], where

Substituting Eqs (2.9) and (2.10) into Eq (2.8) and using

the right-hand of Eq (2.8) becomes

where $NST$ stands for the terms that do not produce secular terms, $cc$ denotes the complex conjugate of the preceding terms.

In order that $x_1$ is periodic, the secular terms with ${\rm e}^{{\rm i}\omega_0T_0}$ must be zero, namely

To analyze the solution of Eq (2.11), it is convenient to express $A(T_1)$ in the polar form as

In addition, the Euler formula gives

Let $\label{eq14}\varphi\overset{{\text{def}}}{ = }\sigma T_1-\theta$. By separating the real and imaginary parts of Eq (2.11), the differential equations governing amplitude $a(T_1)$ and $\varphi(T_1)$ of $A(T_1)$ are as follows respectively

We are interested in the steady solutions satisfying $D_1a = 0$, $D_1\varphi = 0$ in (2.13a) and (2.13b), namely

Because $\sin^2\varphi+\cos^2\varphi = 1$, by eliminating $\varphi$ from Eqs (2.14a) and (2.14b), the amplitude-frequency response equation can be obtained as follows

The real solution $a$ of Eq (2.15) determines the dominant part $x_0$ of the primary resonance response amplitude. Similar results have been discussed in ^[21] for the case when $q = 1$, but the calculation error in the derivation process of this paper leads to inaccurate amplitude-frequency response equation.

2.1.1.   Feedback control

The fractional Duffing-van der Pol oscillator with negative acceleration feedback considered is as follows

in which $-u$ is the feedback gain reflecting the magnitude of a feedback force. Similar to the previous discussion, the method of multiple time scales is used to obtain a uniformly valid, asymptotic expansion of the solution for Eq (2.16). Let $T_n = \varepsilon^n t, (n = 0, 1)$ be two independent time variables, an approximate solution of Eq (2.16) with small amplitudes in three time scales can be represented by

With such a small parameter $\varepsilon$, in order to balance the effect of the nonlinearity, damping and excitation, the following scaling is used

Then, Eq (2.16) becomes

Substituting (2.17), (2.5b) and (2.5c) into (2.19) leads to the following equation

Similar to the previous discussion, the following differential equations about amplitude and phase can be obtained,

The steady solutions satisfying $D_1a = 0$, $D_1\varphi = 0$ in (2.20a) and (2.20b), namely

By eliminating $\varphi$ from Eqs (2.21a) and (2.21b), the amplitude-frequency response equation can be obtained as follows

It should be noted that a fractional dynamical system does not have an exact periodic solution, see ^[27,28], and it may have one if the fractional derivative is defined over $(-\infty, t]$ ^[28]. When the asymptotic behaviors of the fractional system are addressed, a fractional derivative over $[0, t]$ with large $t$ equals approximately the corresponding one over $(-\infty, t]$, according to the "Short Memory Principle" of fractional derivatives ^[26]. Thus, a periodic solution of a fractional system should be understood as an asymptotic solution with transients before a finite time neglected. All the periodic solutions mentioned in this paper are understood in this sense.

2.2.   Stability of the steady resonant solutions

To determine the stability of the steady resonant solutions, suppose that $(a, \varphi) = (a^*, \varphi^*)$ is a steady solution of Eqs (2.20a) and (2.20b), let

are perturbations which are assumed to be small compared to $a^*$ and $\varphi^*$. Substituting Eq (2.23) into Eqs (2.20a) and (2.20b), using Eqs (2.21a) and (2.21b), then the linearized differential equations governing ${\rm\Delta}a$ and ${\rm\Delta}\varphi$ are

Let $M$$=$$\sigma$$+$$\frac{2-(a^*)^2}{4}\mu\omega_0^{q-1}\cos\frac{q\pi}{2}$$-$$\frac{3\lambda}{8\omega_0} (a^*)^2$$-$$\frac{5\beta}{16\omega_0} (a^*)^4$$+$$\frac{u\omega_0}{2}$, $N$$=$$\sigma$$+$$\frac{2-3(a^*)^2}{4}\mu\omega_0^{q-1}\cos\frac{q\pi}{2}$$-$$\frac{9\lambda}{8\omega_0}$$(a^*)^2-$$\frac{25\beta}{8\omega_0}$$(a^*)^4$$+$$\frac{u\omega_0}{2}$. According to Eqs (2.21a) and (2.21b), the characteristic equation can be written as

By expanding the determinant, we can get

in which $\Lambda = \frac{[2-3(a^*)^2][2-(a^*)^2]}{16}(\mu\omega_0^{q-1}\sin\frac{q\pi}{2})^2+MN$.

Thus, the steady solution $(a, \varphi) = (a^*, \varphi^*)$ is asymptotically stable if and only if the real parts of all roots of the characteristic Eq (2.25) are negative.

2.3.   Numerical simulation

In this section, time history diagrams are used to briefly illustrate the impact of feedback control on the model. Figure 1 shows the time history of the system without controller. The various parameters of the system in Figure 1 are $\omega_0 = 1$, $q = 0.9$, $\mu = 0.01$, $\lambda = 0.04$, $\beta = 0.01$ and $f = 0.06$. It can be observed that the amplitude fluctuates significantly in the first 150s and the steady-state amplitude of the system is about 0.28. The phase diagram shows a limit cycle, indicating that there is no chaos in the system. Different initial conditions were tested and it was found that the steady-state amplitude of the system is not sensitive to the initial conditions.

Figure 2 shows the time history of the system with feedback controller. The various parameters of the system in Figure 2 are $\omega_0 = 1$, $q = 0.9$, $\mu = 0.01$, $\lambda = 0.04$, $\beta = 0.01$, $f = 0.06$ and $u = 0.1$. From Figure 2, it can be seen that the steady-state amplitude of the system has decreased to approximately 0.115, this indicates the effectiveness of the controller. The steady-state value appears after 150s without control, while it appears before 50 seconds with control, which also reflects the effect of adding control.

3.   Analysis of the resonant solutions

Rewritten Eqs (2.15) and (2.22) in terms of the original system parameters as follows

and

The effect of different parameters on the amplitude of the resonant solutions is investigated numerically. Here, the difference between the model with quintic nonlinear term and the model without quintic power is given by the amplitude-frequency response curve, as shown in Figure 3.

It can be seen from the figure that the bending degree, resonance region and jump phenomena of the model have obvious changes after considering the quintic nonlinear terms, indicating that it is very important to consider the role of higher order terms for studying the nonlinear characteristics of the model.

3.1.   Parameter influence on the amplitudes of the resonant solutions

This section discusses the influence of parameters on the hardening/softening characteristics and resonance peak value based on Eq (3.1). As shown in Figures 4–7, typical hardening/softening characteristics can be observed also in the curves of the resonant amplitudes with respect to the changes of $\lambda$, $\beta$, $q$, $\mu$ and $f$, respectively.

The influence of the nonlinear terms on the amplitude-frequency response curve can be identified with the changes of $\lambda$, $\beta$. Typical softening characteristic can be observed in Figure 4 with $\lambda = -0.04$ and $\beta = -0.01$, and the frequency response curve bends to low frequency. When $\lambda = 0.04$ and $\beta = 0.01$, hardening characteristic can be observed and the frequency response curve bends to high frequency.

The effect of the fractional-order $q$ on the amplitude-frequency curves is shown in Figure 5. If the fractional order $q$ is closer to 0, the fractional differential term degenerates to linear stiffness, while if $q$ is closer to 1, it degenerates to linear damping. It is further known that the smaller the order $q$ is, the larger the maximum amplitude is. The integer-order model is the lowest. Compared with the integer order case when $q = 1$, the bending degree, resonance peak and resonance region of the amplitude-frequency curve of the fractional-order model change accordingly with the decrease of the fractional-order $q$.

Figure 6 shows the effect of the damping coefficient $\mu$ on the amplitude-frequency curves. With the increase of $\mu$, the nonlinear jump of the system weakens and the resonance amplitude of the system decreases. In other words, as $\mu$ increases, the unstable portions decrease.

The effect of the excitation amplitude $f$ on the amplitude-frequency curves can be seen in Figure 7. With the increase of the excitation amplitude $f$, the frequency response curve will generate a multivalued region, the nonlinear jump of the system is more obvious and the resonance range and resonance amplitude of the system increase.

Figures 5–7 present parameter effects on resonant amplitudes with respect to excitation frequency. With fixed parameter values, all the figures exhibit typical characteristics of hardening spring. Let us take a different view on the resonant solutions. Now, the excitation frequency is fixed, and then with the change of a parameter, the amplitude of the primary resonant solution changes in different ways. In Figures 8–10, the images reflect characteristics similar to hardening spring, where the peak of amplitude is bent to the right as $f$ varies, as $\lambda$ varies or as $\beta$ varies, respectively.

3.2.   Parameter influence on the resonant solutions under feedback control

Through Eq (3.2), this section discusses the influence of parameters on the resonant solutions under feedback control. Figure 11 considers the influence of the gain value of the feedback control on the vibration state of the system. The red line in Figure 11 shows the primary resonance amplitude frequency response curve of the system without feedback control. It can be seen from Figure 11 that compared with the amplitude frequency curve of the primary resonance of the system without feedback control, after adding the feedback control, the resonance curve of the system moves to the left, and the resonance domain of the system changes. When the amplitude of the value of the gain of the control is increased, the deviation of the amplitude frequency response curve becomes larger.

Figure 12 shows the influence of the frequency parameters on the vibration state of the system. The closer the external excitation frequency $\omega$ is to the primary resonance frequency $\omega_0$, the smaller the offset of amplitude frequency curve is, and the smaller the change of resonance domain is. On the contrary, the more the phase difference, the larger the offset of the amplitude frequency curve. The peak amplitude of the primary resonance in Figure 12 does not change, but changes the position where it appears.

4.   Comparison between approximate analytical solution and numerical solution

According to Eq (3.1), the primary resonance amplitude-frequency response curve of the system can be drawn. For comparison, this paper adopts the power series method introduced in references [3,29], and its calculation formula is

where $t_n = nh$ is the sample points, $h$ is the sample step and $C_j^q$ is the fractional binomial coefficient with the iterative relationship as

According to Eqs (4.1) and (4.2), the numerical scheme for Eq (2.4) can be expressed as

The numerical amplitude-frequency curve marked with circle in Figure 13, where the stepsize of time is $h = 0.005$, and the total computation time is 100s with the first 25s neglected. It shows that the resonant amplitude calculated from Eq (3.1) is in good agreement with the numerical results, especially when the $\omega\approx\omega_0$.

5.   Conclusions

We study the primary resonance and feedback control of the fractional Duffing-van der Pol oscillator with quintic nonlinear-restoring force. The approximate analytical solution and the amplitude-frequency equation are obtained by the multiple scale method. The stability condition for steady-state solution is obtained based on the Lyapunov theory. The influence of parameters on system characteristics was studied using the amplitude frequency response equation and phase plane technology, and bifurcation analysis was conducted to verify the stability of the system. It can also be concluded that fractional orders and damping coefficients are very important in the fractional Duffing-van der Pol oscillator. For example, larger fractional orders and larger damping coefficients can reduce the effective amplitude of resonance and change the resonance frequency.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author expresses gratitude to the anonymous reviewers for their helpful comments and suggestions, as well as to the financial support from High Level Talent Research Launch Fund of Jiangsu Second Normal University (No. 928201/058).

Conflict of interest

The author declares that the research was conducted in the absence of any conflict of interest.