On the Stability of Duffing Type Fractional Differential Equation with Cubic Nonlinearity

Nnaji Daniel Ugochukwu¹, Makuochukwu Felix Oguagbaka^1*, Ezenwobodo Somkene Samuel², Nnaemeka Stanley Aguegboh³, Onyiaji Netochukwu Ebube¹

¹Department of Mathematics, University of Nigeria, Nsukka, Nigeria.
²Department of Mathematics, Nwafor Orizu College of Education, Nsugbe, Nigeria.
³Department of Mathematics, Veritas University, Abuja, Nigeria.

DOI: 10.4236/oalib.1106184                

Abstract

The purpose of this paper is to study the stability of nonlinear fractional Duffing equation where , by analysing the eigenvalues generated from the system of the given differential equation. Graphical results furthermore show the effect of the damping and nonlinear parameter on the system. Our contribution relies on its application to the choice of hard/soft spring in the mechanism of shock absorbers.

Keywords

Duffing Equation, Eigenvalues, Fractional Derivative, Hard Spring, Jacobian Matrix, Soft Spring, Stability

1. Introduction

The Duffing equation has received remarkable attention in the recent scientific research years. A typical example of its application to science, engineering and nature can be seen in [1] [2] where it modeled the non-linear mass-spring system, as well as the motion of a classical particle in a double-well potential [3]. Some Duffing equations as proposed by Correig in [4] as a model of microseism time series have been used in [5] [6] to model the prediction of earthquake occurrence. It was also used to model the transverse oscillation of non-linear beams in [7]. Many approaches have been developed for the study of the solution of the differential equation, many also have extended from integer to fractional case, seeking a better understanding (interpretation) on the differential equation.

Most problems that occur in real life are non-linear in nature and may not have analytic solutions by approximations or simulations and so the journey of trying to find an explicit solution may sometimes be complicated and sometimes impossible [8]. The Duffing’s equation is a second-order non-linear differential equation whose model application has been found useful in the model of damped and driven oscillators, modeling of brain [9], prediction of earthquake occurrences [10], signal processing [11] and crash analysis [12]. The integer type Duffing can be displayed by

$D^{2}x\left(t\right)+\delta D^{1}x\left(t\right)+f\left(t\mathrm{,}x\right)=g\left(t\right)$ (*)

where often time, $f\left(t\mathrm{,}x\right)=\rho x\left(t\right)+\mu x^3$ and $g\left(t\right)=\gamma sin\left(\omega t\right)$ with $\delta >0$. For $\rho >0$, $\mu >0$, the Equation (*) represents a “hard spring” and for $\rho >0$, $\mu <0$, the Equation (*) represents a “soft spring”.

On the other hand, fractional calculus is an old mathematical concept since the 17^th century. The fractional integral-differential operators are a generalization of integration and derivation to non-integer order (fractional) operators [13]. The idea of fractional calculus has been known since the development of the regular calculus, with the first reference probably being associated with Leibniz and L’Hopital in 1695. In the last two decades, fractional differential equations (FDEs) have been used to model various stable physical phenomena with anomalous decay, say that is not of exponential type [14]. Moreso, most differential systems used to describe physical phenomena are integer-order systems. With the development of fractional calculus, it has been found that the behaviour of many systems can be described using that fractional calculus system [15] [16] [17] [18] [19]. It is worth mentioning that many physical phenomena having memory and genetic characteristics can be described by using the fractional differential systems. The fractional type Duffing equation can be displayed by

$D^{\alpha }x\left(t\right)+\delta D^{\beta }x\left(t\right)+f\left(t\mathrm{,}x\right)=g\left(t\right)\mathrm{,}$ (**)

$f\left(t\mathrm{,}x\right)=\rho x\left(t\right)+\mu x^3$ and $g\left(t\right)=\gamma sin\left(\omega t\right)$ and $\delta \mathrm{,}\rho \mathrm{,}\mu \mathrm{,}\gamma \mathrm{,}\omega$ has its usual meaning.

Stability analysis is a central task in the study of fractional differential systems. The stability analysis of FDEs is more complex than that of classical differential equations, since fractional derivatives are nonlocal and have weakly singular kernels. Most of the known results on stability analysis of fractional differential systems concentrate on the stability of linear fractional differential systems. As stated by [20] in [21], Matignon who gave a well-known stability criterion for a

linear fractional autonomous differential system with constant coefficient matrix A. The criterion is that the stability is guaranteed if and only if the roots of the eigenfunction of the system lies outside the closed angular sector $\left|\arg \left(\lambda \left(A\right)\right)\right|<\frac{\alpha \pi }{2}$,

where $\alpha$ is the order of the fractional differential equation. The stability of fractional order nonlinear time delay systems for Caputo’s derivative was cited by [13] studied in [22], and two theorems for Mittage-Leffler stability of the fractional nonlinear time delay systems were proved.

In [13], stability analysis of Caputo fractional order nonlinear systems was studied using Bihari’s inequalities, Bellmann-Gronwalls inequality and the proof of comparison theorem for fractional order was proposed. In [23], the stability of fractional order nonlinear dynamic system was studied using Lyapunov direct method with the introductions of Mittage-Leffler stability and generalized Mittage-Leffler stability notions. Therefore, we derive our motivation of study from [8] [13] [23].

2. Preliminaries

Two major kinds of fractional derivatives, i.e., the Riemann-Liouville derivative and Caputo derivative, have often been used in fractional differential systems. We briefly state the two definitions of fractional derivatives as in [24]. But in this paper, we will adopt mainly the Caputo’s definition which is a modification of the Riemann-Liouville definition and has the advantage of dealing properly with initial value problem [13].

Definition 2.1. The Riemann-Liouville derivative with fractional order $\alpha \in {ℝ}_{+}$ of function $x\left(t\right)$ is defined below

${}_{RL}{D}_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)=\frac{{\text{d}}^m}{\text{d}{t}^m}{D}^{-\left(m-\alpha \right)}x\left(t\right)=\frac{1}{\Gamma \left(m-\alpha \right)}\frac{{\text{d}}^m}{\text{d}{t}^m}{\displaystyle {\int }_{t_0}^t}{\left(t-\tau \right)}^{m-\alpha -1}x\left(\tau \right)\text{d}\tau$ (a)

where $m-1<\alpha \le m\in {\mathbb{Z}}_{+}$.

Definition 2.2. The Caputo derivative with order $\alpha \in {ℝ}_{+}$ of function $x\left(t\right)$ is defined

${}_C{D}_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)=D^{-\left(m-\alpha \right)}\frac{{\text{d}}^m}{\text{d}{t}^m}x\left(t\right)=\frac{1}{\Gamma \left(m-\alpha \right)}{\displaystyle {\int }_{t_0}^t}{\left(t-\tau \right)}^{m-\alpha -1}{x}^{\left(m\right)}\left(\tau \right)\text{d}\tau$ (b)

where $m-1<\alpha \le m\in {\mathbb{Z}}_{+}$.

Definition 2.3. The Mittag-Leffler function is defined by

${\text{E}}_{\alpha }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{z^k}{\Gamma \left(\alpha k+1\right)}\mathrm{,}$ (c1)

where $\alpha >\mathrm{0,}z\in ℂ$. The two-parameter type Mittag-Leffler function is defined by

${\text{E}}_{\alpha \mathrm{,}\beta }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{z^k}{\Gamma \left(\alpha k+\beta \right)}\mathrm{,}$ (c2)

where $\alpha \mathrm{,}\beta >\mathrm{0,}z\in ℂ$.

Clearly, ${\text{E}}_{\alpha }\left(z\right)={\text{E}}_{\alpha \mathrm{,1}}\left(z\right)$ from the above equations.

Lemma 2.4. [24] If $0<\alpha <2$ and $\beta$ is an arbitrary complex number, then for an arbitrary integer $p\ge 1$ the following expansion holds:

${\text{E}}_{\alpha \mathrm{,}\beta }\left(z\right)=\frac{1}{\alpha }{z}^{\frac{1-\beta }{\alpha }}exp\left(z^{\frac{1}{\alpha }}\right)-{\displaystyle \underset{k=1}{\overset{p}{\sum }}}\frac{1}{\Gamma \left(\beta -\alpha k\right)}\cdot \frac{1}{z^k}+O\left(\frac{1}{{\left|z\right|}^{p+1}}\right)$ (d1)

with $\left|z\right|\to \infty$, $\left|arg\left(z\right)\right|\le \frac{\alpha \pi }{2}$ and

${\text{E}}_{\alpha \mathrm{,}\beta }\left(z\right)=-{\displaystyle \underset{k=1}{\overset{p}{\sum }}}\frac{1}{\Gamma \left(\beta -\alpha k\right)}\cdot \frac{1}{z^k}+O\left(\frac{1}{{\left|z\right|}^{p+1}}\right)$ (d2)

with $\left|z\right|\to \infty$, $\left|arg\left(z\right)\right|>\frac{\alpha \pi }{2}$. These results were proved in [20].

Definition 2.5. Consider the following fractional differential system

${}_{RL}{D}_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)=Ax\left(t\right)\mathrm{,}\left(0<\alpha <1\right)\mathrm{,}$ (1)

where $x\left(t\right)={\left(x_1\left(t\right)\mathrm{,}{x}_2\left(t\right)\mathrm{,}\cdots \mathrm{,}{x}_n\left(t\right)\right)}^{\text{T}}\in ℝ$, and $A\in {ℝ}^{n\times n}$.

We recall the following theorem from [25].

Theorem 2.6. The system (1) with initial value ${{}_{RL}{D}_{t_0\mathrm{,}t}^{\alpha -1}x\left(t\right)|}_{t=t_0}$, where $0<\alpha <1$ and $t_0=0$, is:

1) asymptotically stable if all the nonzero eigenvalues of A satisfy $\left|arg\left(spec\left(A\right)\right)\right|>\frac{\alpha \pi }{2}$, or A has K-multiple zero eigenvalues corresponding to a Jordan block:

$diag\left(J_1\mathrm{,}{J}_2\mathrm{,}\cdots \mathrm{,}{J}_i\right)$, where $J_l$ is a Jordan canonical form with order $n_l\times n_l$, ${\displaystyle \underset{l=1}{\overset{i}{\sum }}{n}_l}=k$, and $n_l\alpha <1$ for each $1\le l\le i$.

2) stable if all the non-zero eigenvalues of A satisfy $\left|arg\left(spec\left(A\right)\right)\right|\ge \frac{\alpha \pi }{2}$ and the critical eigenvalues satisfying $\left|arg\left(spec\left(A\right)\right)\right|=\frac{\alpha \pi }{2}$ have the same algebraic and geometric multiplicities, or A has k-multiple zero eigenvalues corresponding to a Jordan block matrix $diag\left(J_1\mathrm{,}{J}_2\mathrm{,}\cdots \mathrm{,}{J}_i\right)$ where $J_l$ is a Jordan canonical form with order $n_l\times n_l$, ${\displaystyle \underset{l=1}{\overset{i}{\sum }}{n}_l}=k$, and $n_l\alpha <1$ for each $1\le l\le i$.

Lemma 2.7. [14] Consider the system of differential equation

$D_{t_0\mathrm{,}t}^{\bar{\alpha }}x\left(t\right)=f\left(t\mathrm{,}x\right)$ (2)

with suitable initial values $x_k={\left[x_{k_1}\mathrm{,}{x}_{k_2}\mathrm{,}\cdots \mathrm{,}{x}_{k_n}\right]}^{\text{T}}\in {ℝ}^n$ $\left(k=\mathrm{0,1,}\cdots \mathrm{,}m-1\right)$, where $x\left(t\right)={\left[x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_n\right]}^{\text{T}}\in {ℝ}^n$, $\bar{\alpha }={\left[{\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}\cdots \mathrm{,}{\alpha }_n\right]}^{\text{T}}$, $m-1<{\alpha }_i<m\in {\mathbb{Z}}_{+}$ $\left(i=\mathrm{1,2,}\cdots \mathrm{,}n\right)$ and $D_{t_0\mathrm{,}t}^{{\alpha }_i}$ denotes either ${}_C{D}_{t_0\mathrm{,}t}^{{\alpha }_i}$ or ${}_{RL}{D}_{t_0\mathrm{,}t}^{{\alpha }_i}$.

In particular, if ${\alpha }_1={\alpha }_2=\cdots ={\alpha }_n=\alpha$, then (2) can be written as

$D_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)=f\left(t\mathrm{,}x\right)$ (3)

Then the following definitions are associated with the stability problem:

Definition 2.8. [23] The constant vector $x_{eq}$ is an equilibrium point of the fractional differential system (2), if and only if $f\left(t\mathrm{,}{x}_{eq}\right)={{}_{RL}{D}_{t_0\mathrm{,}t}^{\bar{\alpha }}x\left(t\right)|}_{x\left(t\right)=x_{eq}},\forall t>t_0$.

without loss of generality, let the equilibrium point be $x_{eq}=0$, we introduce the following definition [14].

Definition 2.9. The zero solution of fractional differential system (2) is said to be stable if for any initial value $x_k={\left[x_{k_1}\mathrm{,}{x}_{k_2}\mathrm{,}\cdots \mathrm{,}{x}_{k_n}\right]}^{\text{T}}\in {ℝ}^n$ $\left(k=\mathrm{0,1,}\cdots \mathrm{,}m-1\right)$, there exists $\epsilon >0$ such that any solution $x\left(t\right)$ of (2) satisfies $\Vert x\left(t\right)\Vert <\epsilon$ $\forall t>t_0$. The zero solution is said to be asymptotically stable if in addition to being stable, $\Vert x\left(t\right)\Vert \to 0$ as $t\to +\infty$.

Alidousti et al. [20] proposed the stability of linear fractional differential equations as follows:

Consider the following fractional differential system defined in the Caputo sense

${}_C{D}_{0,t}^{\alpha }x\left(t\right)=Ax\left(t\right)+f\left(t,x\left(t\right)\right),\left(0<\alpha <1\right)$ (4)

under the initial condition $x\left(0\right)=x_0$, where $x\left(t\right)={\left(x_1\left(t\right)\mathrm{,}{x}_2\left(t\right)\mathrm{,}\cdots \mathrm{,}{x}_n\left(t\right)\right)}^{\text{T}}\in {ℝ}^n$ and $A\in {ℝ}^{n\times n}$, $f\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\infty \right)\times {ℝ}^n\to {ℝ}^n$. We can get the solution of (4), by using the Laplace and inverse Laplace transforms, as

$x\left(t\right)={\text{E}}_{\alpha }\left(A{t}^{\alpha }\right)x_0+{\displaystyle {\int }_0^t}{\left(t-\tau \right)}^{\alpha -1}{\text{E}}_{\alpha \mathrm{,}\alpha }\left(A{\left(t-\tau \right)}^{\alpha }\right)f\left(\tau \mathrm{,}x\left(\tau \right)\right)\text{d}\tau$ (5)

and then the following stability results where established.

Theorem 2.10. Suppose f is a continuous vector function for which there exists $p>\frac{1}{\alpha }$ such that $\Vert f\left(\mathrm{.,}x\left(t\right)\right)\Vert \in L^p\left({ℝ}^{+}\right)$. Then the system (4) is stable if all the eigenvalues of A satisfy

$\left|arg\left(eig\left(A\right)\right)\right|>\frac{\alpha \pi }{2}$ (6)

especially if $\Vert f\left(\mathrm{.,}x\left(t\right)\right)\Vert \in L^p\left({ℝ}^{+}\right)$ the stability holds.

The prove of the above Theorem is in [20].

From the concept of the theorem stated in ( [26], p. 169), the following corollary was stated.

Corollary 2.11. Consider the equation

${}^C{D}_{t_0\mathrm{,}t}^{N+\alpha }x\left(t\right)=f\left(t\mathrm{,}x\left(t\right)\mathrm{,}{D}_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)\right)\mathrm{,}0<\alpha \le \mathrm{1,}$ (7a)

with the initial condition

$x^{\left(j\right)}\left(0\right)=x_0^{\left(j\right)}\mathrm{,}j=\mathrm{0,1,}\cdots \mathrm{,}\left[N\right]-1$ (7b)

where $f\left(t\mathrm{,}x\left(t\right)\mathrm{,}{D}_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)\right)=-\delta D_{t_0\mathrm{,}t}^{\alpha }x\left(t\right)-\rho x\left(t\right)-\mu x^3\left(t\right)$.

the augmented equivalent system can be written as in [27].

$(\begin{array}{l}{D}^{\alpha }x\left(t\right)=D^{\alpha }{x}_0\left(t\right)=x_1\left(t\right)\\ D^{2\alpha }x\left(t\right)=D^{\alpha }\left(D^{\alpha }{x}_0\left(t\right)\right)=D^{\alpha }{x}_1\left(t\right)=x_2\left(t\right)\\ D^{3\alpha }x\left(t\right)=D^{\alpha }\left(D^{2\alpha }{x}_0\left(t\right)\right)=D^{\alpha }{x}_2\left(t\right)=x_3\left(t\right)\\ ⋮\\ D^{\left(N-1\right)\alpha }x\left(t\right)=D^{\alpha }\left(D^{\left(N-2\right)\alpha }{x}_0\left(t\right)\right)=D^{\alpha }{x}_{N-2}\left(t\right)=x_{N-1}\left(t\right)\\ D^{N\alpha }x\left(t\right)=D^{\alpha }\left(D^{\left(N-1\right)\alpha }{x}_0\left(t\right)\right)=D^{\alpha }{x}_{N-1}\left(t\right)=-\delta x_1\left(t\right)-\rho x_0\left(t\right)-\mu x_0^3\left(t\right)\end{array}$ (8)

which can further be simplified to

$(\begin{array}{l}{D}^{\alpha }{x}_0\left(t\right)=x_1\left(t\right)\\ D^{\alpha }{x}_1\left(t\right)=x_2\left(t\right)\\ D^{\alpha }{x}_2\left(t\right)=x_3\left(t\right)\\ ⋮\\ D^{\alpha }{x}_{N-2}\left(t\right)=x_{N-1}\left(t\right)\\ D^{\alpha }{x}_{N-1}\left(t\right)=-\delta x_1-\rho x_0-\mu x_0^3\end{array}$ (9)

where the function x also satisfies the initial condition (7b).

The eigenvalues of system (9), ${\lambda }_i$ for $i=\mathrm{1,2,}\cdots \mathrm{,}n$ is written as follows

$D^{\alpha }\left(\begin{array}{c}{x}_0\\ x_1\\ x_2\\ ⋮\\ x_{N-1}\end{array}\right)=\left|\begin{array}{ccccc}0-\lambda & 1& 0& 0& \cdots \\ 0& 0-\lambda & 1& 0& \cdots \\ 0& 0& 0-\lambda & 1& \cdots \\ ⋮& \ddots & \ddots & \ddots & ⋮\\ -\rho -3\mu x_0^2& -\delta & \cdots & \cdots & 0-\lambda \end{array}\right|=0$ (10)

which gives

${\lambda }^n+\delta \lambda +\left(\rho +3\mu x_0^2\right)=0$ (11)

We therefore solve for $\lambda$ at different values of $n\in \mathbb{N}$, $\delta \mathrm{,}\rho \in {ℝ}_{+}$ and $\mu \in ℝ-\left\{0\right\}$.

Note: $x_i$ ’s can be evaluated at equilibrium point(s) i.e.

$(\begin{array}{l}{D}^{\alpha }{x}_0\left(t\right)=x_1\left(t\right)=0\\ D^{\alpha }{x}_1\left(t\right)=x_2\left(t\right)=0\\ D^{\alpha }{x}_2\left(t\right)=x_3\left(t\right)=0\\ ⋮\\ D^{\alpha }{x}_{N-2}\left(t\right)=x_{N-1}\left(t\right)=0\\ D^{\alpha }{x}_{N-1}\left(t\right)=-\delta x_1-\rho x_0-\mu x_0^3=0\end{array}$ (12)

which results in:

・ Trivial equilibrium point

$x_i=\left(x_0,x_1,\cdots ,x_{N-1}\right)=\left(0,0,\cdots ,0\right)$

・ Other equilibrium points

$x_i=\left(x_0\mathrm{,}{x}_1\mathrm{,}\cdots \mathrm{,}{x}_{N-1}\right)=\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,}\cdots \mathrm{,0}\right)\mathrm{,}\mu \ne 0$

Then Equation (11) will be evaluated and analysed at each of these equilibrium points.

3. Illustrative Examples

Example 3.1. Consider the homogeneous (unforced) non-linear fractional Duffing oscillator given by

$D^{\frac{3}{2}}x\left(t\right)+\delta D^{\frac{1}{2}}x\left(t\right)+\rho x\left(t\right)+\mu x^3=0$ (13)

whose equivalent system is gotten by letting $x=x_0$, so that

$D^{\frac{1}{2}}{x}_0=x_1$

$D^{\frac{1}{2}}{x}_1=x_2$ (14)

$D^{\frac{1}{2}}{x}_2=-\delta x_1-\rho x_0-\mu x_0^3$

The Jacobian matrix for Equation (14) is given as follows

$J=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -\rho -3\mu x_0^2& -\delta & 0\end{array}\right)$ (15)

The eigenvalues, ${\lambda }_i$ for $i=\mathrm{0,1,2}$ are obtained as

$\left|J-\lambda I\right|=\left|\begin{array}{ccc}-\lambda & 1& 0\\ 0& -\lambda & 1\\ -\rho -3\mu x_0^2& -\delta & -\lambda \end{array}\right|=0$ (16)

which gives

${\lambda }^3+\delta \lambda +\left(\rho +3\mu x_0^2\right)=0$ (17)

At equilibrium, $D^{\frac{1}{2}}{x}_i=0$ for $i=\mathrm{0,1,2}$. Therefore from system (14)

$x_1=0$ (18)

$x_2=0$ (19)

$-\delta x_1-\rho x_0-\mu x_0^3=0$ (20)

From (20)

$-\delta x_1-x_0\left(\rho +\mu x_0^2\right)=0$

$\text{at}{x}_1=0\Rightarrow x_0\left(\rho +\mu x_0^2\right)=0$

$x_0=0,\rho +\mu x_0^2=0$

$x_0=0,x_0^2=\frac{-\rho }{\mu }$

$x_0=0,x_0=\pm \sqrt{\frac{-\rho }{\mu }}$

The possible equilibrium points are:

$\begin{array}{c}\left(x_0,x_1,x_2\right)=\left(0,0,0\right),\left(\sqrt{\frac{-\rho }{\mu }},0,0\right),\left(-\sqrt{\frac{-\rho }{\mu }},0,0\right)\\ =\left(\mathrm{0,0,0}\right)\mathrm{,}\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,0}\right)\end{array}$ (21)

We evaluate Equation (17) for each of the equilibrium points in (21).

Case 1

For the trivial equilibrium point,

${\left|J-\lambda I\right|}_{\left(\mathrm{0,0,0}\right)}={\lambda }^3+\delta \lambda +\left(\rho +3\mu x_0^2\right)=0$

${\lambda }^3+\delta \lambda +\rho =0$ (22)

to investigate the stability of the analytic solution, we therefore solve for $\lambda$ for different values of $\delta$ and $\rho$ (where $\alpha =\frac{1}{2}$ and $\frac{\alpha \pi }{2}\approx 0.7854$ ):

Table 1 below shows that for any positive values of $\delta$ and $\rho$, the equilibrium point $\left(x_0\mathrm{,}{x}_1\mathrm{,}{x}_2\right)=\left(\mathrm{0,0,0}\right)$ of the system (13) by Theorem (2.6) is asymptotically stable. Garrappa also asserted in her work that for negative values of $\rho$, the system (13) is unstable.

Case 2

For

${\left|J-\lambda I\right|}_{\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,0}\right)}={\lambda }^3+\delta \lambda +\left(\rho +3\mu x_0^2\right)=0$

${\lambda }^3+\delta \lambda +\left[\rho +3\mu {\left(\pm \sqrt{\frac{-\rho }{\mu }}\right)}^2\right]=0$

${\lambda }^3+\delta \lambda +\left[\rho +3\mu {\left(\sqrt{\frac{-\rho }{\mu }}\right)}^2\right]=0$ (23)

・ If $\mu >0$ (hard spring), (23) becomes

${\lambda }^3+\delta \lambda -2\rho =0$ (24)

・ also for $\mu <0$ (soft spring), (23) becomes

${\lambda }^3+\delta \lambda -2\rho =0$ (25)

As shown in Table 2 below, for positive value (s) of $\delta$ and $\rho$ the equilibrium point $\left(x_0\mathrm{,}{x}_1\mathrm{,}{x}_2\right)=\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,0}\right)$ by Theorem (2.6) is unstable since $\left|arg\lambda \right|\cancel{>}\frac{\pi }{4}\mathrm{,}\forall {\lambda }_i$.

Example 3.2. Consider the homogeneous non-linear fractional Duffing oscillator given by

$D^{\frac{4}{3}}x\left(t\right)+\delta D^{\frac{1}{3}}x\left(t\right)+\rho x\left(t\right)+\mu x^3=0$ (26)

whose equivalent system is

$D^{\frac{1}{3}}{x}_0=x_1$

$D^{\frac{1}{3}}{x}_1=x_2$ (27)

$D^{\frac{1}{3}}{x}_2=x_3$

Table 1. Showing the stability results for Equation (22) for different values of $\delta$ and $\rho$.

Table 2. Showing the stability results for Equations (24, 25) for different values of $\delta$ and $\rho$.

$D^{\frac{1}{3}}{x}_3=-\delta x_1-\rho x_0-\mu x_0^3$

The Jacobian matrix for Equation (27) is given as follows

$J=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ -\rho -3\mu x_0^2& -\delta & 0& 0\end{array}\right)$ (28)

The eigenvalues, ${\lambda }_i$ for $i=\mathrm{0,1,2,3}$ are obtained as

$\left|J-\lambda I\right|=\left|\begin{array}{cccc}-\lambda & 1& 0& 0\\ 0& -\lambda & 1& 0\\ 0& 0& -\lambda & 1\\ -\rho -3\mu x_0^2& -\delta & 0& -\lambda \end{array}\right|=0$ (29)

which gives

${\lambda }^4+\delta \lambda +\left(\rho +3\mu x_0^2\right)=0$ (30)

As usual if proper arithmetic is done, we note that the equilibrium points of the system (27) are:

$\left(x_0\mathrm{,}{x}_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)=\left(\mathrm{0,0,0,0}\right)\mathrm{,}\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,0,0}\right)$ (31)

Case 1

For the trivial equilibrium point $\left(x_0\mathrm{,}{x}_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)=\left(\mathrm{0,0,0,0}\right)$, we have

${\lambda }^4+\delta \lambda +\rho =0$ (32)

to investigate the stability of the analytic solution, we solve for $\lambda$ for different values of $\delta$ and $\rho$ (where $\alpha =\frac{1}{3}$ and $\frac{\alpha \pi }{2}\approx 0.5236$ ).

From Table 3 below, it can be seen that the stability of the equilibrium point $\left(x_0\mathrm{,}{x}_1\mathrm{,}{x}_2\mathrm{,}{x}_3\right)=\left(\mathrm{0,0,0,0}\right)$ of the analytic solution was inconclusive (indecisive) since the stability/instability of the system does not depend on any specific value of $\delta$ or $\rho$.

Case 2

For the remaining equilibrium points, we compute:

${\left|J-\lambda I\right|}_{\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,0,0}\right)}={\lambda }^4+\delta \lambda +\left(\rho +3\mu x_0^2\right)=0$

${\lambda }^4+\delta \lambda +\left[\rho +3\mu {\left(\pm \sqrt{\frac{-\rho }{\mu }}\right)}^2\right]=0$

${\lambda }^4+\delta \lambda +\left[\rho +3\mu {\left(\sqrt{\frac{-\rho }{\mu }}\right)}^2\right]=0$ (33)

・ If $\mu >0$ (hard spring), we have: ${\lambda }^4+\delta \lambda -2\rho =0$.

・ If $\mu <0$ (soft spring), we also have: ${\lambda }^4+\delta \lambda -2\rho =0$.

It can be observed that the soft and hard spring appear to respond the same way. Then we investigate the stability for different values of $\delta$ and $\rho$ (Table 4).

Therefore, the equilibrium points $\left(\pm \sqrt{\frac{-\rho }{\mu }}\mathrm{,0,0,0}\right)$ for system (26) is unstable for both hard and soft spring for all positive value (s) of $\delta$ and $\rho$ (Figures 1-3).

Table 3. Showing the stability results for Equation (32) for different values of $\delta$ and $\rho$.

Table 4. Showing the stability results for Equation (33) for different values of $\delta$ and $\rho$.

(a) $\alpha =1.25$, $\beta =0.25$, $\rho =1$, $\mu =1$, $c=1$ (b) $\alpha =\frac{4}{3}$, $\beta =\frac{1}{3}$, $\rho =1$, $\mu =1$, $c=1$

Figure 1. FDE plot for displacement (x) against time (t), for increasing $\alpha$ and $\beta$ with variation in damping parameter.

(a) $\alpha =\frac{3}{2}$, $\beta =\frac{1}{2}$, $\rho =1$, $\mu =1$, $c=1$ (b) $\alpha =1.85$, $\beta =0.85$, $\rho =1$, $\mu =1$, $c=1$

Figure 2. FDE plot for displacement (x) against time (t), as $\alpha \to 2$, $\beta \to 1$ with variation in damping parameter.

(a) $\alpha =\frac{3}{2}$, $\beta =\frac{1}{2}$, $\rho =1$, $\mu =1$, $c=1$ (b) $\alpha =\frac{8}{5}$, $\beta =\frac{3}{5}$, $\rho =1$, $\mu =1$, $c=1$

Figure 3. FDE plot for displacement (x) against time (t), with variation in the non-linear term.

4. Conclusion

In this paper, the survey on the stability of Duffing type FDEs with cubic nonlinearity was studied. The concept of spring motion in the shock absorber of a vehicle in motion was adopted to analyse the stability of the equations used. In the above concept, stability can be discussed as follows:

・ Stability, in this case, means that the spring returns back to normal (almost) immediately after disturbance (disturbance due to bumps or stip slope) as $t\to \infty$.

・ Instability, in this case, means the spring may not quickly return back to normal after a disturbance, which in turn may lead to the vehicle moving off track as $t\to \infty$.

Acknowledgements

We wish to thank Prof. M. O. Oyesanya and my co-authors for their kind suggestions towards the success of this work.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hassan, A. (1973) Perturbation Method. John Wiley & Sons Inc., Hoboken.
[2]	Amazigo, J.C. (2011) Postgraduate Lecture Notes on Perturbation Methods. Department of Mathematics, University of Nigeria, Nsukka.
[3]	Rand, R.H. (2003) Lecture Notes on Nonlinear Vibration. http://www.tam.cornell.edu/randdocs/nlvibe54.pdf
[4]	Correig, A.M. and Urquizu, M. (2002) Some Dynamical Aspects of Microseism Time Series. Geophysical Journal International, 149, 589-598. https://doi.org/10.1046/j.1365-246X.2002.01602.x
[5]	Oyesanya, M.O. (2008) Duffing Oscillator as a Model for Predicting Earthquake Occurrence. Journal of Nigerian Association of Mathematical Physics, 12, 133-142. https://doi.org/10.4314/jonamp.v12i1.45495
[6]	Oyesanya, M.O. and Nwamba, J.I. (2013) Stability Analysis of Damped Cubic-Quintic Duffing Oscillator. World Journal of Mechanics (Scientific Research), 3, 43-57. https://doi.org/10.4236/wjm.2013.31003
[7]	Sedighi, H.M., Shirazi, K.H. and Zare, J. (2012) An Analytic Solution of Transversal Oscillation of Quintic Nonlinear Beam with Homotopy Analysis Method. International Journal of Nonlinear Mechanics, 47, 777-784. https://doi.org/10.1016/j.ijnonlinmec.2012.04.008
[8]	Eze, E.O., Obasi, U.E. and Agwu, E.U. (2019) Stability Analysis of Periodic Solution of Some Duffing’s Equations. Open Journal of Applied Sciences, 9, 198-214. https://doi.org/10.4236/ojapps.2019.94017
[9]	Zeeman, E. (1976) Duffing Equation in Brain Modelling. Bulletin Institute of Mathematics and Its Applications, 12, 207-214.
[10]	Zeeman, E. (2008) Duffing Oscillator as Model for Predicting Earthquake Occurrence I. Journal of Nigerian Association of Mathematical Physics, 12, 133-142. https://doi.org/10.4314/jonamp.v12i1.45495
[11]	Wang, G., Zheng, W. and He, S. (2002) Estimation of Amplitude and Phase of a Weak Signal by Using the Property of Sensitive Dependence on Initial Condition of a Non-Linear Oscillator. Signal Processing, 82, 103-115. https://doi.org/10.1016/S0165-1684(01)00166-9
[12]	Yang, L. and Li, Y.K. (2014) Existence and Global Exponential Stability of Almost Periodic Solution for a Class of Delay Duffing Equation on Time Scale. Abstract and Applied Analysis, 2014, Article ID: 857161. https://doi.org/10.1155/2014/857161
[13]	Hadi, D., Dumitru, B. and Jalil, S. (2012) Stability Analysis of Caputo Fractional-Order Nonlinear Systems Revisited. Nonlinear Dynamics, 67, 2433-2439. https://doi.org/10.1007/s11071-011-0157-5
[14]	Li, C.P. and Zhang, F.R. (2011) A Survey on the Stability of Fractional Differential Equations. The European Physical Journal Special Topics, 193, 27-47. https://doi.org/10.1140/epjst/e2011-01379-1
[15]	Bagley, R.L. and Calico, R.A. (1991) Fractional Order State Equations for the Control of Viscoelastically Damped Structures. Journal of Guidance, Control, and Dynamics, 14, 304. https://doi.org/10.2514/3.20641
[16]	Sun, H.H., Abdelwahab, A.A. and Onaral, B. (1984) Linear Approximation of Transfer Function with a Pole of Fractional Power. IEEE Transactions on Automatic Control, 29, 441-444. https://doi.org/10.1109/TAC.1984.1103551
[17]	Ichisea, M., Nagayanagia, Y. and Kojima, T.J. (1971) An Analog Simulation of Non-Integer Order Transfer Functions for Analysis of Electrode Processes. Journal of Electroanalytical Chemistry, 33, 253.
[18]	Laskin, N. (2000) Fractional Market Dynamics. Physica A: Statistical Mechanics and Its Applications, 287, 482-492. https://doi.org/10.1016/S0378-4371(00)00387-3
[19]	Kusnezov, D., Bulgac, A. and Dang, G.D. (1999) Quantum Lévy Processes and Fractional Kinetics. Physical Review Letters, 82, 1136. https://doi.org/10.1103/PhysRevLett.82.1136
[20]	Alidousti, J., Ghaziani, R.K. and Bayati, A.E. (2017) Stability Analysis of Nonlinear Fractional Differential Order Systems with Caputo and Riemann-Liouville Derivatives. Turkish Journal of Mathematics, 41, 1260-1278. https://doi.org/10.3906/mat-1510-5
[21]	Matignon, D. (1996) Stability Results for Fractional Differential Equations with Applications to Control Processing. Computational Engineering in Systems and Application Multiconference, Vol. 2, 963-968.
[22]	Sadati, S.J., Baleanu, D., Ranjbar, A., Ghaderi, R. and Abdeljawad, T. (2010) Mittag-Leffler Stability Theorem for Fractional Nonlinear Systems with Delay. Abstract and Applied Analysis, 2010, Article ID: 108651. https://doi.org/10.1155/2011/213485
[23]	Li, Y., Chen, Y.Q. and Podlubny, I. (2009) Mittag-Leffler Stability of Fractional Order Nonlinear Dynamics Systems. Automatica, 45, 1965-1969. https://doi.org/10.1016/j.automatica.2009.04.003
[24]	Podlubny, I. (1999) Fractional Differential Equations. Academic Press, New York.
[25]	Qian, D., Li, C., Agarwal, R.P. and Wong, P.J. (2010) Stability Analysis of Fractional Differential System with Riemann-Liouville Derivative. Mathematical and Computer Modelling, 52, 862-874. https://doi.org/10.1016/j.mcm.2010.05.016
[26]	Diethelm, K. (2010) The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type. Springer, Berlin. https://doi.org/10.1007/978-3-642-14574-2
[27]	Garrappa, R. (2010) On Linear Stability of Predictor-Corrector Algorithms for Fractional Differential Equations. International Journal of Computer Mathematics, 87, 2281-2290. https://doi.org/10.1080/00207160802624331