Immediate solution for fractional nonlinear oscillators using the equivalent linearized method

Abstract

The existence of the derivative with a fractional-order in a class of differential equations could lead to complicating the analysis. In this paper, a novel approach has been introduced to facilitate the analysis of the oscillation having fractional-order derivatives and to obtain the analytical solution easily. The present technique is formulated to provide an easy way of understanding. The suggested technique has been utilized to study two examples for illustration. The similitude between the analytical and numerical solution verifies and gives satisfactory precision to the equivalent solution. The new technique is represented to be the best tool for solving the nonlinear oscillation problems in physics and engineering which have a fractional-order.

Introduction

The dynamics of real-world complicated problems are mostly related to fractional calculus. Due to the non-local nature of fractional operators, several natural phenomena are described systematically.^1–3 In this regard, fractional-order differential equations are the generalization of the classical integer-order differential equations.^4–6 The fractional approach is now considered to be the most powerful tool for modeling and it is applied to many physical and engineering problems.^7,8 It plays a very important role in modeling the frequency-dependent damping behavior of many applications in biomedical engineering, finance, probability theory and hydrology, electrochemistry, computational biology, physics, fluid mechanics, and many other fields.^9,10 Much literature focused on simple fractal oscillators.^11,12

In many nonlinear dynamical systems, such as van der-Pol, Duffing oscillators, Toda oscillators, Klein-Gordon oscillators, and other systems, a large number of results can be found in the literature.^13–16 The numerical methods ^17–20 and analytical perturbation solutions^21–24 of fractional oscillators are found to be hard work of research.^25–28 Therefore, the researchers have given full attention to developing strong and novel techniques to handle this complicated class of fractional mathematical tools. ²⁹ In this regard, a new technique is urgent to simplify and reduce the hard work required for obtaining an asymptotic solution that is closer to the exact numerical solution.

In the current work, a new simple technique will be introduced to reduce the hard work. The approach depends on replacing the original fractional equation with an equivalent linear oscillation with the integer derivatives. Based on the awareness of the non-perturbative method,^30–32 one could obtain a good asymptotic solution for a fractional-order dynamical system, and then analyze some important properties such as the amplitude-frequency equation. Accordingly, it may present more information about fractional-order differential equations. Galerkin’s method for fractional expression has been introduced. The non-perturbative approach has been applied to two examples for illustration. The fractional damping Duffing oscillator ³³ or oscillator with full fractional-order derivatives ³⁴ is analytically studied by the equivalent linearized method. Also, two kinds of damping fractional-order derivatives of the van-der-Pol-Duffing ³⁵ oscillator can be easily handled by using the present technique.

Fractional problem statement

The goal is to try to obtain the solution of the fractional oscillator by a simple approach based on the equivalent linearized method. In this section, a new approach is introduced to convert the fractional nonlinear oscillation to an equivalent linear oscillator with an integer derivative. The established method leads to finding a quasi-exact analytic solution for the fractional nonlinear oscillator. Thus, the following class of fractional nonlinear oscillators

$$f\left(D^{\alpha +1}y,D^{\alpha }y,y\right)=0;\hspace{1em}0<\alpha <1$$

(1)

is aimed to turn into the following ordinary linear differential equation of integer-order concerning the same independent variable

$$g\left(\ddot{y},\dot{y},y\right)=0$$

(2)

where $f$ is a polynomial in $y\left(t\right)$ and the various derivatives of the fractional-order of $y\left(t\right)$ besides the nonlinearity part. The function $g$ is the equivalent linear polynomial $y\left(t\right)$ with derivatives of integer-order represented by the dots. The fractional derivative operator $D^{\alpha }$ refers to the Riemann-Liouville definition.^1,4 The following fractional derivatives of trigonometric functions are used

$$D^{\alpha }\cos \mathrm{\Omega }t={\mathrm{\Omega }}^{\alpha }\cos \left(\mathrm{\Omega }t+\frac{1}{2}\pi \alpha \right)$$

(3)

$$D^{\alpha }\sin \mathrm{\Omega }t={\mathrm{\Omega }}^{\alpha }\sin \left(\mathrm{\Omega }t+\frac{1}{2}\pi \alpha \right)$$

(4)

Since equation (2) is assumed to be a linear second-order damping oscillator, then it is convent to introduce the trail solution in the form

$$\widehat{y}\left(t\right)=A\cos \left(\mathrm{\Omega }t\right)$$

(5)

where $\mathrm{\Omega }$ is assumed to represent the total frequency that is controlled by the fractional nonlinear oscillator, $A$ is the amplitude of the oscillation (2). It is noted that in the limit case for $\alpha \to 1$ , then $D^{\alpha }y$ will contribute to the damping forces while $\alpha \to 0$ should produce a contribution to the stiffens forces. In general, the fractional derivative $D^{\alpha }y$ has contributions to both damping and stiffening forces.

Since the trial solution (5) is not exact, then there will be a residual or error remaining. Based on the principle of minimum mean-square error, ³⁶ the located displacement point, and the velocity point are estimated as given below

$${\overline{y}}^2=\frac{1}{2T}{\displaystyle {\int }_0^T{\widehat{y}}^2\left(t\right)dt}=\frac{1}{4}{A}^2;\hspace{1em}T=\frac{2\pi }{\mathrm{\Omega }}.$$

(6)

$${\overline{\dot{y}}}^2=\frac{1}{2T}{\displaystyle {\int }_0^T{\widehat{\dot{y}}}^2\left(t\right)dt}=\frac{1}{4}{\mathrm{\Omega }}^2{A}^2.$$

(7)

When the oscillator is described by the fractional derivative $D^{\alpha }y$ , then the located displacement and velocity points due to the fractional contributions are denoted by ${\overline{y}}^{\alpha }$ and ${\overline{\dot{y}}}^{\alpha }$ which could be estimated based on the weight residual method ³⁶ as

$${\overline{y}}^{\alpha }=\frac{1}{AT}{\displaystyle {\int }_0^T\widehat{y}\left(t\right)D^{\alpha }\widehat{y}\left(t\right)dt}=\frac{1}{2}A{\mathrm{\Omega }}^{\alpha }\cos \left(\frac{1}{2}\pi \alpha \right),$$

(8)

$${\overline{\dot{y}}}^{\alpha }=\frac{1}{A\mathrm{\Omega }T}{\displaystyle {\int }_0^T\dot{y}\left(t\right)D^{\alpha }y\left(t\right)dt}=\frac{1}{2}A{\mathrm{\Omega }}^{\alpha }\sin \left(\frac{1}{2}\pi \alpha \right).$$

(9)

To highlight the main results of this work and to show the effectiveness of solving nonlinear fractional oscillators, two illustrative examples are studied to illustrate the aims: the fractional-Duffing oscillator ³³ and the fractional Van-der-Pol-Duffing oscillator. ³⁴

Example 1: Duffing oscillator with the fractional-order

This example is concerned with how to solve the following class of second-order fractional-Duffing equations ³⁷ using the non-perturbative approach

$$\ddot{y}+\mu D^{\alpha }y+{\omega }_0^{2}y+Q{y}^3=0;\hspace{1em}0<\alpha <1$$

(10)

where $\mathit{\mu },{\mathit{\omega }}_0^2$ and Q are constants. Assuming that the above oscillator, subject to the initial conditions $y\left(0\right)=A,\hspace{0.17em}\dot{y}\left(0\right)=0$ . The fractional derivative obeys the definition of the Riemann-Liouville time-fractional derivative (3) and (4). As mentioned above the fractional $D^{\alpha }y$ will be replaced by

$$D^{\alpha }y=\left(\frac{D^{\alpha }y}{\dot{y}}\right)\dot{y}+\left(\frac{D^{\alpha }y}{y}\right)y.$$

(11)

Therefore, equation (10) may be rewritten in the form

$$\ddot{y}+\mu \left(\frac{D^{\alpha }y}{\dot{y}}\right)\dot{y}+\mu \left(\frac{D^{\alpha }y}{y}\right)y+{\omega }_0^{2}y+Q{y}^3=0.$$

(12)

The equivalent integer-order can be formulated using the principle of the minimum mean-square error to yield

$$\ddot{y}+\mu a_{eq}\dot{y}+\left(\mu b_{eq}+{\omega }_0^2\right)y+Q{y}^3=0$$

(13)

where $a_{eq}$ and $b_{eq}$ are estimated in terms of the trail function (5) to be

$$a_{eq}=\frac{2}{A\hspace{0.17em}\mathrm{\Omega }T}{\displaystyle {\int }_0^T\dot{y}\left(t\right)(\frac{D^{\alpha }y\left(t\right)}{\overline{\dot{y}}})dt}={\mathrm{\Omega }}^{\alpha -1}\sin \left(\frac{1}{2}\pi \alpha \right),$$

(14)

$$b_{eq}=\frac{2}{AT}{\displaystyle {\int }_0^{T}y\left(t\right)\left(\frac{D^{\alpha }y\left(t\right)}{\overline{y}}\right)dt}={\mathrm{\Omega }}^{\alpha }\cos \left(\frac{1}{2}\pi \alpha \right).$$

(15)

At this end, equation (13) can be linearized and solved with the non-perturbative approach.^31–33

In the simplest approach, the equivalent linearized form of equation (10) can be sought in the form

$$\ddot{y}+\mu a_{eq}\dot{y}+{\omega }_{eq}^{2}y=0$$

(16)

To estimate the coefficients $a_{eq}$ and ${\omega }_{eq}^2$ , we may assume that there is a total frequency $\mathrm{\Omega }$ governing the fractional influence of the above oscillator. Accordingly, $a_{eq}=a_{eq}\left(\mathrm{\Omega }\right)$ and ${\omega }_{eq}^2={\omega }_{eq}^2\left(\mathrm{\Omega }\right)$ . Employing equations (7) and (9) the coefficient $a_{eq}$ is given by

$$a_{eq}\left(\mathrm{\Omega }\right)={\left(\frac{D^{\alpha }y}{\dot{y}}\right)|}_{\begin{array}{l}\dot{y}=\frac{1}{2}\mathrm{\Omega }A,\\ {\dot{y}}^{\alpha }=\frac{1}{2}A\hspace{0.17em}{\mathrm{\Omega }}^{\alpha }\hspace{0.17em}\sin \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)\end{array}}={\mathrm{\Omega }}^{\alpha -1}\hspace{0.17em}\sin \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right).$$

(17)

In terms of the restoring force

$$f\left(y\right)=\left({\omega }_0^2+\mu \left(\frac{D^{\alpha }y}{y}\right)\right)y+Q{y}^3$$

(18)

and based on He’s frequency formula,^38–41 ${\omega }_{eq}^2\left(\mathrm{\Omega }\right)$ could be derived using equations (6) and (8) in the form

$${\omega }_{eq}^2\left(\mathrm{\Omega }\right)={\frac{df\left(y\right)}{dy}|}_{\begin{array}{l}y=\frac{1}{2}A\\ y^{\alpha }=\frac{1}{2}A\hspace{0.17em}{\mathrm{\Omega }}^{\alpha }\hspace{0.17em}\cos \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)\end{array}}={\omega }_0^2+\frac{3}{4}Q{A}^2+\mu {\mathrm{\Omega }}^{\alpha }\hspace{0.17em}\cos \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right).$$

(19)

The linear damping equation (16) has the following solution

$$y\left(t\right)=A{e}^{-\frac{1}{2}\mu a_{eq}t}\cos \mathrm{\Omega }t$$

(20)

where $A$ is the amplitude of the oscillator The solution (20) is called the quasi-exact solution for the fractional nonlinear Duffing oscillator (10), the total frequency $\mathrm{\Omega }$ is given by

$${\mathrm{\Omega }}^2={\omega }_{eq}^2\left(\mathrm{\Omega }\right)-\frac{1}{4}{\mu }^2{a}_{eq}^2\left(\mathrm{\Omega }\right)$$

(21)

Employing equations (17) and (19) into equation (21) gives the explicit frequency in the form

$${\mathrm{\Omega }}^2=\left({\omega }_0^2+\frac{3}{4}Q{A}^2\right)+\mu {\mathrm{\Omega }}^{\alpha }\hspace{0.17em}\cos \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)-\frac{1}{4}{\mu }^2{\mathrm{\Omega }}^{2\alpha -2}\hspace{0.17em}{\sin }^2\left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right).$$

(22)

The above equation is the frequency-amplitude relationship which contains the frequency $\mathrm{\Omega }$ with the fraction power. The perturbation technique is urgent to solve equation (22). Therefore, introducing a small parameter $\delta$ into equation (22) becomes

$${\mathrm{\Omega }}^2=\left({\omega }_0^2+\frac{3}{4}Q{A}^2\right)+\delta \left(\mu {\mathrm{\Omega }}^{\alpha }\hspace{0.17em}\cos \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)-\frac{1}{4}{\mu }^2{\mathrm{\Omega }}^{2\alpha -2}\hspace{0.17em}{\sin }^2\left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)\right).$$

(23)

The frequency $\mathrm{\Omega }\left(\delta \right)$ is expanded as

$${\mathrm{\Omega }}^2\left(\delta \right)={\mathrm{\Omega }}_0^2+\delta {\mathrm{\Omega }}_1+\mathrm{\dots }$$

(24)

where ${\mathrm{\Omega }}_0\hspace{0.17em}\&\hspace{0.17em}\hspace{0.17em}{\mathrm{\Omega }}_1$ are unknowns to be determined by inserting equation (24) into equation (23), and equating the identical powers of $\delta$ to zero yields

$${\mathrm{\Omega }}_0^2={\omega }^2;\hspace{1em}{\omega }^2=\left({\omega }_0^2+\frac{3}{4}Q{A}^2\right)\mathrm{and}$$

(25)

$${\mathrm{\Omega }}_1=\mu {\mathrm{\Omega }}^{\alpha }\hspace{0.17em}\cos \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)-\frac{1}{4}{\mu }^2{\mathrm{\Omega }}_0^{2\alpha -2}\hspace{0.17em}{\sin }^2\left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right).$$

(26)

Substituting equations (25) and (26) into equation (24) yields the approximate frequency formula that could be expressed as

$${\mathrm{\Omega }}^2={\omega }^2+\mu {\omega }^{\alpha }\hspace{0.17em}\cos \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)-\frac{1}{4}{\mu }^2{\omega }^{2\alpha -2}\hspace{0.17em}{\sin }^2\left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right).$$

(27)

Now, the function (20) and the frequency formula (27) represent the approximate solution of the fractional damping Duffing equation (10).

Remark

In the current example when the operators in Duffing equation (10) become of full fractional forms as

$$D^{\mathrm{\alpha }+1}y+\mathit{\mu }{D}^{\mathrm{\alpha }}y+{\mathit{\omega }}_{\mathit{0}}^{\mathit{2}}y+Q{y}^3=0$$

(28)

then it becomes difficult to solve equation (28) directly by the traditional method. There are no suitable analytical solutions that use perturbation techniques that can be successfully asymptotic to their numerical solution. It is worthwhile to note that the knowledge of the equivalent linearized method will lead to obtaining an approximate solution close to the numerical solution. This will be illustrated in what follows:

It is noted that the solution of equation (28) as a first-order differential equation is not suitable because the original equation has order greater than the first-order. Without the fractional-order derivatives and based on the methods mentioned above, it is convenient to convert equation (28) into the following second-order equation

$$a_{eq}\ddot{y}+\left(\mu a_{eq}+b_{eq}\right)\dot{y}+{\omega }_{eq}^{2}y=0$$

(29)

where the following substitution is used

$$D^{\alpha +1}y=\left(\frac{D^{\alpha }y}{\ddot{y}}\right)\ddot{y}+\left(\frac{D^{\alpha }y}{\dot{y}}\right)\dot{y}$$

(30)

The parameters in equation (29) are evaluated with those of equations (14) and (15). equation (29) can be simply rewritten in the form

$$\ddot{y}+\left(\mu +\frac{b_{eq}}{a_{eq}}\right)\dot{y}+\frac{{\omega }_{eq}^2}{a_{eq}}y=0,$$

(31)

The solution of the equivalent equation (31) and its corresponding frequency in terms of the values of $a_{eq}$ , $b_{eq}$ and ${\omega }_{eq}^2$ has the form

$$y\left(t\right)=A{e}^{-\frac{1}{2}\left(\mu +\mathrm{\Omega }\cot \left(\frac{1}{2}\pi \alpha \right)\right)t}\cos \mathrm{\Omega }t$$

(32)

and

$${\mathrm{\Omega }}^2=\frac{{\omega }^2}{{\mathrm{\Omega }}^{\alpha -1}\sin \left(\frac{1}{2}\pi \alpha \right)}+\mu \mathrm{\Omega }\hspace{0.17em}\cot \left(\frac{1}{2}\pi \hspace{0.17em}\alpha \right)-\frac{1}{4}{\left(\mu +\mathrm{\Omega }\cot \left(\frac{1}{2}\pi \alpha \right)\right)}^2.$$

(33)

The obtained solution (32) with its frequency formula (33) will be numerically illustrated.

Comparison with the numerical solution of the traditional integer-order

Comparing the solution of the fractional-order with the traditional integer-order can be accomplished, especially, for the two cases of $\alpha =1$ or $\alpha =0$ into the original equation (10) and the resulting solution (20). The comparison with the solution obtained by HPM could be considered for the case of $0<\alpha <1$ . Considering $\alpha =1$ yields equation (10) in the form

$$\ddot{y}+\mu \dot{y}+{\omega }_0^{2}y+Q{y}^3=0$$

(34)

Also, the limiting case of solution (20) is reduced to

$$y\left(t\right)=A{e}^{-\frac{1}{2}\mu t}\cos \left(\sqrt{{\omega }^2-\frac{1}{4}{\mu }^2}t\right)$$

(35)

In addition, considering $\alpha =0$ will reduce equation (10) to the un-damping Duffing equation

$$\ddot{y}+\left(\mu +{\omega }_0^2\right)y+Q{y}^3=0$$

Consequently, the periodic solution arises. In what follows, a numerical illustration is made to compare the obtained quasi-exact solution (20) with the numerical solution. The numerical simulation is performed using the Mathematica software. The calculations are displayed in Figures 1 and 2. In Figure 1, the comparison is made for a system having $\mu =0.3,\hspace{0.17em}Q=0.1,\hspace{0.17em}{\omega }_0=1.2,\hspace{0.17em}A=1$ , and $\alpha \to 1$ . The graph shows that there is an excellent agreement. In Figure 2, the parameter $\alpha$ is selected to be zero. This graph shows that the oscillation has the same amplitude. This is expected because $\alpha =0$ makes equation (10) is coincident with the un-damped Duffing equation. This case represents non-damped behavior in which the system oscillates at its natural resonant frequency. From graphs (1) and (2), it is confirmed that the exact and derived results are in close contact with each other. Thus, the proposed method has provided an accurate solution for equation (10). In the case of the parameter $0<\alpha <1$ , the examination will be done by comparison with the solution of the HPM.

Figure 1.

Comparison of the quasi-exact solution (20) (Red-dashing curve) with the numerical solution of equation (10) (Green-solid curve) in the case of $\alpha =1$ , the system of $\mu =0.3,\hspace{0.17em}\hspace{0.17em}Q=0.1,\hspace{0.17em}\hspace{0.17em}{\omega }_0=1.2,\hspace{0.17em}\mathrm{and}\hspace{0.17em}A=1$ .

Figure 2.

Comparison of the quasi-exact solution with the numerical solution in the case of $\alpha =0$ , the same system given in Figure 1.

To examine the influence of the variation in the parameter $0<\alpha <1$ , the comparison may occur with the solution resulting by applying a perturbation solution such as the homotopy perturbation method. The solution of this oscillator was obtained before, by using HPM, in Ref. 33 which is found to be

$$y\left(t\right)=A{e}^{-\hspace{0.17em}\frac{1}{2}t\mu {\mathrm{\Omega }}^{\alpha -1}\sin \left(\frac{1}{2}\pi \alpha \right)}\left(\cos \left(\mathrm{\Omega }t\right)+\frac{1}{32{\mathrm{\Omega }}^2}{A}^{2}Q\cos \left(3\mathrm{\Omega }t\right)\right)$$

(36)

where the corresponding total frequency $\mathrm{\Omega }$ was given by

$${\mathrm{\Omega }}^2={\omega }^2+\mu {\mathrm{\Omega }}^{\alpha }\cos \left(\frac{1}{2}\pi \alpha \right)$$

(37)

The comparison between the solutions (20) and their frequency (22) with the solution (36) and its frequency (37) shows that the damping factor in both solutions is identical, while the frequency (22) is more accurate than the frequency (37).

To demonstrate the accuracy to achieve the solution (20), the comparison with the result obtained by the HPM is plotted at various fractional-orders of the derivatives in Figures 3–7. The response diagram in these graphs shows that the decrease in the parameter $\alpha$ leads to suppressing the damping effect until this damping influence vanishes as shown in Figure 2. The relative error between the numerical solution of equation (10) and the analytical solution (20) has been displayed in Table 1 for Figure 1 and 2. Also, the relative error between the analytical solution (20) and the HPM-solution (36) has been demonstrated in Table 1 corresponding to the Figures 3–6.

Figure 3.

Comparison of the solution (20) with the solution (36) for the same system given in Figure 1 in the case of $\alpha =0.9$ .

Figure 4.

Comparison of the solution (20) with the solution (36) for the same system given in Figure 1 in the case of $\alpha =0.7$ .

Figure 5.

Comparison of the solution (20) with the solution (36) for the same system given in Figure 1 in the case of $\alpha =0.5$ .

Figure 6.

Comparison of the solution (20) with the solution (36) for the same system given in Figure 1 in the case of $\alpha =0.3$ .

Figure 7.

Comparison of the solution (20) with the solution (36) for the same system given in Figure 1 in the case of $\alpha =0.1$ .

Table 1.

Illustrate the relative error of the analytical solution (20).

$t$	Relative error Figure 1	Relative error Figure 2	Relative error Figure 3	Relative error Figure 4	Relative error Figure 5	Relative error Figure 6
0	0	0	0.01366	0.001865	0.001782	0.001734
1	0.2522	7.422 × 10−8	−0.1323	−0.01271	0.01974	0.06032
2	−0.1149	−3.633 × 10−8	0.01443	0.002636	−0.006961	−0.01507
3	0.06617	2.521 × 10−8	−0.1285	−0.007275	0.01793	0.05304
4	−2.302	−1.853 × 10−7	0.01673	0.01381	−0.04221	−0.08234
5	−0.02216	3.817 × 10−8	−0.1212	−0.00003	0.01122	0.03565
6	0.1753	1.967 × 10−7	0.02122	−0.1441	−0.4549	−0.4692
7	−0.1716	2.905 × 10−8	− 0.1108	0.00739	−0.006209	0.0006607
8	0.04273	1.771 × 10−7	0.02908	−0.02528	0.1108	0.4522
9	0.9963	−1.018 × 10−6	−0.09788	0.02409	−0.04506	−0.06803
10	−0.04686	1.688 × 10−7	0.04229	−0.006585	0.04798	0.1781
11	0.1275	3.855 × 10−7	−0.08295	0.2669	−0.1612	−0.2294
12	−0.2706	6.071 × 10−8	0.06412	0.007692	0.01456	0.08705
13	0.02111	2.811 × 10−7	− 0.0666	−0.05944	1.022	−1.221
14	0.4095	3.129×10−6	0.1004	0.02793	−0.02796	0.003561
15	−0.07644	2.419 × 10−7	−0.04932	−0.01872	0.1436	0.7583
16	0.09326	5.016 × 10−7	0.163	0.1224	−0.109	−0.1189
17	−0.5112	4.2 × 10−8	−0.03148	0.00383	0.0604	0.2973
18	−0.0001	4.159 × 10−7	0.2821	−0.1386	−0.4431	−0.3869
19	0.2517	1.473 × 10−6	−0.01328	0.02783	0.005262	0.1408
20	−0.1151	3.037 × 10−7	0.5688	−0.03853	0.5065	−2.376
Averring of the total error is	−0.0569	3.002 × 10−7	0.02293	0.002623	0.0317	−0.138

In Figure 8, when the order changes from integer to fractional, the amplitude of the curve of the fractional differential model increases. From the variation in the curves in Figure 8, it could be seen that the decrease of $\alpha$ the behavior tends to be the periodic solution. Figure 9 illustrates the solution of equation (24) relative to the numerical solution of equation (30) in the case of $\alpha =1$ . The variation of the parameter $\alpha$ of the quasi-exact solution (32) is displayed in Figure 10. As can be seen from Fig.(10), there are different sizes in the amplitude in the time history, where the amplitude value decreases as $\alpha$ decreases. At the smallest $\alpha$ , there is no oscillation observed especially for $\alpha <0.5$ . This is known as the over-damped oscillator.

Figure 8.

The variation of the parameter $\alpha$ in the solution (20) for the same system as given in Figure 1.

Figure 9.

Comparison of the solution (32) with the numerical solution of (28) for the same system given in Figure 1 in the case of $\alpha =1$ .

Figure 10.

The variation of the parameter $\alpha$ in the solution (32) for the same system as given in Figure 1.

Example 2

Application of the suggested method to the fractional Van-der-Pol-Duffing oscillator

In this example, the quasi-exact analytical solution of the problem van-der-Pol Duffing oscillator will be discussed. The problem is selected with two kinds of fractional-order derivatives. The solution will be obtained based on the equivalent approach established in the present paper.

Consider the fractional Van-der-Pol-Duffing oscillator having two kinds of fractional-order ³⁵ given in the following form

$$\ddot{x}+\left(1-x^2\right)\left(\mu D^{\alpha }x+\eta D^{\alpha +1}x\right)+\left({\omega }_0^{2}x+Q{x}^3\right)=0;\hspace{0.17em}\hspace{0.17em}y\left(0\right)=A,\dot{y}\left(0\right)=0$$

(38)

where $\mu ,\hspace{0.17em}\eta ,\hspace{0.17em}{\omega }_0$ and $Q$ are constants. The main step of the above-mentioned method is to obtain the equivalent form of equation (38) with an integer-order. This method could be transformed in equation (38) into the following configuration

$$\ddot{x}+\frac{\left(1-x^2\right)\left(\mu a_{eq}+\eta b_{eq}\right)}{1+\left(1-x^2\right)\eta a_{eq}}\dot{x}+\frac{\left({\omega }_0^2+\mu b_{eq}\right)x+\left(Q-\mu b_{eq}\right)x^3}{1+\left(1-x^2\right)\eta a_{eq}}=0$$

(39)

This configurative represents a nonlinear oscillator having a nonlinear damping force with a restoring force $f\left(x\right)$ given by

$$f\left(x\right)=\frac{\left({\omega }_0^2+\mu b_{eq}\right)x+\left(Q-\mu b_{eq}\right)x^3}{1+(1-x^2)\eta a_{eq}}.$$

(40)

Nonlinear equation (39) could be solved with the non-perturbative approach. The linearization configuration is performed in the form

$$\ddot{x}+\phi \left(\mathrm{\Omega }\right)\hspace{0.17em}\dot{x}+{\sigma }^2\left(\mathrm{\Omega }\right)x=0,$$

(41)

where $\phi \left(\mathrm{\Omega }\right)$ and $\sigma \left(\mathrm{\Omega }\right)$ are estimated using He’s frequency formula ⁴² as

$$\phi \left(\mathrm{\Omega }\right)={\frac{\left(1-x^2\right)\left(\mu a_{eq}\left(\mathrm{\Omega }\right)+\eta b_{eq}\left(\mathrm{\Omega }\right)\right)}{1+\left(1-x^2\right)\eta a_{eq}\left(\mathrm{\Omega }\right)}|}_{x=\frac{1}{2}A}=\frac{\left(1-\frac{1}{4}{A}^2\right)\left(\mu a_{eq}\left(\mathrm{\Omega }\right)+\eta b_{eq}\left(\mathrm{\Omega }\right)\right)}{1+\left(1-\frac{1}{4}{A}^2\right)\eta a_{eq}\left(\mathrm{\Omega }\right)},$$

(42)

$$\begin{array}{l}{\sigma }^2\left(\mathrm{\Omega }\right)={\frac{f\left(x\right)}{x}|}_{x=\frac{1}{2}A}=\underset{x\to \frac{1}{2}A}{\lim }\frac{\frac{d}{dx}\left[\left({\omega }_0^2+\mu b_{eq}\left(\mathrm{\Omega }\right)\right)x+\left(Q-\mu b_{eq}\left(\mathrm{\Omega }\right)\right)x^3\right]}{\frac{d}{dx}\left[x+\left(x-x^3\right)\eta a_{eq}\left(\mathrm{\Omega }\right)\right]}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{\left({\omega }_0^2+\mu b_{eq}\left(\mathrm{\Omega }\right)\right)+\frac{3}{4}\left(Q-\mu b_{eq}\left(\mathrm{\Omega }\right)\right)A^2}{1+\left(1-\frac{3}{4}{A}^2\right)\eta a_{eq}\left(\mathrm{\Omega }\right)},\end{array}$$

(43)

where $\mathrm{\Omega }$ is selected to be the total frequency of equation (38). Look for a solution to equation (41) in the form

$$y\left(t\right)=A{e}^{-\hspace{0.17em}\frac{1}{2}\phi \left(\mathrm{\Omega }\right)t}\cos \mathrm{\Omega }t$$

(44)

where the frequency could be expressed as

$${\mathrm{\Omega }}^2={\sigma }^2\left(\mathrm{\Omega }\right)-\frac{1}{4}{\phi }^2\left(\mathrm{\Omega }\right)$$

(45)

Employing equations (14) and (15) with equations (42) and (43) yields the frequency formula in the form

$${\mathrm{\Omega }}^2=\left(\frac{{\omega }^2+\mu \left(1-\frac{3}{4}{A}^2\right){\mathrm{\Omega }}^{\alpha }\cos \left(\frac{1}{2}\pi \alpha \right)}{\mathrm{\Omega }+\left(1-\frac{3}{4}{A}^2\right)\eta {\mathrm{\Omega }}^{\alpha }\sin \left(\frac{1}{2}\pi \alpha \right)}\right)\mathrm{\Omega }-\frac{{\left(1-\frac{1}{4}{A}^2\right)}^2{\mathrm{\Omega }}^{2\alpha }{\left(\mu \sin \left(\frac{1}{2}\pi \alpha \right)+\eta \mathrm{\Omega }\cos \left(\frac{1}{2}\pi \alpha \right)\right)}^2}{4{\left[\mathrm{\Omega }+\left(1-\frac{1}{4}{A}^2\right)\eta {\mathrm{\Omega }}^{\alpha }\sin \left(\frac{1}{2}\pi \alpha \right)\right]}^2}.$$

(46)

Numerical verification of the quasi-exact solution

To demonstrate the accuracy achieved by using the linearized approach to find the quasi-exact solution (44), a comparison with the numerical solution has been made and displayed in Figures 11–14. In Figure 11, the comparison is made for the case $\alpha =1$ and the solution is obtained in (44) with the numerical solution of equation (38). The system is selected for the behavior of the damping influence. The graph shows good agreement with the numerical solution. Consequently, the examination of the behavior when $0<\alpha <1$ appears in Figure 12. The variation of the parameter $\alpha$ has been plotted in case $\alpha =1$ , to be the reference for the comparison. This examination shows that as $\alpha$ decreases, the damping behavior should oscillate with a larger amplitude. When the parameter $\eta$ is selected to have a large value ( $\eta =1$ ) an opposite behavior for the amplitude of the damping force is found as shown in Figure 13. Decreasing $\alpha$ the amplitude is decreased making supports the damping behavior. The same behavior is observed before in Figure 10.

Figure 11.

The quasi-exact solution (44) compared with the numerical solution of equation (38) in the case of $\alpha =1$ the system of $\mu =0.3,\hspace{0.17em}\eta =0.4,Q=\hspace{0.17em}{\omega }_0=A=1$ .

Figure 12.

The variation of the parameter $\alpha$ of the quasi-exact solution (44) for a system having $\mu =0.3,\hspace{0.17em}\eta =0.01,Q=\hspace{0.17em}{\omega }_0=A=1$ .

Figure 13.

The variation of the parameter $\alpha$ of the quasi-exact solution (44) for the same system given in Figure 12 except that $\eta =1$ .

Figure 14.

The quasi-exact solution (41) compared with the numerical solution of equation (35) in the case of $\alpha =1$ the system of $\mu =-0.1,\hspace{0.17em}\eta =-0.01,\hspace{0.17em}Q=\hspace{0.17em}{\omega }_0=1,\hspace{0.17em}A=.1$ .

The case of the growing behavior that can occur for negative values of both the coefficient $\mu$ and $\eta$ has been studied $\alpha =1$ as displayed in Figure 14, and the variation $\alpha$ has been demonstrated in Figure 15. In the last case, the decreases in the parameter $\alpha$ lead to a decrease in the amplitude of the oscillation. This means that the decrease $\alpha$ plays a suppressing role in the growing behavior.

Figure 15.

The variation of the parameter $\alpha$ of the quasi-exact solution (44) for the system is given in Figure 14.

Conclusion

There are two main methods for research on fractional-order dynamical systems, namely, the numerical simulated and the approximately analytical perturbation methods. In the present work, a new approach is introduced to obtain an approximate solution without utilizing perturbation techniques to solve the fractional-order oscillation. In the new approach, two equivalent system parameters, that is, equivalent damping coefficient and equivalent stiffness coefficient, are defined, which could characterize the impacts of the fractional parameters on the conduct of the fractional-order Duffing oscillator. Also, the same approach has been used to perform the solution of the Van-der-Pol-Duffing oscillator with two fractional-orders is analytically investigated in this article. The confirmation between the traditional integer-order and the equivalent solution when the fraction factor has the value $\alpha =1$ is discussed numerically. The similarity between the numerical solution and the obtained analytical solution is confirmed the validation and satisfactory precision of the analytical solution. The technique and results may contribute to the understanding of similar fractional-order oscillators.