1. Introduction
Many classical and quantum systems with position-dependent mass have been already investigated [5–7]. The Mathews–Lakshmanan oscillator is a prominent example in which a particle with position-dependent mass defined by $m=1/(1+\lambda {x}^{2})$ moves in a harmonic potential [8]. For more details about these models, refer to [9–11]. In addition, for the point moving along the well-known Lemniscate curve in classical mechanics [12, 13], the effective mass can be obtained from $m=1/(1+{x}^{4})$.
Fractional calculus dates from more than 300 years ago [14]. It deals with derivatives and integrals to any order, not only integer. In the last 40 years, fractional calculus has found a wide range of applications in many branches of science and engineering [15–21]. New aspects of complicated dynamics with memory trace in many physical systems can be taken into account by fractional calculus. Nevertheless, the singularity of the classical fractional operators means that nonlocal dynamical behavior of real-world physical systems cannot be displayed exactly. To overcome this difficulty, various kinds of fractional derivatives and integrals have been introduced whose kernel is nonsingular and, hence, they can specify nonlocal dynamics accurately. The fractional differential operator with the Mittag–Leffler (ML) kernel, introduced by Atangana and Baleanu (ABC) [22], is one of the most applicable. The nonlocality of real-world complicated phenomena can be more accurately modeled by the ABC derivative compared with classical fractional derivatives [23–26]. Moreover, since the dynamical behavior of practical physical systems is studied via the Lagrangian or Hamiltonian in a fractional framework, it is important to find a high-performance numerical method for solving the related fractional Euler–Lagrange equations (FELEs) in the ABC sense. Note that, because of the effect of memory, the numerical solving methods in fractional calculus are not a straightforward extension of the related classical methods. Thus, designing numerical schemes for solving real-world dynamics with related fractional models should be considered. In this paper, we investigate the free motion of a harmonic oscillator with position-dependent mass by applying its new fractional formulation. The derived FELEs in the ABC concept and the numerical procedure for solving them that are suggested for the harmonic oscillator are new and include more accurate information compared with their equivalent fractional equations in the standard form.
This paper is organized as follows. The basic definitions and preliminaries are discussed in section 2 . The descriptions of the 1D and 2D systems are illustrated and discussed in section 3 . In section 4 , the system is solved numerically using the fractional derivative operator ABC. The results and discussions are presented in section 5 , and we conclude the paper in section 6 .
2. Basic definitions and preliminaries
In this section, we represent the fractional derivatives in the concept of Caputo [14] and ABC [22]. For a function $f:[a,b]\to {{\mathbb{R}}}^{n},0\lt \alpha \lt 1$, the left and right fractional derivatives of order α in the concept of Caputo are respectively defined by [14]
$\begin{eqnarray}{}_{\hspace{0.02cm}a}^{C}\hspace{-0.05cm}{{\rm{D}}}_{t}^{\alpha }f(t)=\displaystyle \frac{1}{{\rm{\Gamma }}(1-\alpha )}{\int }_{a}^{t}{\left(t-\gamma \right)}^{-\alpha }\dot{f}(\gamma ){\rm{d}}\gamma ,\end{eqnarray}$
$\begin{eqnarray}{}_{\hspace{0.02cm}t}^{C}\hspace{-0.05cm}{{\rm{D}}}_{b}^{\alpha }f(t)=-\displaystyle \frac{1}{{\rm{\Gamma }}(1-\alpha )}{\int }_{t}^{b}{\left(\gamma -t\right)}^{-\alpha }\dot{f}(\gamma ){\rm{d}}\gamma ,\end{eqnarray}$
where Γ(.) is the Euler’s Gamma function. Moreover, for $f\in {H}^{1}(a,b)$ and 0 < α < 1, the left and right fractional derivatives of order α in the concept of ABC are respectively determined as [22] $\begin{eqnarray}{}_{a}^{{\rm{ABC}}}{{\rm{D}}}_{t}^{\alpha }f(t)=\displaystyle \frac{M(\alpha )}{1-\alpha }{\int }_{a}^{t}{{\rm{E}}}_{\alpha }[-\beta {\left(t-\gamma \right)}^{\alpha }]\dot{f}(\gamma ){\rm{d}}\gamma ,\end{eqnarray}$
$\begin{eqnarray}{}_{t}^{{\rm{ABC}}}{{\rm{D}}}_{b}^{\alpha }f(t)=\displaystyle \frac{M(\alpha )}{1-\alpha }{\int }_{t}^{b}{{\rm{E}}}_{\alpha }[-\beta {\left(\gamma -t\right)}^{\alpha }]\dot{f}(\gamma ){\rm{d}}\gamma ,\end{eqnarray}$
where $\beta =\tfrac{\alpha }{1-\alpha },{{\rm{E}}}_{\alpha }(t)={\sum }_{l=0}^{\infty }\tfrac{{t}^{l}}{{\rm{\Gamma }}(\alpha l+1)}$ is the ML function, and M(α) is a normalization function with M(0) = M(1) = 1.3. The physical system
3.1. 1D system
We start this section by introducing a definition to our physical system of interest. As a starting point, let us take a point mass that moves along the Lemniscate of Bernoulli, as shown in figure 1. In 1694, Jacob Bernoulli published an article about a curve shaped like figure 1, or a knot or bow of a ribbon. This curve is called the Lemniscate of Bernoulli, which has the following Cartesian equation9 ), we yield the following classical Euler–Lagrange equation (CELE)9 ) in the sense of the fractional calculus. As a result, equation (9 ) is fractionalized as follows11 ) into $\tfrac{\partial {L}_{f}}{\partial x}{+}_{t}\hspace{-0.05cm}{{\rm{D}}}_{b}^{\alpha }\tfrac{\partial {L}_{f}}{{\partial }_{a}\hspace{-0.05cm}{{\rm{D}}}_{t}^{\alpha }x}\,+{}_{a}{{\rm{D}}}_{t}^{\alpha }\tfrac{\partial {L}_{f}}{{\partial }_{t}\hspace{-0.05cm}{{\rm{D}}}_{b}^{\alpha }x}=0$, the FELE reads12 ) reduces to the CELE (10 ).
$\begin{eqnarray}{\left(X+Y\right)}^{2}=\displaystyle \frac{1}{2}({X}^{2}-{Y}^{2}).\end{eqnarray}$
Assume that the bead is subjected to the harmonic potential given by $\begin{eqnarray}U={\epsilon }^{2}({X}^{2}+{Y}^{2})/2.\end{eqnarray}$
This means that a stretched rubber band sticks the bead to the origin. The kinetic and the potential energies equations for the bead are respectively described by [27] $\begin{eqnarray}T=\displaystyle \frac{1}{2}({\dot{X}}^{2}+{\dot{Y}}^{2})=\displaystyle \frac{1}{2}m(x){\dot{x}}^{2},\end{eqnarray}$
$\begin{eqnarray}U=\displaystyle \frac{{\epsilon }^{2}}{2}m(x){x}^{2}.\end{eqnarray}$
A main feature of the harmonic potential considered here is that it preserves its form after pullback. Therefore, the classical Lagragian of the system reads $\begin{eqnarray}{L}_{c}=\displaystyle \frac{1}{2}m(x){\dot{x}}^{2}-\displaystyle \frac{{\epsilon }^{2}}{2}m(x){x}^{2},\end{eqnarray}$
where $m(x)=\tfrac{1}{1+{x}^{4}}$ [27]. Applying $\tfrac{\partial {L}_{c}}{\partial x}-\tfrac{{\rm{d}}}{{\rm{d}}t}\tfrac{\partial {L}_{c}}{\partial \dot{x}}=0$ to equation ( $\begin{eqnarray}\ddot{x}-\displaystyle \frac{2{x}^{3}}{1+{x}^{4}}{\dot{x}}^{2}+{\epsilon }^{2}\displaystyle \frac{1-{x}^{4}}{1+{x}^{4}}x=0.\end{eqnarray}$
Note that the last equation is classified as the quadratic Lienard-type equation. However, as was mentioned in [28], the theory of the calculus of variations cannot capture many laws of the natural phenomena; for instance, the behavior of nonconservative systems cannot be described properly by dynamical equations derived based on the traditional energy approach. On the contrary, the fractional calculus has overcome this limitation as it can characterize the behavior of many complex physical systems, including hereditary effects. Therefore, following the procedures explained in [29, 30], we can generalize the classical Lagrangain ( $\begin{eqnarray}{L}_{f}=\displaystyle \frac{1}{2}{\left({}_{a}{{\rm{D}}}_{t}^{\alpha }x\right)}^{2}m(x)-\displaystyle \frac{{\epsilon }^{2}}{2}m(x){x}^{2}.\end{eqnarray}$
By substituting equation ( $\begin{eqnarray}{}_{t}{{\rm{D}}}_{b}^{\alpha }{\displaystyle \frac{1}{1+{x}^{4}}}_{a}{{\rm{D}}}_{t}^{\alpha }x-\displaystyle \frac{2{x}^{3}}{{\left(1+{x}^{4}\right)}^{2}}{\left({}_{a}{{\rm{D}}}_{t}^{\alpha }x\right)}^{2}-{\epsilon }^{2}\displaystyle \frac{1-{x}^{4}}{{\left(1+{x}^{4}\right)}^{2}}x=0.\end{eqnarray}$
As $\alpha \to 1$, the FELE (
Figure 1. A bead on a Lemniscate. |
3.2. 2D system
Following the derivation given in [27], we can write the classical Lagrangian of the 2D systems as13 ) for q1 = r and q2 = φ; thus, we obtain, respectively13 ) can be generalized as17 )–(18 ) reduce to the CELEs (14 )–(15 ).
$\begin{eqnarray}{L}_{c}=\displaystyle \frac{1}{2}m(r)[{\dot{r}}^{2}+{r}^{2}{\dot{\phi }}^{2}]-\displaystyle \frac{{\epsilon }^{2}}{2}m(r){r}^{2},\end{eqnarray}$
where the mass is again position-dependent as $m(r)=\tfrac{1}{1+{r}^{4}}$. Next, we apply $\tfrac{\partial {L}_{c}}{\partial {q}_{i}}-\tfrac{{\rm{d}}}{{\rm{d}}t}\tfrac{\partial {L}_{c}}{\partial \dot{{q}_{i}}}=0$ to equation ( $\begin{eqnarray}\ddot{r}-\displaystyle \frac{2{r}^{3}}{1+{r}^{4}}{\dot{r}}^{2}+{\epsilon }^{2}\displaystyle \frac{1-{r}^{4}}{1+{r}^{4}}r-r\displaystyle \frac{1-{r}^{4}}{1+{r}^{4}}{\dot{\phi }}^{2}=0,\end{eqnarray}$
$\begin{eqnarray}\ddot{\phi }+\displaystyle \frac{2}{r}\displaystyle \frac{1-{r}^{4}}{1+{r}^{4}}\dot{r}\dot{\phi }=0.\end{eqnarray}$
Again, equation ( $\begin{eqnarray}{L}_{f}=\displaystyle \frac{1}{2}m(r)[{\left({}_{a}{{\rm{D}}}_{t}^{\alpha }r\right)}^{2}+{r}^{2}{\left({}_{a}{{\rm{D}}}_{t}^{\alpha }\phi \right)}^{2}]-\displaystyle \frac{{\epsilon }^{2}}{2}m(r){r}^{2}.\end{eqnarray}$
As was the case in the 1D system, the FELEs read $\begin{eqnarray}\begin{array}{l}{}_{t}{{\rm{D}}}_{b}^{\alpha }\displaystyle \frac{1}{1+{r}^{4}}{}_{a}{{\rm{D}}}_{t}^{\alpha }r+\displaystyle \frac{r-{r}^{5}}{{\left(1+{r}^{4}\right)}^{2}}{\left({}_{a}{{\rm{D}}}_{t}^{\alpha }\phi \right)}^{2}\\ \quad -\,\displaystyle \frac{2{r}^{3}}{{\left(1+{r}^{4}\right)}^{2}}{\left({}_{a}{{\rm{D}}}_{t}^{\alpha }r\right)}^{2}-{\epsilon }^{2}\displaystyle \frac{r-{r}^{5}}{{\left(1+{r}^{4}\right)}^{2}}=0,\end{array}\end{eqnarray}$
$\begin{eqnarray}{}_{t}{{\rm{D}}}_{b}^{\alpha }\displaystyle \frac{{r}^{2}}{1+{r}^{4}}{}_{a}{{\rm{D}}}_{t}^{\alpha }\phi =0.\end{eqnarray}$
Again, as $\alpha \to 1$, the FELEs (4. Numerical procedure
In this section, we propose an impressive numerical method which solves the FELEs in equation (12 ) and equations (17 )–(18 ) for the 1D and 2D systems, respectively. The fractional derivative operator ABC with several values of α is taken into account. Also, we rewrite equation (12 ), and similarly equations (17 )–(18 ), by considering the new variables ${x}_{1}\,=\,x,{x}_{2}{=}_{a}{{\rm{D}}}_{t}^{\alpha }x,{r}_{1}=r,{r}_{2}{=}_{a}{{\rm{D}}}_{t}^{\alpha }r,{\phi }_{1}=\phi $ and ${\phi }_{2}{=}_{a}{{\rm{D}}}_{t}^{\alpha }\phi $ as follows:
| • | for the 1D system: $\begin{eqnarray}\begin{array}{c}\left\{\begin{array}{l}{}_{a}^{{\rm{ABC}}}{{\rm{D}}}_{t}^{\alpha }{x}_{1}={x}_{2},\\ {}_{a}^{{\rm{ABC}}}{{\rm{D}}}_{t}^{\alpha }\tfrac{1}{1+{x}_{1}^{4}}{x}_{2}=\tfrac{2{x}_{1}^{3}}{{\left(1+{x}_{1}^{4}\right)}^{2}}{x}_{2}^{2}+{\epsilon }^{2}\tfrac{1-{x}_{1}^{4}}{{\left(1+{x}_{1}^{4}\right)}^{2}}{x}_{1},\end{array}\right.\end{array}\end{eqnarray}$ |
| • | for the 2D system: $\begin{eqnarray}\begin{array}{c}\left\{\begin{array}{l}{}_{a}^{{\rm{ABC}}}{{\rm{D}}}_{t}^{\alpha }{r}_{1}={r}_{2},\\ {}_{t}^{{\rm{ABC}}}{{\rm{D}}}_{b}^{\alpha }\tfrac{1}{1+{r}_{1}^{4}}{r}_{2}=-\tfrac{{r}_{1}-{r}_{1}^{5}}{{\left(1+{r}_{1}^{4}\right)}^{2}}{\phi }_{2}^{2}+\tfrac{2{r}_{1}^{3}}{{\left(1+{r}_{1}^{4}\right)}^{2}}{r}_{2}^{2}+{\epsilon }^{2}\tfrac{{r}_{1}-{r}_{1}^{5}}{{\left(1+{r}_{1}^{4}\right)}^{2}},\\ {}_{a}^{{\rm{ABC}}}{{\rm{D}}}_{t}^{\alpha }{\phi }_{1}={\phi }_{2},\\ {}_{t}^{{\rm{ABC}}}{{\rm{D}}}_{b}^{\alpha }\tfrac{{r}_{1}^{2}}{1+{r}_{1}^{4}}{\phi }_{2}=0.\end{array}\right.\end{array}\,\end{eqnarray}$ Exerting the ABC integral operator [22] into equations ( $\begin{eqnarray}\left\{\begin{array}{l}{x}_{1}(t)={x}_{1}(a)+\tfrac{1-\alpha }{N(\alpha )}{x}_{2}(t)+\tfrac{\alpha }{{\rm{\Gamma }}(\alpha )N(\alpha )}{\displaystyle \int }_{a}^{t}{x}_{2}(\gamma )\\ \times \,{\left(t-\gamma \right)}^{\alpha -1}{\rm{d}}\gamma ,\\ \\ {x}_{2}(t)=\tfrac{1-\alpha }{N(\alpha )}\left(\tfrac{2{x}_{1}^{3}(t)}{1+{x}_{1}^{4}(t)}{x}_{2}^{2}(t)+{\epsilon }^{2}\tfrac{1-{x}_{1}^{4}(t)}{1+{x}_{1}^{4}(t)}{x}_{1}(t\right)\\ +\tfrac{\alpha (1+{x}_{1}^{4}(t))}{{\rm{\Gamma }}(\alpha )N(\alpha )}{\displaystyle \int }_{t}^{b}\left(\tfrac{2{x}_{1}^{3}(\gamma )}{{\left(1+{x}_{1}^{4}(\gamma \right)}^{2}}{x}_{2}^{2}(\gamma )+{\epsilon }^{2}\tfrac{1-{x}_{1}^{4}(\gamma )}{{\left(1+{x}_{1}^{4}(\gamma \right)}^{2}}{x}_{1}(\gamma \right)\\ \times \,{\left(\gamma -t\right)}^{\alpha -1}{\rm{d}}\gamma ,\end{array}\right.\end{eqnarray}$ $\begin{eqnarray}\left\{\begin{array}{l}{r}_{1}(t)={r}_{1}(a)+\tfrac{1-\alpha }{N(\alpha )}{r}_{2}(t)+\tfrac{\alpha }{{\rm{\Gamma }}(\alpha )N(\alpha )}{\displaystyle \int }_{a}^{t}{r}_{2}(\gamma )\\ {\left(t-\gamma \right)}^{\alpha -1}{\rm{d}}\gamma ,\\ \\ {r}_{2}(t)=\tfrac{1-\alpha }{N(\alpha )}\left(-\tfrac{{r}_{1}(t)-{r}_{1}^{5}(t)}{1+{r}_{1}^{4}(t)}{\phi }_{2}^{2}(t)+\tfrac{2{r}_{1}^{3}(t)}{1+{r}_{1}^{4}(t)}{r}_{2}^{2}(t)\right.\\ +\left.{\epsilon }^{2}\tfrac{{r}_{1}(t)-{r}_{1}^{5}(t)}{1+{r}_{1}^{4}(t)}\right)\\ +\tfrac{\alpha (1+{r}_{1}^{4}(t))}{{\rm{\Gamma }}(\alpha )N(\alpha )}{\displaystyle \int }_{t}^{b}\left(-\tfrac{{r}_{1}(\gamma )-{r}_{1}^{5}(\gamma )}{{\left(1+{r}_{1}^{4}(\gamma \right)}^{2}}{\phi }_{2}^{2}(\gamma )+\tfrac{2{r}_{1}^{3}(\gamma )}{{\left(1+{r}_{1}^{4}(\gamma \right)}^{2}}{r}_{2}^{2}(\gamma )\right.\\ +\left.{\epsilon }^{2}\tfrac{{r}_{1}(\gamma )-{r}_{1}^{5}(\gamma )}{{\left(1+{r}_{1}^{4}(\gamma \right)}^{2}}\right){\left(\gamma -t\right)}^{\alpha -1}{\rm{d}}\gamma ,\\ {\phi }_{1}(t)={\phi }_{1}(a)+\tfrac{1-\alpha }{N(\alpha )}{\phi }_{2}(t)+\tfrac{\alpha }{{\rm{\Gamma }}(\alpha )N(\alpha )}{\displaystyle \int }_{a}^{t}{\phi }_{2}(\gamma )\\ \times \,{\left(t-\gamma \right)}^{\alpha -1}{\rm{d}}\gamma ,\\ \\ {\phi }_{2}(t)=\tfrac{\alpha }{{\rm{\Gamma }}(\alpha )N(\alpha )}\left({r}_{1}^{2}(t)+\tfrac{1}{{r}_{1}^{2}(t)}\right)(b-t),\end{array}\right.\end{eqnarray}$ where x1(a), r1(a) and φ1(a) are the initial values and x2(b) = r2(b) = φ2(b) = 0. Now, we consider the length of time step ${{\ell }}_{M}=\tfrac{b-a}{M}$ for a uniform partition on [a, b] where M is a positive integer that can be selected optionally. Moreover, we represent the numerical approximations of xp(tq), rp(tq), and φp(tq) by ${x}_{p,q},{r}_{p,q}$, and ${\phi }_{p,q}$, respectively, where p = 1, 2, and ${t}_{q}\,=a+q{{\ell }}_{M},0\leqslant q\leqslant M$ is the time instant at the q-th node. Afterwards, the fractional Euler approach [31] is applied to discretize the convolution integrals in equations ( $\begin{eqnarray}\left\{\begin{array}{l}{X}_{1}-\tfrac{1-\alpha }{N(\alpha )}{X}_{2}-\tfrac{\alpha }{N(\alpha )}{A}_{M,\alpha }{X}_{2}={X}_{\mathrm{1,0}},\\ \\ {X}_{2}-\tfrac{1-\alpha }{N(\alpha )}{{\rm{\Lambda }}}_{1}({X}_{1},{X}_{2})-\tfrac{\alpha }{N(\alpha )}{{\rm{\Lambda }}}_{2}({X}_{1}){B}_{M,\alpha }{{\rm{\Lambda }}}_{3}({X}_{1},{X}_{2})=0,\end{array}\right.\end{eqnarray}$ $\begin{eqnarray}\left\{\begin{array}{l}{R}_{1}-\tfrac{1-\alpha }{N(\alpha )}{R}_{2}-\tfrac{\alpha }{N(\alpha )}{A}_{M,\alpha }{R}_{2}={R}_{\mathrm{1,0}},\\ \\ {R}_{2}-\tfrac{1-\alpha }{N(\alpha )}{{\rm{\Lambda }}}_{4}({R}_{1},{R}_{2},{{\rm{\Theta }}}_{2})-\tfrac{\alpha }{N(\alpha )}{{\rm{\Lambda }}}_{5}({R}_{1}){B}_{M,\alpha }{{\rm{\Lambda }}}_{6}\\ ({R}_{1},{R}_{2},{{\rm{\Theta }}}_{2})=0,\\ \\ {{\rm{\Theta }}}_{1}-\tfrac{1-\alpha }{N(\alpha )}{{\rm{\Theta }}}_{2}-\tfrac{\alpha }{N(\alpha )}{A}_{M,\alpha }{{\rm{\Theta }}}_{2}={{\rm{\Theta }}}_{\mathrm{1,0}},\\ \\ {{\rm{\Theta }}}_{2}-\tfrac{\alpha }{N(\alpha )}{{\rm{\Lambda }}}_{7}({R}_{1}){B}_{M,\alpha }=0,\end{array}\right.\end{eqnarray}$ where $\begin{eqnarray}\begin{array}{rcl}{A}_{M,\alpha } & = & {\left({B}_{M,\alpha }\right)}^{{\rm{T}}}={{\ell }}_{M}\left[\begin{array}{cccc}{\tau }_{0,\alpha } & 0 & \ldots & 0\\ {\tau }_{1,\alpha } & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0\\ {\tau }_{M,\alpha } & \ldots & {\tau }_{1,\alpha } & {\tau }_{0,\alpha }\end{array}\right],\\ {\tau }_{0,\alpha } & = & 1,\hspace{0.2cm}{\tau }_{q,\alpha }=\left(1+\displaystyle \frac{\alpha -1}{q}\right){\tau }_{q-1,\alpha },\hspace{0.2cm}q=1,2,\ldots ,\end{array}\end{eqnarray}$ $\begin{eqnarray}\begin{array}{rcl}{X}_{p} & = & \left[\begin{array}{c}{x}_{p,0}\\ \vdots \\ {x}_{p,M}\end{array}\right],\hspace{0.1cm}{R}_{p}=\left[\begin{array}{c}{R}_{p,0}\\ \vdots \\ {R}_{p,M}\end{array}\right],\hspace{0.1cm}{{\rm{\Theta }}}_{p}=\left[\begin{array}{c}{{\rm{\Theta }}}_{p,0}\\ \vdots \\ {{\rm{\Theta }}}_{p,M}\end{array}\right],\\ {X}_{p,0} & = & \left[\begin{array}{c}{x}_{p,0}\\ \vdots \\ {x}_{p,0}\end{array}\right],\hspace{0.1cm}{R}_{p,0}=\left[\begin{array}{c}{R}_{p,0}\\ \vdots \\ {R}_{p,0}\end{array}\right],\hspace{0.1cm}{{\rm{\Theta }}}_{p,0}=\left[\begin{array}{c}{{\rm{\Theta }}}_{p,0}\\ \vdots \\ {{\rm{\Theta }}}_{p,0}\end{array}\right],\hspace{0.1cm}p=1,2,\end{array}\end{eqnarray}$ $\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{{\rm{\Lambda }}}_{1}({X}_{1},{X}_{2}) & = & \left[\begin{array}{c}\tfrac{2{x}_{1,0}^{3}}{1+{x}_{1,0}^{4}}{x}_{2,0}^{2}+{\epsilon }^{2}\tfrac{1-{x}_{1,0}^{4}}{1+{x}_{1,0}^{4}}{x}_{1,0}\\ \vdots \\ \tfrac{2{x}_{1,M}^{3}}{1+{x}_{1,M}^{4}}{x}_{2,M}^{2}+{\epsilon }^{2}\tfrac{1-{x}_{1,M}^{4}}{1+{x}_{1,M}^{4}}{x}_{1,M}\end{array}\right],\\ {{\rm{\Lambda }}}_{3}({X}_{1},{X}_{2}) & = & \left[\begin{array}{c}\tfrac{2{x}_{1,0}^{3}}{{\left(1+{x}_{1,0}^{4}\right)}^{2}}{x}_{2,0}^{2}+{\epsilon }^{2}\tfrac{1-{x}_{1,0}^{4}}{{\left(1+{x}_{1,0}^{4}\right)}^{2}}{x}_{1,0}\\ \vdots \\ \tfrac{2{x}_{1,M}^{3}}{{\left(1+{x}_{1,M}^{4}\right)}^{2}}{x}_{2,M}^{2}+{\epsilon }^{2}\tfrac{1-{x}_{1,M}^{4}}{{\left(1+{x}_{1,M}^{4}\right)}^{2}}{x}_{1,M}\end{array}\right],\\ {{\rm{\Lambda }}}_{2}({X}_{1},{X}_{2}) & = & {\rm{diag}}\left[\begin{array}{ccc}1+{x}_{1,0}^{4} & \ldots & 1+{x}_{1,M}^{4}\end{array}\right],\end{array}\end{array}\end{eqnarray}$ as well as $\begin{eqnarray}\begin{array}{l}{{\rm{\Lambda }}}_{4}({R}_{1},{R}_{2},{{\rm{\Theta }}}_{2})\\ =\,\left[\begin{array}{c}-\tfrac{{r}_{\mathrm{1,0}}-{r}_{1,0}^{5}}{1+{r}_{1,0}^{4}}{\phi }_{2,0}^{2}+\tfrac{2{r}_{1,0}^{3}}{1+{r}_{1,0}^{4}}{r}_{2,0}^{2}+{\epsilon }^{2}\tfrac{{r}_{\mathrm{1,0}}-{r}_{1,0}^{5}}{1+{r}_{1,0}^{4}}\\ \vdots \\ -\tfrac{{r}_{1,M}-{r}_{1,M}^{5}}{1+{r}_{1,M}^{4}}{\phi }_{2,M}^{2}+\tfrac{2{r}_{1,M}^{3}}{1+{r}_{1,M}^{4}}{r}_{2,M}^{2}+{\epsilon }^{2}\tfrac{{r}_{1,M}-{r}_{1,M}^{5}}{1+{r}_{1,M}^{4}}\end{array}\right],\\ {{\rm{\Lambda }}}_{6}({R}_{1},{R}_{2},{{\rm{\Theta }}}_{2})\\ =\,\left[\begin{array}{c}-\tfrac{{r}_{\mathrm{1,0}}-{r}_{1,0}^{5}}{{\left(1+{r}_{\mathrm{1,0}}^{4}\right)}^{2}}{\phi }_{2,0}^{2}+\tfrac{2{r}_{1,0}^{3}}{{\left(1+{r}_{\mathrm{1,0}}^{4}\right)}^{2}}{r}_{2,0}^{2}+{\epsilon }^{2}\tfrac{{r}_{\mathrm{1,0}}-{r}_{1,0}^{5}}{{\left(1+{r}_{\mathrm{1,0}}^{4}\right)}^{2}}\\ \vdots \\ -\tfrac{{r}_{1,M}-{r}_{1,M}^{5}}{{\left(1+{r}_{1,M}^{4}\right)}^{2}}{\phi }_{2,M}^{2}+\tfrac{2{r}_{1,M}^{3}}{{\left(1+{r}_{1,M}^{4}\right)}^{2}}{r}_{2,M}^{2}+{\epsilon }^{2}\tfrac{{r}_{1,M}-{r}_{1,M}^{5}}{{\left(1+{r}_{1,M}^{4}\right)}^{2}}\end{array}\right],\end{array}\end{eqnarray}$ $\begin{eqnarray*}\begin{array}{c}\begin{array}{rcl}{{\rm{\Lambda }}}_{5}({R}_{1}) & = & {\rm{diag}}\left[\begin{array}{ccc}1+{r}_{1,0}^{4} & \ldots & 1+{r}_{1,M}^{4}\end{array}\right],\\ {{\rm{\Lambda }}}_{7}({R}_{1}) & = & {\rm{diag}}\left[\begin{array}{ccc}\left({r}_{1,0}^{2}+\tfrac{1}{{r}_{1,0}^{2}}\right)(b-0) & \ldots & \left({r}_{1,M}^{2}+\tfrac{1}{{r}_{1,M}^{2}}\right)(b-M)\end{array}\right].\end{array}\end{array}\end{eqnarray*}$ |
5. Simulation results and discussion
In this section, the behavior of the FELEs of motion for the harmonic oscillator with position-dependent mass for both the 1D and 2D cases are investigated by considering different values of the fractional order α. To this aim, we consider the following cases for the physical system under investigation.
Case one (1D system with x(0) = 0.1 and ϵ = 1): simulation results for x(t) against time t are shown in figure 2 for the fractional orders α = 0.85, 0.88, 0.91, 0.94, 0.97, and 1.0. It is clear from figure 2 that, as the fractional order α approaches 1, the fractional simulations become closer to the classical simulation, and they coincide with each other for the case α = 1. Also, figure 2 indicates that the simulation results for x(t) nearly look like a damping oscillator for the fractional values of α, while behaving like a harmonic oscillator when α = 1.0.
Figure 2. Simulation results of x(t) for the 1D case when x(0) = 0.1 and ϵ = 1. |
Case two (2D system with r(0) = 0.1, φ(0) = 0 and ϵ = 1): simulation results for r(t) against time t are shown in figure 3 for the fractional orders α = 0.85, 0.88, 0.91, 0.94, 0.97, and 1.0. Again, it is clear from figure 3 that, as the fractional order α approaches 1, the fractional simulations become closer to the classical simulation, and they coincide with each other for the case α = 1. Also, it is apparent from figure 3 that the simulation results for x(t) nearly look like a damping oscillator for the fractional values of α, whereas they behave like a harmonic oscillator for α = 1.0. Furthermore, simulation results for φ(t) versus time t are shown in figure 4 for the fractional orders α = 0.85, 0.88, 0.91, 0.94, 0.97, and 1.0. Figures 2–4 indicate that the FELEs numerical solution demonstrates various asymptotic behaviors considering different values of the fractional order α. Moreover, the fractional solution becomes closer to the classical integer-order related solution as α approaches 1. In addition, taking the new fractional operators such as the ABC derivative can present quite different and hidden aspects of the considered physical system compared to the classical Caputo derivative or the ordinary time-derivatives. In other words, the complicated behaviors of real-world dynamical phenomena can be detected more accurately through the new fractional operators.
Figure 3. Simulation results of r(t) for the 2D case when r(0) = 0.1, φ(0) = 0 and ϵ = 1. |
Figure 4. Simulation results of φ(t) for the 2D case when r(0) = 0.1, φ(0) = 0 and ϵ = 1. |
6. Conclusion
This research considered the significance of the integrals and derivatives in fractional form to investigate the motion of a mass in a system called a harmonic oscillator with a position-dependent mass. We established the Lagrangian equation in both classical and fractional form and concluded the FELEs in the concept of the newly introduced fractional operator with an ML kernel. Then, for solving the aforesaid fractional equations, we designed an efficacious approximation approach wherein the Euler convolution quadrature rule was utilized for discretizing the related convolution integral. Simulation results demonstrated that the behavior of the FELEs is modified based on the various values of the fractional order as well as the types of fractional differential operators. Moreover, the relevant classical solution was recovered as $\alpha \to 1$. Accordingly, the fractional calculus reveals the new aspects of the harmonic oscillator with mass dependent on position, which are hidden when the relative ordinary differential equations are considered. Thus, the new fractional operators can present more accurate and flexible models of real-world dynamical systems.