Bateman–Feshbach Tikochinsky and Caldirola–Kanai Oscillators with New Fractional Differentiation

Abstract

In this work, the study of the fractional behavior of the Bateman–Feshbach–Tikochinsky and Caldirola–Kanai oscillators by using different fractional derivatives is presented. We obtained the Euler–Lagrange and the Hamiltonian formalisms in order to represent the dynamic models based on the Liouville–Caputo, Caputo–Fabrizio–Caputo and the new fractional derivative based on the Mittag–Leffler kernel with arbitrary order α. Simulation results are presented in order to show the fractional behavior of the oscillators, and the classical behavior is recovered when α is equal to 1.

1. Introduction

Several phenomenological models of dissipative systems have been proposed, such as the Bateman–Feshbach–Tikochinsky (BFT) or Caldirola–Kanai (CK) oscillators, the first model consists of a damped and an amplified oscillator, and this one-dimensional system exhibits an exponentially increasing mass with a Lagrangian given by Bateman [1,2,3,4,5]. Both quantum damped oscillators have been studied as a model to understand dissipation in quantum theory [6]. Bateman suggested the time-dependent Hamiltonian [2] and Caldirola the time dependent Hamiltonian to describe damped oscillations [4]. The Caldirola–Kanai oscillator is an open system whose parameters such as mass and frequency are all time dependent, while the Bateman–Feshbach–Tikochinsky oscillator is a closed system whose total energy is conserved and the dissipated energy from the damped oscillator is transferred to amplified one [7,8]. The fractional Hamiltonians are non-local and they are associated with dissipative systems [8]. There are few definitions of operators with fractional order, the Liouville–Caputo fractional derivative involving a kernel with singularity, and this definition is based on the power law and present singularity at the origin [9]. Recently, in order to solve the problem of singularity at the origin, Caputo and Fabrizio used the exponential decay law to construct a derivative with no singularity; however, the used kernel was local [10,11,12,13,14,15,16,17,18]. Thus, Atangana and Baleanu used the generalized Mittag–Leffler function to construct a derivative with no-singular and non-local kernel [19,20,21,22]. In this paper, we obtain alternative representations of the BFT and CK oscillators by using the Liouville–Caputo, Caputo–Fabrizio–Caputo and the new fractional derivative based in Mittag–Leffler kernel with arbitrary order α. Numerical solutions are based in a Crank–Nicholson scheme.

2. Fractional Operators

The Adams method is a multi-step method, and this method uses the information of all the previous values, $y_i$, $y_{i-1}$, $y_{i-m+1}$, in order to calculate $y_{i+1}$. This is the difference between the Adams method and the single-step methods, such as the Heun, Taylor and Runge–Kutta numerical schemes, which use only the last value to calculate the next one. There are two types of Adams methods, the Adams–Bashforth and the Adams–Moulton. The combination of these methods results in the predictor–corrector Adams–Bashforth–Moulton Method [23,24,25,26].

The generalization of this method for any order of derivative is called the fractional Adams–Bashforth Method [23]

$${}_0^C{\mathcal{D}}_t^{\alpha }f\left(t\right)=g(t,f\left(t\right)),f^w\left(0\right)=f_0^w,w=0,1,...,n-1,$$

(1)

where $\alpha >0$ and ${}_0^C{D}_t^{\alpha }$ is the Liouville–Caputo operator

$${}_0^C{\mathcal{D}}_t^{\alpha }f\left(t\right)=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^t\frac{f^{\left(n\right)}\left(\eta \right)}{{(t-\eta )}^{\alpha -n+1}}d\eta .$$

(2)

Equation (1) satisfies the following Volterra integral equation

$$f\left(t\right)=\sum _{w=0}^{n-1}{f}_0^{\left(w\right)}{\displaystyle \frac{t^w}{w!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}\underset{0}{\overset{t}{\int }}{(t-\kappa )}^{\alpha -1}g(\kappa ,f\left(\kappa \right))d\kappa \kappa ,t<T.$$

(3)

The fractional Adams method to solve (1) has been studied firstly by Diethelm, Ford and Freed [24], and this solution scheme is derived as follows:

$$\begin{array}{c}{f}_{w+1}^P={\displaystyle \sum _{j=0}^{n-1}}{\displaystyle \frac{t_{w+1}^j}{j!}}{f}_0^{\left(j\right)}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\displaystyle \sum _{j=0}^w}{b}_{j,w+1}g(t_j,f_j),\hfill \\ f_{w+1}={\displaystyle \sum _{j=0}^{n-1}}{\displaystyle \frac{t_{w+1}^j}{j!}}{f}_0^{\left(j\right)}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}\left({\displaystyle \sum _{j=0}^w}{a}_{j,w+1}g(t_j,f_j)+a_{w+1,w+1}g(t_{w+1},f_{k+1}^P)\right).\hfill \end{array}$$

(4)

The fractional operator proposed by Caputo and Fabrizio in Liouville–Caputo sense (CFC) is expressed as follows [10]:

$${}_0^{CFC}{\mathcal{D}}_t^{\alpha }f\left(t\right)=\frac{(2-\alpha )B\left(\alpha \right)}{2(1-\alpha )}\underset{0}{\overset{t}{\int }}exp\left(\frac{-\alpha }{1-\alpha }(t-\varsigma )\right)f^{\left(n\right)}\left(\varsigma \right)d\varsigma ,$$

(5)

where $B\left(\alpha \right)=B\left(0\right)=B\left(1\right)=1$ (is a normalization function). In this sense, the Laplace transform is given by

$$\begin{array}{cc}\hfill \mathcal{L}\left[{}_0^{CFC}{D}_t^{\left(n+\alpha \right)}f\left(t\right)\right]\left(s\right)& =\frac{s^{n+1}\mathcal{L}\left[f\left(t\right)\right]-s^{n}f\left(0\right)-s^{n-1}{f}^{\prime }\left(0\right)\dots -f^{\left(n\right)}\left(0\right)}{s+\alpha \left(1-s\right)}.\hfill \end{array}$$

(6)

The fractional derivative based in Mittag–Leffler kernel (Atangana–Baleanu fractional operator in Liouville–Caputo sense, ABC) is given as

$${}_0^{ABC}{\mathcal{D}}_t^{\alpha }f\left(t\right)=\frac{B\left(\alpha \right)}{1-\alpha }{\int }_0^t\dot{f}\left(\theta \right)E_{\alpha }\left[-\alpha \frac{{(t-\theta )}^{\alpha }}{1-\alpha }\right]d\theta ,$$

(7)

where $E_{\alpha }$ is a Mittag–Leffler function [19]. The fractional integral is defined as

$${}_a^{AB}{I}_t^{\alpha }f\left(t\right)=\frac{1-\alpha }{B\left(\alpha \right)}f\left(t\right)+\frac{\alpha }{B\left(\alpha \right)\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^{t}f\left(\varsigma \right){(t-\varsigma )}^{\alpha -1}d\varsigma .$$

(8)

The Laplace transform of (7) produces

$$\begin{array}{c}\hfill \mathcal{L}\left[{}_0^{ABC}{\mathcal{D}}_t^{\alpha }f\left(t\right)\right]\left(s\right)=\frac{B\left(\alpha \right)}{1-\alpha }\frac{s^{\alpha }\mathcal{L}\left[f\left(t\right)\right]\left(s\right)-s^{\alpha -1}f\left(0\right)}{s^{\alpha }+\frac{\alpha }{1-\alpha }}.\end{array}$$

(9)

3. Applications

3.1. Bateman–Feshbach–Tikochinsky Oscillator

The classical Lagrangian of the BFT oscillator is given by

$$L=m{\dot{q}}_1{\dot{q}}_2+\rho (q_1{\dot{q}}_2-{\dot{q}}_1{q}_2)-K{q}_1{q}_2,$$

(10)

where $q_1$ is the damped harmonic oscillator coordinate and $q_2$ corresponds to the time-reversed counterpart, and the parameters m, ρ, K are time independent.

The fractional Lagrangian (10) is given by

$$L^F=m{}_{a}D_t^{\alpha }{q}_1{}_{a}D_t^{\alpha }{q}_2+\rho (q_1{}_{a}D_t^{\alpha }{q}_2-{}_{a}D_t^{\alpha }{q}_1{q}_2)-K{q}_1{q}_2,$$

(11)

and the Lagrange model of fractional order is

$$\begin{array}{c}m{}_{a}D_t^{\alpha }{}_{a}D_t^{\alpha }{q}_1+\rho {}_{a}D_t^{\alpha }{q}_1+K{q}_1=0,\hfill \\ m{}_{a}D_t^{\alpha }{}_{a}D_t^{\alpha }{q}_2-\rho {}_{a}D_t^{\alpha }{q}_2+K{q}_2=0.\hfill \end{array}$$

(12)

Now, we can get the generalized momentum as follows:

$$p_i={\displaystyle \frac{\partial L^F}{\partial {}_{a}D_t^{\alpha }{q}_i}},$$

(13)

where $L^F$ is the Lagrangian of fractional order and $i=1,2$.

The two generalized momentums are given by

$$\begin{array}{c}{p}_1={\displaystyle \frac{\partial L^F}{\partial {}_{a}D_t^{\alpha }{q}_1}}=m{}_{a}D_t^{\alpha }{q}_2-{\displaystyle \frac{\rho }{2}}{q}_2,\hfill \\ p_2={\displaystyle \frac{\partial L^F}{\partial {}_{a}D_t^{\alpha }{q}_2}}=m{}_{a}D_t^{\alpha }{q}_1+{\displaystyle \frac{\rho }{2}}{q}_1.\hfill \end{array}$$

(14)

Applying the Legendre transformation, we obtain the Hamiltonian of fractional order

$$H^F(t,q_i,p_i)=\sum _i{p}_i{}_{a}D_t^{\alpha }{q}_i(q_i,p_i)-L(t,q_i,{}_{a}D_t^{\alpha }{q}_i(q_i,p_i)).$$

(15)

Using the Equation (15), we have

$$H^F=(K-{\displaystyle \frac{{\rho }^2}{4m}})q_1{q}_2+{\displaystyle \frac{\rho }{2m}}(q_2{p}_2-q_1{p}_1)+{\displaystyle \frac{p_1{p}_2}{m}}.$$

(16)

We define $\omega =\sqrt{K-\frac{{\rho }^2}{4m}}$ and the Hamiltonian takes the form

$$H^F={\omega }^2{q}_1{q}_2+{\displaystyle \frac{\rho }{2m}}(q_2{p}_2-q_1{p}_1)+{\displaystyle \frac{p_1{p}_2}{m}}.$$

(17)

The fractional Hamilton model of the BFT oscillator is given by

$$\begin{array}{c}{}_{a}D_t^{\alpha }{q}_1=-{\displaystyle \frac{\rho q_1}{2m}}+{\displaystyle \frac{p_2}{m}},\hfill \\ {}_{a}D_t^{\alpha }{q}_2={\displaystyle \frac{\rho q_2}{2m}}+{\displaystyle \frac{p_1}{m}},\hfill \\ {}_{a}D_t^{\alpha }{p}_1={\displaystyle \frac{{\rho }^2{q}_2}{4m}}+{\displaystyle \frac{\rho p_1}{2m}}-K{q}_2,\hfill \\ {}_{a}D_t^{\alpha }{p}_2={\displaystyle \frac{{\rho }^2{q}_1}{4m}}-{\displaystyle \frac{\rho p_2}{2m}}-K{q}_1.\hfill \end{array}$$

(18)

Now, we consider the fractional operators of Liouville–Caputo, Caputo–Fabrizio–Caputo and the fractional derivative based in the Mittag–Leffler kernel.

• First case. In the Liouville–Caputo sense, we have

$$\begin{array}{c}{q}_1\left(t\right)={\displaystyle \sum _{i=0}^{n-1}}{q}_1{\left(0\right)}^{\left(i\right)}{\displaystyle \frac{t^i}{i!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\displaystyle \underset{0}{\overset{t}{\int }}}{(t-\kappa )}^{\alpha -1}\left(-{\displaystyle \frac{\chi q_1\left(\kappa \right)}{2m}}+{\displaystyle \frac{p_2\left(\kappa \right)}{m}}\right)d\kappa ,\hfill \\ q_2\left(t\right)={\displaystyle \sum _{i=0}^{n-1}}{q}_2{\left(0\right)}^{\left(i\right)}{\displaystyle \frac{t^i}{i!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\displaystyle \underset{0}{\overset{t}{\int }}}{(t-\kappa )}^{\alpha -1}\left({\displaystyle \frac{\chi q_2\left(\kappa \right)}{2m}}+{\displaystyle \frac{p_1\left(\kappa \right)}{m}}\right)d\kappa ,t<T,\hfill \\ p_1\left(t\right)={\displaystyle \sum _{i=0}^{n-1}}{p}_1{\left(0\right)}^{\left(i\right)}{\displaystyle \frac{t^i}{i!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\displaystyle \underset{0}{\overset{t}{\int }}}{(t-\kappa )}^{\alpha -1}\left({\displaystyle \frac{{\chi }^2{q}_2\left(\kappa \right)}{4m}}+{\displaystyle \frac{\chi p_1\left(\kappa \right)}{2m}}-K{q}_2\left(\kappa \right)\right)d\kappa ,\hfill \\ p_2\left(t\right)={\displaystyle \sum _{i=0}^{n-1}}{p}_2{\left(0\right)}^{\left(i\right)}{\displaystyle \frac{t^i}{i!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\displaystyle \underset{0}{\overset{t}{\int }}}{(t-\kappa )}^{\alpha -1}\left({\displaystyle \frac{{\chi }^2{q}_1\left(\kappa \right)}{4m}}-{\displaystyle \frac{\chi p_2\left(\kappa \right)}{2m}}-K{q}_1\left(\kappa \right)\right)d\kappa .\hfill \end{array}$$

(19)

The numerical approximation of (19) is obtained using the algorithm (4).

• Second case. In the Caputo–Fabrizio–Caputo sense,

$$q_{1(l+1)}\left(t\right)=q_{\left(1\right)}\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[-\left(\frac{\rho }{2m}\right)q_{1(l+1)}\left(t\right)+\left(\frac{1}{m}\right)p_{2(l+1)}\left(t\right)\right]\right\}+$$

$$+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{\epsilon }_{1,z,l}·\left[-\left(\frac{\rho }{2m}\right)q_{1\left(l\right)}\left(t\right)+\left(\frac{1}{m}\right)p_{2\left(l\right)}\left(t\right)\right],$$

$$q_{2(l+1)}\left(t\right)=q_{\left(2\right)}\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(\frac{\rho }{2m}\right)q_{2(l+1)}\left(t\right)+\left(\frac{1}{m}\right)p_{1(l+1)}\left(t\right)\right]\right\}+$$

$$+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{\epsilon }_{2,z,l}·\left[\left(\frac{\rho }{2m}\right)q_{2\left(l\right)}\left(t\right)+\left(\frac{1}{m}\right)p_{1\left(l\right)}\left(t\right)\right],$$

$$p_{1(l+1)}\left(t\right)=p_{\left(1\right)}\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(\frac{{\rho }^2}{4m}\right)q_{2(l+1)}\left(t\right)+\left(\frac{\rho }{2m}\right)p_{1(l+1)}\left(t\right)-Z{q}_{2(l+1)}\left(t\right)\right]\right\}+$$

$$+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{\epsilon }_{3,z,l}·\left[\left(\frac{{\rho }^2}{4m}\right)q_{2\left(l\right)}\left(t\right)+\left(\frac{\rho }{2m}\right)p_{1\left(l\right)}\left(t\right)-Z{q}_{2\left(l\right)}\left(t\right)\right],$$

$$p_{2(l+1)}\left(t\right)=p_{\left(2\right)}\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(\frac{{\rho }^2}{4m}\right)q_{1(l+1)}\left(t\right)-\left(\frac{\rho }{2m}\right)p_{2(l+1)}\left(t\right)-Z{q}_{1(l+1)}\left(t\right)\right]\right\}+$$

$$+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{\epsilon }_{4,z,l}·\left[\left(\frac{{\rho }^2}{4m}\right)q_{1\left(l\right)}\left(t\right)-\left(\frac{\rho }{2m}\right)p_{2\left(l\right)}\left(t\right)-Z{q}_{1\left(l\right)}\left(t\right)\right],$$

(20)

where

$${\epsilon }_{(1,2,3,4),z,l+1}\left\{\begin{array}{c}{l}^{\alpha }-(l-\alpha ){(l+1)}^{\alpha },z=0,\hfill \\ {(l-z+2)}^{\alpha +1}+{(l-z)}^{\alpha +1}-2{(l-z+1)}^{\alpha +1},0\le z\le l.\hfill \end{array}\right.$$

• Third case. For the fractional derivative based on the Mittag–Leffler kernel, we used the numerical approximation scheme developed in [20]

$${}_0^{AB}{I}_t^{\alpha }\left[f\left(t_{l+1}\right)\right]=\frac{1-\alpha }{B\left(\alpha \right)}\left[\frac{f\left(t_{l+1}\right)-f\left(t_l\right)}{2}\right]+\frac{\alpha }{\mathsf{\Gamma }\left(\alpha \right)}\sum _{z=0}^{\infty }\left[\frac{f\left(t_{z+1}\right)-f\left(t_z\right)}{2}\right]b_z^{\alpha },$$

(21)

where

$$b_z^{\alpha }={(z+1)}^{1-\alpha }-{\left(z\right)}^{1-\alpha },$$

(22)

and the system (18) is represented by

$$q_{1(l+1)}\left(t\right)-q_{1\left(l\right)}\left(t\right)=q_{\left(1\right)}^l\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\right[-\left(\frac{\chi }{2m}\right)\left(\frac{q_{1(l+1)}\left(t\right)-q_{1\left(l\right)}\left(t\right)}{2}\right)+$$

$$+\left(\frac{1}{m}\right)\left(\frac{p_{2(l+1)}\left(t\right)-p_{2\left(l\right)}\left(t\right)}{2}\right)\left]\right\}+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{b}_z^{\alpha }·[-\left(\frac{\rho }{2m}\right)\left(\frac{q_{1(z+1)}\left(t\right)-q_{1\left(z\right)}\left(t\right)}{2}\right)+$$

$$+\left(\frac{1}{m}\right)\left(\frac{p_{2(z+1)}\left(t\right)-p_{2\left(z\right)}\left(t\right)}{2}\right)],$$

$$q_{2(l+1)}\left(t\right)-q_{2\left(l\right)}\left(t\right)=q_{\left(2\right)}^l\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\right[\left(\frac{\rho }{2m}\right)\left(\frac{q_{2(l+1)}\left(t\right)-q_{2\left(l\right)}\left(t\right)}{2}\right)+$$

$$+\left(\frac{1}{m}\right)\left(\frac{p_{1(l+1)}\left(t\right)-p_{1\left(l\right)}\left(t\right)}{2}\right)\left]\right\}+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{b}_z^{\alpha }·[\left(\frac{\rho }{2m}\right)\left(\frac{q_{2(z+1)}\left(t\right)-q_{2\left(z\right)}\left(t\right)}{2}\right)+$$

$$+\left(\frac{1}{m}\right)\left(\frac{p_{1(z+1)}\left(t\right)-p_{1\left(z\right)}\left(t\right)}{2}\right)],$$

$$p_{1(l+1)}\left(t\right)-p_{1\left(l\right)}\left(t\right)=p_{\left(1\right)}^l\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\right[\left(\frac{\rho }{2m}\right)\left(\frac{q_{2(l+1)}\left(t\right)-q_{2\left(l\right)}\left(t\right)}{2}\right)+$$

$$+\left(\frac{1}{m}\right)\left(\frac{p_{1(l+1)}\left(t\right)-p_{1\left(l\right)}\left(t\right)}{2}\right)-Z\left(\frac{q_{2(l+1)}\left(t\right)-q_{2\left(l\right)}\left(t\right)}{2}\right)\left]\right\}+$$

$$+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{b}_z^{\alpha }·[\left(\frac{\rho }{2m}\right)\left(\frac{q_{2(z+1)}\left(t\right)-q_{2\left(z\right)}\left(t\right)}{2}\right)+$$

$$+\left(\frac{1}{m}\right)\left(\frac{p_{1(z+1)}\left(t\right)-p_{1\left(z\right)}\left(t\right)}{2}\right)-Z\left(\frac{q_{2(z+1)}\left(t\right)-q_{2\left(z\right)}\left(t\right)}{2}\right)],$$

$$p_{2(l+1)}\left(t\right)-p_{2\left(l\right)}\left(t\right)=p_{\left(2\right)}^l\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\right[\left(\frac{\rho }{2m}\right)\left(\frac{q_{1(l+1)}\left(t\right)-q_{1\left(l\right)}\left(t\right)}{2}\right)-$$

$$-\left(\frac{\rho }{m}\right)\left(\frac{p_{2(l+1)}\left(t\right)-p_{2\left(l\right)}\left(t\right)}{2}\right)-Z\left(\frac{q_{1(l+1)}\left(t\right)-q_{1\left(l\right)}\left(t\right)}{2}\right)\left]\right\}+$$

$$+\frac{\alpha }{B\left(\alpha \right)}\sum _{z=0}^{\infty }{b}_z^{\alpha }·[\left(\frac{\rho }{2m}\right)\left(\frac{q_{1(z+1)}\left(t\right)-q_{1\left(z\right)}\left(t\right)}{2}\right)-$$

$$-\left(\frac{\rho }{m}\right)\left(\frac{p_{2(z+1)}\left(t\right)-p_{2\left(z\right)}\left(t\right)}{2}\right)-Z\left(\frac{q_{1(z+1)}\left(t\right)-q_{1\left(z\right)}\left(t\right)}{2}\right)].$$

(23)

Numerical Simulations

Figure 1, Figure 2 and Figure 3 shows the position $q_1=x_2\left(t\right)$, $q_2=x_1\left(t\right)$, ${}_{a}D_t^{\alpha }{x}_1\left(t\right)=x_3\left(t\right)$ and ${}_{a}D_t^{\alpha }{x}_2\left(t\right)=x_4\left(t\right)$ for systems (19), (20) and (23), respectively. For the simulation, the following values were considered: $m=5$, $\rho =2$, $K=0.1$ and different values of α, the total simulation time considered is 5 s, and the computational step $1\times {10}^{-3}$. The initial conditions $x_1\left(0\right)=1$, $x_2\left(0\right)=0.1$, $x_3\left(0\right)=1$ and $x_4\left(0\right)=0.5$ were considered. The results show that by keeping the parameters constant and by varying α, we obtain different results. The reported results illustrate that the fractional approach is more suitable to describe the complex dynamics of the investigated model.

3.2. Caldirola–Kanai Oscillator

We consider a harmonic CK oscillator whose mass depends on time such that $m\left(t\right)=mexp(sin\beta \gamma t)$, in this case, the Lagrangian is given by

$$L=exp(sin\beta \gamma t)[{\displaystyle \frac{1}{2}}m{\dot{q}}^2-{\displaystyle \frac{1}{2}}m{\omega }^2\left(t\right)q^2],$$

(24)

where m depends explicitly on time, and β and γ are variable parameter and damping factors.

The fractional Lagrangian (24) is given by

$$L^F=E_{\alpha ,1}(sin\beta \gamma t)[{\displaystyle \frac{1}{2}}m{(}_a{D}_t^{\alpha }{q}^2)-{\displaystyle \frac{1}{2}}m{\omega }^2\left(t\right)q^2],$$

(25)

and

$${}_{a}D_t^{\alpha }\left(E_{\alpha ,1}(sin\beta \gamma t){}_{a}D_t^{\alpha }q\right)-E_{\alpha ,1}(sin\beta \gamma t){\omega }^2\left(t\right)q=0.$$

(26)

The generalized momentum is

$$p_i={\displaystyle \frac{\partial L^F}{\partial {}_{a}D_t^{\alpha }{q}_i}},$$

(27)

$$\begin{array}{c}p={\displaystyle \frac{\partial L^F}{\partial {}_{a}D_t^{\alpha }q}}=E_{\alpha ,1}(sin\beta \gamma t)\left[m\left({}_{a}D_t^{\alpha }q\right)\right],\hfill \end{array}$$

(28)

where $L^F$ is the Lagrangian of fractional order of (24) with $i=1$, $q_1=q$ and $p_1=p$.

The Hamiltonian of fractional order is obtained using the Legendre transformation

$$H^F(t,q_i,p_i)=\sum _i{p}_i{}_{a}D_t^{\alpha }{q}_i(q_i,p_i)-\mathcal{L}(t,q_i{,}_a{D}_t^{\alpha }{q}_i(q_i,p_i)),$$

(29)

where

$${\mathcal{H}}^F={\displaystyle \frac{p^2}{2m}}{E}_{\alpha ,1}(-sin\beta \gamma t)+{\displaystyle \frac{m}{2}}{\omega }^2\left(t\right)q^2{E}_{\alpha ,1}(sin\beta \gamma t).$$

(30)

The fractional Hamilton model of the CK oscillator is given by

$$\begin{array}{c}{}_{a}D_t^{\alpha }q={\displaystyle \frac{p}{m}}{E}_{\alpha ,1}(-sin\beta \gamma t),\hfill \\ {}_{a}D_t^{\alpha }p=mq{\omega }^2\left(t\right)E_{\alpha ,1}(sin\beta \gamma t).\hfill \end{array}$$

(31)

Now, we consider the fractional operators of Liouville–Caputo, Caputo–Fabrizio–Caputo and the fractional derivative based on the Mittag–Leffler kernel.

• First case. In the Liouville–Caputo sense, we have

$$\begin{array}{c}q\left(t\right)={\displaystyle \sum _{i=0}^{n-1}}q{\left(0\right)}^{\left(i\right)}{\displaystyle \frac{t^i}{i!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}\underset{0}{\overset{t}{\int }}{(t-\kappa )}^{\alpha -1}({\displaystyle \frac{p\left(\kappa \right)}{m}}{E}_{\alpha ,1}(-sin\beta \gamma \kappa ))d\kappa ,\hfill \\ p\left(t\right)={\displaystyle \sum _{i=0}^{n-1}}p{\left(0\right)}^{\left(i\right)}{\displaystyle \frac{t^i}{i!}}+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}\underset{0}{\overset{t}{\int }}{(t-\kappa )}^{\alpha -1}\left(mq\left(\kappa \right){\omega }^2\left(\kappa \right)E_{\alpha ,1}(sin\beta \gamma \kappa )\right)d\kappa ,t<T.\hfill \end{array}$$

(32)

The numerical approximation of (32) is obtained using algorithm (4).

• Second case. In the Caputo–Fabrizio–Caputo sense, the Adams–Moulton rule for system (31) is given by

$$\begin{array}{c}{q}_{1(l+1)}\left(t\right)=q_{\left(1\right)}\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(\frac{1}{m}\right)p_{1(l+1)}\left(t\right)E_{\alpha ,1}(-sin\beta \gamma t)\right]\right\}+\\ +\frac{\alpha }{B\left(\alpha \right)}{\displaystyle \sum _{z=0}^{\infty }}{\epsilon }_{1,z,l}·\left[\left(\frac{1}{m}\right)p_{1\left(l\right)}\left(t\right)E_{\alpha ,1}(-sin\beta \gamma t)\right],\\ p_{1(l+1)}\left(t\right)=p_{\left(1\right)}\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(m{\omega }^2\left(t\right)\right)q_{1(l+1)}\left(t\right)E_{\alpha ,1}(sin\beta \gamma t)\right]\right\}+\\ +\frac{\alpha }{B\left(\alpha \right)}{\displaystyle \sum _{z=0}^{\infty }}{\epsilon }_{2,z,l}·\left[\left(m{\omega }^2\left(t\right)\right)q_{1\left(l\right)}\left(t\right)E_{\alpha ,1}(sin\beta \gamma t)\right],\end{array}$$

(33)

where

$${\epsilon }_{(1,2),z,l+1}\left\{\begin{array}{c}{l}^{\alpha }-(l-\alpha ){(l+1)}^{\alpha },z=0,\hfill \\ {(l-z+2)}^{\alpha +1}+{(l-z)}^{\alpha +1}-2{(l-z+1)}^{\alpha +1},0\le z\le l.\hfill \end{array}\right.$$

• Third case. For the fractional derivative based on the Mittag–Leffler kernel, we have

$$\begin{array}{c}{q}_{1(l+1)}\left(t\right)-q_{1\left(l\right)}\left(t\right)=q_{\left(1\right)}^l\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(\frac{1}{m}\right)\left(\frac{p_{1(l+1)}\left(t\right)-p_{1\left(l\right)}\left(t\right)}{2}\right)·(-sin\beta \gamma t)\right]\right\}+\\ +\frac{\alpha }{B\left(\alpha \right)}{\displaystyle \sum _{z=0}^{\infty }}{b}_z^{\alpha }·\left[\left(\frac{1}{m}\right)\left(\frac{p_{1(z+1)}\left(t\right)-p_{1\left(z\right)}\left(t\right)}{2}\right)·E_{\alpha ,1}(-sin\beta \gamma t)\right],\\ p_{1(l+1)}\left(t\right)-p_{1\left(l\right)}\left(t\right)=p_{\left(1\right)}^l\left(t\right)+\left\{\frac{1-\alpha }{B\left(\alpha \right)}\left[\left(m{\omega }^2\left(t\right)\right)\left(\frac{p_{1(l+1)}\left(t\right)-p_{1\left(l\right)}\left(t\right)}{2}\right)·(sin\beta \gamma t)\right]\right\}+\\ +\frac{\alpha }{B\left(\alpha \right)}{\displaystyle \sum _{z=0}^{\infty }}{b}_z^{\alpha }·\left[\left(m{\omega }^2\left(t\right)\right)\left(\frac{p_{1(z+1)}\left(t\right)-p_{1\left(z\right)}\left(t\right)}{2}\right)·E_{\alpha ,1}(sin\beta \gamma t)\right].\end{array}$$

(34)

Numerical Simulations

Figure 4, Figure 5 and Figure 6 depicted the numerical evaluation of (32)–(34) in Liouville–Caputo, Caputo–Fabrizio–Caputo and the fractional derivative based on the Mittag–Leffler kernel, respectively, considering different values of $\omega \left(t\right)$ and fractional order γ, for all cases $a=0$ and $b=1$, and the total simulation time considered is one second and computational step $1\times {10}^{-5}$. It is clear from the figures that the behaviors of the fractional equations strongly depend on the order α of the fractional derivatives, in addition to the form of the function $w\left(t\right)$.

4. Conclusions

Alternative representations of the Bateman–Feshbach–Tikochinsky and Caldirola–Kanai oscillators were studied using fractional operators of Liouville–Caputo type. We derive new solutions of these models using an iterative scheme and via a Crank–Nicholson scheme. The Liouville–Caputo fractional derivative involves a kernel with singularity, and this definition is based on the power law and present singularity at the origin. Recently, Caputo and Fabrizio solved the problem of singularity at the origin and used the exponential decay law to construct a derivative with no singularity; however, the used kernel is local. This derivative therefore has an advantage over the Liouville–Caputo derivative because the full effect of the memory can be portrayed. Atangana and Baleanu suggested two fractional derivatives based on the generalized Mittag–Leffler function. These derivatives with fractional order in Liouville–Caputo and Riemann–Liouville sense have non-singular and non-local kernel and preserve the benefits of the Riemann–Liouville, Liouville–Caputo and Caputo–Fabrizio operators.

Using these fractional operators, the results show that, by keeping the parameters constant and by varying α, we obtain different behaviors. The reported results illustrate that the fractional approach is more suitable to describe the complex dynamics of the investigated models. Finally, we observe novel behaviors that cannot be obtained with standard models and using local derivatives.