Application of the Adams-Bashfort-Mowlton Method to the Numerical Study of Linear Fractional Oscillators Models

Models of oscillatory systems (oscillators) are used in various fields of knowledge from mechanics to economics and biology [1, 2, 3]. From the point of view of mathematics, these models are usually described using ordinary differential equations of the second order and the corresponding initial conditions, i.e. the Cauchy problem is posed [1]. It should be noted that such mathematical models cannot take into account the properties of the environment, for example, heredity (heredity) or memory effects. This effect is characterized by the fact that the oscillating medium “remembers” the impact on it for a long time.

For the first time, the model of the hereditary oscillator was presented in his work by the Italian mathematician V. Volterra [4]. He proposed to take into account heredity in the linear oscillator model using an integro-differential equation with a difference kernel, which he later called the function of heredity or memory. It should be noted that V. Volterra derived the total energy conservation law for this generalized oscillator, in the formula of which an additional term appeared, which is responsible for the dissipation of its energy. This important result was confirmed in subsequent works on this topic.

If we choose a power-law memory function, then we can, using the mathematical apparatus of fractional calculus, go to other model equations that contain derivatives of fractional orders. In this case, the orders of fractional derivatives, as shown by the results of [5, 6], will be responsible for the intensity of energy dissipation and are related to the Q-factor of the oscillator. Oscillators with such a description are usually called fractional.

Research methods for mathematical models of fractional oscillators can be divided into exact and numerical. The exact methods, for example, include integral transformations [7] or decomposition methods [8], and numerical methods include the theory of finite-difference schemes [9], variational-iterative methods [8].

In this paper, we will carry out a numerical analysis of mathematical models of fractional linear oscillators (FLO) using elements of the theory of finite-difference schemes. As methods of numerical analysis, we will choose a method based on an explicit nonlocal finite-difference scheme studied in the author’s work [10] and the fractional ABM method, which was investigated in [11, 12, 13]. Let us analyze the errors of the methods using Runge’s rule.

Here we will consider the main definitions from the theory of fractional calculus, in more detail its aspects can be studied in the books [14, 15, 16].

Definition 1. Fractional Riemann-Liouville integral of order $\alpha$:

$$I_{0t}^{\alpha}x\left(\tau\right)=\frac{1}{{\Gamma\left(\alpha\right)}}\int\limits_{0}^{t}{\frac{{x\left(\tau\right)d\tau}}{{{{\left({t-\tau}\right)}^{1-\alpha}}}}},\alpha>0,t>0,$$

(1)

here $\Gamma\left(.\right)$  is the gamma function.

The operator (1) has the following properties:

$$I_{0t}^{0}x\left(\tau\right)=x\left(t\right),I_{0t}^{\alpha}I_{0}^{\beta}x\left(\tau\right)=I_{0t}^{\alpha+\beta}x\left(\tau\right),I_{0t}^{\alpha}I_{0}^{\beta}x\left(\tau\right)=I_{0}^{\beta}I_{0t}^{\alpha}x\left(\tau\right).$$

Definition 2. The fractional Gerasimov-Caputo derivative of order $\alpha$ has the form:

$$\partial_{0t}^{\alpha}x\left(\tau\right)=\left\{\begin{array}[]{l}\dfrac{1}{{\Gamma\left({m-\alpha}\right)}}\int\limits_{0}^{t}{\dfrac{{{x^{(m)}}\left(\tau\right)d\tau}}{{{{\left({t-\tau}\right)}^{\alpha+1-m}}}},0\leq m-1<\alpha<m,}\\ \dfrac{{{d^{m}}x\left(t\right)}}{{d{t^{m}}}},m\in N.\end{array}\right.$$

(2)

The operator (2) has the following properties[14]:

$$\partial_{0t}^{\alpha}I_{0t}^{\alpha}x\left(\tau\right)=x\left(t\right),I_{0t}^{\alpha}\partial_{0t}^{\alpha}x\left(\tau\right)-\sum\limits_{k=0}^{i-1}{\frac{{{x^{\left(k\right)}}\left(0\right){t^{k}}}}{{k!}}},t>0.$$

Consider the following Cauchy problem:

$$\partial_{0t}^{\beta}x\left(\tau\right)+\lambda\partial_{0t}^{\gamma}x\left(\tau\right)+\omega\left(t\right)x\left(t\right)=f\left(t\right),x\left(0\right)={\alpha_{1}},\dot{x}\left(0\right)={\alpha_{2}},$$

(3)

where $x\left(t\right)\in{C^{2}}\left[{0,T}\right]$ is the solution (displacement) function; $t\in\left[{0,T}\right]$ is the coordinate, responsible for time, $T>0$ is the constant, simulation time; $\omega\left(t\right)$ is the function responsible for the frequency of free oscillations and determines the type of FLO; $\lambda$ is the coefficient of friction; $f\left(t\right)$ is the function responsible for external influence; ${\alpha_{1}},{\alpha_{2}}$ is the given constants defining the initial condition, fractional differentiation operators are understood in the Gerasimov-Caputo sense (2) of orders $1<\beta<2$ and $0<\gamma<1$.

Remark 1. The Cauchy problem (3) describes a wide class of FLO and in the case of $\beta=2$ and $\gamma=1$ becomes the class of ordinary linear oscillators [1].

Definition 3. The Cauchy problem (3) will be called the fractional linear oscillator model (FLO).

Consider two numerical methods for solving the Cauchy problem (3): ENFDS and ABM method. The ABM method uses a combination of the explicit Adams-Bashforth method and the implicit Adams-Moulton method.

Explicit non-local finite-difference scheme. Let $x\left(t\right)\in C^{3}\left[0,T\right]$ to achieve the required smoothness. Divide the time interval $[0,T]$ into $N$ equal parts with a constant step $\tau={T\mathord{\left/{\vphantom{TN}}\right.\kern-1.2pt}N}$. The solution function $x\left(t\right)$ will go to the grid function $x\left({{t_{k}}}\right)={x_{k}}$, $k=1,\ldots,\,N-1$.

The Gerasimov-Caputo operators in the equation (3) are approximated as follows. For the second operator from (3):

$$\partial_{0t}^{\gamma}x\left(\tau\right)=\frac{1}{{\Gamma\left({1-\gamma}\right)}}\int\limits_{0}^{t}{\frac{{\dot{x}\left(\tau\right)d\tau}}{{{{\left({t-\tau}\right)}^{\gamma}}}}}\approx\frac{1}{{\Gamma\left({1-\gamma}\right)}}\sum\limits_{j=0}^{k}{\int\limits_{j\tau}^{\left({j+1}\right)\tau}{\frac{{\dot{x}\left(\tau\right)d\tau}}{{{{\left({{t_{k+1}}-\tau}\right)}^{\gamma}}}}}}=$$

$$=\frac{1}{{\Gamma\left({1-\gamma}\right)}}\sum\limits_{j=0}^{k}{\frac{{x\left({{t_{j+1}}}\right)-x\left({{t_{j}}}\right)}}{\tau}\int\limits_{j\tau}^{\left({j+1}\right)\tau}{\frac{{d\tau}}{{{{\left({{t_{k+1}}-\tau}\right)}^{\gamma}}}}}}=$$

(4)

$$\frac{1}{{\Gamma\left({1-\gamma}\right)}}\sum\limits_{j=0}^{k}{\frac{{x\left({{t_{j+1}}}\right)-x\left({{t_{j}}}\right)}}{\tau}\int\limits_{\left({k-j}\right)\tau}^{\left({k-j+1}\right)\tau}{\frac{{d\eta}}{{{\eta^{\gamma}}}}}}=\frac{1}{{\Gamma\left({1-\gamma}\right)}}\sum\limits_{j=0}^{k}{\frac{{x\left({{t_{k-j+1}}}\right)-x\left({{t_{k-j}}}\right)}}{\tau}\int\limits_{j\tau}^{\left({j+1}\right)\tau}{\frac{{d\eta}}{{{\eta^{\gamma}}}}}}=$$

$$=\frac{{{\tau^{-\gamma}}}}{{\Gamma\left({2-\gamma}\right)}}\sum\limits_{j=0}^{k}{\left({{{\left({j+1}\right)}^{2-\beta}}-{j^{2-\beta}}}\right)\left({x\left({{t_{k-j+1}}}\right)-x\left({{t_{k-j}}}\right)}\right)}.$$

Therefore, we have according to (4):

$$\partial_{0t}^{\gamma}x\left(\tau\right)\approx B\sum\limits_{j=0}^{k-1}{{b_{j}}\left({{x_{k-j+1}}-{x_{k-j-1}}}\right)},{b_{j}}={\left({j+1}\right)^{1-\gamma}}-{j^{1-\gamma}},B=\frac{{\lambda{\tau^{-\gamma}}}}{{\Gamma\left({2-\gamma}\right)}},$$

(5)

Similarly, for the first operator from (3) we will have:

$$\partial_{0t}^{\beta}x\left(\eta\right)\approx A\sum\limits_{j=0}^{k-1}{{a_{j}}\left({{x_{k-j+1}}-2{x_{k-j}}+{x_{k-j-1}}}\right)},{a_{j}}={\left({j+1}\right)^{2-\beta}}-{j^{2-\beta}},$$

(6)

Substituting the approximations (5) and (6) in the original equation (3), we arrive at the following discrete Cauchy problem:

$$A\sum\limits_{j=0}^{k-1}{{a_{j}}\left({{x_{k-j+1}}-2{x_{k-j}}+{x_{k-j-1}}}\right)}+B\sum\limits_{j=0}^{k-1}{{b_{j}}\left({{x_{k-j+1}}-{x_{k-j}}}\right)}+{\omega_{k}}{x_{k}}={f_{k}},$$

(7)

$${x_{0}}={\alpha_{1}},{x_{1}}={\alpha_{2}}+\tau{\alpha_{1}}.$$

For the discrete Cauchy problem (7), the following lemma holds.

Lemma 1. The coefficients of the discrete Cauchy problem (7) have the following properties:

$$1)\sum\limits_{j=0}^{k-1}{{a_{j}}={k^{2-\beta}}},\,\sum\limits_{j=0}^{k-1}{{b_{j}}={k^{1-\gamma}}},\,2)\,1=a_{0}>a_{1}>\ldots>0,1=b_{0}>b_{1}>\ldots>0,\,3)\,A\geq 0,B\geq 0$$

(8)

Proof. The first property follows from the definition:

$$\sum\limits_{j=0}^{k-1}{{a_{j}}=\sum\limits_{j=0}^{k-1}{\left[{{{\left({j+1}\right)}^{2-\beta}}-{j^{2-\beta}}}\right]}}=1-0+{2^{2-\beta}}-1+{3^{2-\beta}}-{2^{2-\beta}}+...+$$

$$+{\left({k-1}\right)^{2-\beta}}+{k^{2-\beta}}-{\left({k-1}\right)^{2-\beta}}={k^{2-\beta}}.$$

$$\sum\limits_{j=0}^{k-1}{{b_{j}}=\sum\limits_{j=0}^{k-1}{\left[{{{\left({j+1}\right)}^{1-\gamma}}-{j^{1-\gamma}}}\right]}}=1-0+{2^{1-\gamma}}-1+{3^{1-\gamma}}-{2^{1-\gamma}}+...+$$

$$+{\left({k-1}\right)^{1-\gamma}}+{k^{1-\gamma}}-{\left({k-1}\right)^{1-\gamma}}={k^{1-\gamma}}.$$

We prove the second condition as follows. Let’s introduce the following functions: $\phi\left(z\right)=(z+1)^{2-\beta}-z^{2-\beta}$, $\eta\left(z\right)=(z+1)^{1-\gamma}-z^{1-\gamma}$ where $z>0$. These functions are monotonically decreasing. Indeed, let us take the derivatives with respect to $z$ of these functions. We get:

$$\phi^{\prime}\left(z\right)=\left(2-\beta\right)(z+1)^{1-\beta}-z^{1-\beta},\,\eta^{\prime}\left(z\right)=\left(1-\gamma\right)(z+1)^{1-\gamma}-z^{1-\gamma}$$

The third property follows from the definition of the gamma function. It is known that the gamma function $\Gamma\left(z\right)>0$ for $\forall z$. Since $\tau>0,\lambda>0$, then we come to the conclusion that $A\geq 0$, $B\geq 0$. $\square$

Let us now investigate the order of approximation of fractional operators $\partial_{0t}^{\beta}x\left(\tau\right)$ and $\partial_{0t}^{\gamma}x\left(\tau\right)$. Let $\overline{\partial}_{0t}^{\beta}x\left(\tau\right)$ and $\overline{\partial}_{0t}^{\gamma}x\left(\tau\right)$ is the approximation operators. Then we have the following lemma.

Lemma 2. Approximations $\overline{\partial}_{0t}^{\beta}x\left(\tau\right)$ and $\overline{\partial}_{0t}^{\gamma}x\left(\tau\right)$ Gerasimov-Caputo operators (2) $\partial_{0t}^{\beta}x\left(\tau\right)$ and $\partial_{0t}^{\gamma}x\left(\tau\right)$ satisfy the following estimates:

$$\partial_{0t}^{\beta}x\left(\tau\right)=\overline{\partial}_{0t}^{\beta}x\left(\tau\right)+O\left(\tau^{4-\beta}\right),\partial_{0t}^{\gamma}x\left(\tau\right)=\overline{\partial}_{0t}^{\gamma}x\left(\tau\right)+O\left(\tau^{2-\gamma}\right),$$

(9)

where $O\left(\cdot\right)$ is the Landau symbol.

Proof. Using the first property (8) of Lemma 1 and the relations (5) and (6), we obtain

$$\overline{\partial}_{0t}^{\beta}x\left(\tau\right)=\frac{\tau^{2-\beta}}{\Gamma\left(3-\beta\right)}\sum\limits_{j=0}^{k-1}a_{j}\left[\ddot{x}\left(t-j\tau\right)+O\left(\tau^{2}\right)\right]=\frac{\tau^{2-\beta}}{\Gamma\left(3-\beta\right)}\sum\limits_{j=0}^{k-1}a_{j}\ddot{x}\left(t-j\tau\right)+\frac{k^{2-\beta}O\left(\tau^{4-\beta}\right)}{\Gamma\left(3-\beta\right)}=$$

$$=\frac{\tau^{2-\beta}}{\Gamma\left(3-\beta\right)}\sum\limits_{j=0}^{k-1}{{a_{j}}\ddot{x}\left({t-j\tau}\right)+O\left({{\tau^{4-\beta}}}\right)}.$$

$$\partial_{0t}^{\beta}x\left(\tau\right)=\frac{{{\tau^{2-\beta}}}}{{\Gamma\left({3-\beta}\right)}}\sum\limits_{j=0}^{k-1}{{a_{j}}\left[{\ddot{x}\left({t-{\xi_{j}}}\right)}\right]},{\xi_{j}}\in\left[{j\tau,\left({j+1}\right)\tau}\right].$$

$$|\overline{\partial}_{0t}^{\beta}x\left(\tau\right)-\partial_{0t}^{\beta}x\left(\tau\right)|=\left|\frac{\tau^{2-\beta}}{\Gamma\left(3-\beta\right)}\sum\limits_{j=0}^{k-1}a_{j}\left[\ddot{x}\left(t-j\tau\right)-\ddot{x}\left(t-\xi_{j}\right)\right]+O\left(\tau^{4-\beta}\right)\right|=$$

$$=\left|\frac{\tau^{2-\beta}}{\Gamma\left(3-\beta\right)}\sum\limits_{j=0}^{k-1}a_{j}O\left(\tau^{2}\right)+O\left(\tau^{4-\beta}\right)\right|=\left|\frac{k^{2-\beta}O\left(\tau^{4-\beta}\right)}{\Gamma\left(3-\beta\right)}+O\left(\tau^{4-\beta}\right)\right|=O\left(\tau^{4-\beta}\right).$$

Likewise, we can show the second estimate (9). The lemma is proved. $\square$

Remark 2. In the case when in relations (9) $\beta=2$ and $\gamma=1$, then we obtain approximations of the derivatives of the first and second orders with the corresponding first and second orders of approximation.

Remark 3. Note that at the internal nodes of the computational grid, taking into account Lemma 2, the explicit nonlocal finite-difference scheme (7) approximates the differential Cauchy problem with order $2-\gamma$, however, the general order of approximation of scheme (7) is the first due to the approximation of the initial conditions.

We rewrite the discrete Cauchy problem (7) as follows:

$$\left\{\begin{array}[]{l}{x_{0}}={\alpha_{1}},\,j=0,\\ {x_{1}}={\alpha_{1}}+\tau{\alpha_{2}},\,j=1,\\ {x_{k+1}}=\dfrac{{{f_{k}}+\left({2A-B-{\omega_{k}}}\right){x_{k}}-A{x_{k-1}}}}{{A+B}}-\\ -\dfrac{A}{{A+B}}\sum\limits_{j=1}^{k-1}{{a_{j}}\left({{x_{k-j+1}}-2{x_{k-j}}+{x_{k-j-1}}}\right)}-\\ -\dfrac{B}{{A+B}}\sum\limits_{j=1}^{k-1}{{b_{j}}\left({{x_{k-j+1}}-{x_{k-j}}}\right)},\,j=2,\ldots,\,N-1.\end{array}\right.$$

(10)

or in matrix form:

$${X_{k+1}}=M{X_{k}}+{F_{k}},$$

(11)

$${X_{k+1}}={\left({{x_{1}},{x_{2}},\ldots,{x_{k}}}\right)^{T}},\,{X_{k}}={\left({{x_{0}},{x_{1}},\ldots,{x_{k-1}}}\right)^{T}},$$

$${F_{k}}={\left({{f_{0}},\frac{{{f_{1}}}}{{A+{\lambda_{1}}B}},\ldots,\frac{{{f_{k-1}}}}{{A+{\lambda_{k-1}}B}}}\right)^{T}},\,{f_{0}}=\tau{\alpha_{2}},$$

where the Hessenberg matrix $M$ in (11) has the form:

$${m_{ij}}=\left\{\begin{array}[]{l}0,\,j\geq i+1,\\ \dfrac{{A\left({2-{a_{1}}}\right)-B{\lambda_{i-1}}{b_{1}}}}{{A+{\lambda_{i-1}}B}},\,j=i=3,\ldots,\,N-1,\\ \dfrac{{-A\left({{a_{i-j+1}}-2{a_{i-j}}+{a_{i-j-1}}}\right)-B{\lambda_{i-1}}\left({{b_{i-j+1}}-{b_{i-j-1}}}\right)}}{{A+{\lambda_{i-1}}B}},\,j\leq i-1,\end{array}\right.$$

(12)

$${m_{1,1}}=1,\,{m_{2,2}}=\dfrac{{2A}}{{A+{\lambda_{1}}B}},\,{m_{i,1}}=\dfrac{{B{\lambda_{i-1}}{b_{i-2}}-A{a_{i-2}}}}{{A+{\lambda_{i-1}}B}},\,i=2,\ldots,\,N-1,$$

$${m_{i,2}}=\dfrac{{A\left({2{a_{i-2}}-{a_{i-3}}}\right)+{\lambda_{i-1}}B{b_{i-3}}}}{{A+{\lambda_{i-1}}B}},\,i=3,\ldots,\,N-1.$$

The following theorems are true.

Theorem 1 [10]. A nonlocal explicit finite-difference scheme (10) converges with the first order if the following condition is satisfied:

$$\tau\leq{\tau_{0}}=\min\left({1,{{\left({\dfrac{{\Gamma\left({2-\gamma}\right)}}{{\lambda\Gamma\left({3-\beta}\right)}}}\right)}^{\dfrac{1}{{\beta-\gamma}}}}}\right).$$

(13)

Let ${X_{k}},{Y_{k}}$ be two different solutions of the matrix equation (11) with the initial conditions ${X_{0}},{Y_{0}}$. Then the scheme stability theorem is valid.

Theorem 2 [10]. A nonlocal explicit finite-difference scheme (10) is conditionally stable if the condition (13) is satisfied and the estimate $\left|{{Y_{k}}-{X_{k}}}\right|\leq C\left|{{Y_{0}}-{X_{0}}}\right|$ for any $k$, where $C>0$ is a constant independent of the step $\tau$.

The proof of Theorems 1 and 2 is based on the results of Lemmas 1 and 2.

Adams-Bashfort-Moulton method. The ABM method is a type of numerical predictor-corrector method for solving differential equations. It has been studied and discussed in detail in the papers [11, 12, 13]. Let’s generalize this method for solving the Cauchy problem (3). Taking into account Definitions 1 and 2, as well as the corresponding properties of the operators of fractional integration and differentiation, we write it in the form of a system:

$$\left\{\begin{array}[]{l}\partial_{0t}^{{\sigma_{1}}}x\left(\tau\right)=y\left(t\right),\\ \partial_{0t}^{{\sigma_{2}}}y\left(\tau\right)=f\left(t\right)-\lambda y\left(t\right)-\omega\left(t\right)x\left(t\right),\\ x\left(0\right)={\alpha_{1}},y\left(0\right)={\alpha_{2}}\end{array}\right.$$

(14)

where ${\sigma_{1}}=\gamma$, ${\sigma_{2}}=\beta-\gamma$.

Let $0<{\sigma_{2}}\leq 1$, i.e. $\left\{\beta\right\}<\left\{\gamma\right\}$,   $\left\{\cdot\right\}$ is the fractional part of the number.

Remark 4. Note that if the system (14) $\left\{\beta\right\}>\left\{\gamma\right\}$, i.e. $1<{\sigma_{2}}<2$, then we can always transform it by adding one more equation.

On a uniform grid, we introduce the functions $x_{n+1}^{p},\,y_{n+1}^{p}$, $n=0,\ldots,\,N-1$, which will be determined by the Adams-Bashforth formula (predictor):

$$\left\{\begin{array}[]{l}x_{n+1}^{p}={x_{0}}+\dfrac{{{\tau^{{\sigma_{1}}}}}}{{\Gamma\left({{\sigma_{1}}+1}\right)}}\sum\limits_{j=0}^{n}{\theta_{j,n+1}^{1}{y_{j}},}\\ y_{n+1}^{p}={y_{0}}+\dfrac{{{\tau^{{\sigma_{2}}}}}}{{\Gamma\left({{\sigma_{2}}+1}\right)}}\sum\limits_{j=0}^{n}{\theta_{j,n+1}^{2}\left({-\lambda{y_{j}}-{\omega_{j}}{x_{j}}+{f_{j}}}\right),}\\ \theta_{j,n+1}^{i}={\left({n-j+1}\right)^{{\sigma_{i}}}}-{\left({n-j}\right)^{{\sigma_{i}}}},i=1,2.\end{array}\right.$$

(15)

Then, using the Adams-Moulton formula for the corrector, we get:

$$\left\{\begin{array}[]{l}{x_{n+1}}={x_{0}}+\dfrac{{{\tau^{{\sigma_{1}}}}}}{{\Gamma\left({{\sigma_{1}}+2}\right)}}\left({y_{n+1}^{p}+\sum\limits_{j=0}^{n}{\rho_{j,n+1}^{1}{y_{j}}}}\right),\\ {y_{n+1}}={y_{0}}+\dfrac{{{\tau^{{\sigma_{2}}}}}}{{\Gamma\left({{\sigma_{2}}+2}\right)}}\left({-\lambda y_{n+1}^{p}-{\omega_{n+1}}x_{n+1}^{p}+{f_{n+1}}+\sum\limits_{j=0}^{n}{\rho_{j,n+1}^{2}\left({-\lambda{y_{j}}-{\omega_{j}}{x_{j}}+{f_{j}}}\right)}}\right),\end{array}\right.$$

(16)

where the weight factors in (16) are determined by the formula:

$$\rho_{j,n+1}^{i}=\left\{\begin{array}[]{l}{n^{{\sigma_{i}}+1}}-\left({n-{\sigma_{i}}}\right){\left({n+1}\right)^{{\sigma_{i}}}},\,j=0,\\ {\left({n-j+2}\right)^{{\sigma_{i}}+1}}+{\left({n-j}\right)^{{\sigma_{i}}+1}}-2{\left({n-j+1}\right)^{{\sigma_{i}}+1}},\,1\leq j\leq n,\\ 1,\,j=n+1,\\ i=1,2.\end{array}\right.$$

Theorem 3 [12]. If $\partial_{0t}^{{\sigma_{i}}}{x_{i}}\left(\tau\right)\in{C^{2}}\left[{0,T}\right],\left({{x_{1}}=x\left(t\right),\,{x_{2}}=y\left(t\right),\,i=1,2}\right)$, then

$$\mathop{\max}\limits_{1\leq j\leq n}\left|{{x_{i}}\left({{t_{j}}}\right)-{x_{i,j}}}\right|=O\left({{h^{1+\mathop{\min}\limits_{i}{\sigma_{i}}}}}\right).$$

(17)

Let us examine how the computational accuracy of the methods behaves. To do this, we will use the double recalculation method (Runge’s rule) to estimate the error using the formula:

$$\xi=\frac{{\mathop{\max}\limits_{i}\left|{{x_{i}}-{x_{2i}}}\right|}}{{{2^{\mu}}-1}},$$

(18)

where $\mu$ is the order of approximation of the numerical method, ${x_{2i}}$ is the numerical solution at the step ${\tau\mathord{\left/{\vphantom{\tau 2}}\right.\kern-1.2pt}2}$. The computational accuracy of $p$ is determined from the formula:

$$p={\log_{\frac{{{\tau_{1}}}}{{{\tau_{2}}}}}}\frac{{{\xi_{1}}}}{{{\xi_{2}}}},$$

(19)

${\tau_{1}},\,{\tau_{2}}={{{\tau_{1}}}\mathord{\left/{\vphantom{{{\tau_{1}}}2}}\right.\kern-1.2pt}2}$   — grid steps, ${\xi_{1}},\,{\xi_{2}}$   — errors at step ${\tau_{1}}$ and at step ${\tau_{2}}$.

Example 1. Consider the Cauchy problem (3) with homogeneous initial conditions, choosing the following functions and parameter values: $\omega\left(t\right)=\omega_{0}^{\beta}$,  $f\left(t\right)=\omega_{0}^{\beta}{t^{3}}+\dfrac{{6{t^{3-\beta}}}}{{\Gamma\left({4-\beta}\right)}}+\dfrac{{6\lambda{t^{3-\gamma}}}}{{\Gamma\left({4-\gamma}\right)}},\,\lambda=0.1,\,\beta=1.8,\,\gamma=0.9,\,{\omega_{0}}=1,\,t\in\left[{0,1}\right]$.

$$\partial_{0t}^{1.8}x\left(\tau\right)+0.1\partial_{0t}^{0.9}x\left(\tau\right)+x\left(t\right)={t^{3}}+\frac{{6{t^{1.2}}}}{{\Gamma\left({2.2}\right)}}+\frac{{0.6{t^{2.1}}}}{{\Gamma\left({3.1}\right)}},x\left(0\right)=\dot{x}\left(0\right)=0.$$

(20)

The solution to the Cauchy problem (20), as you can see, is the function:

$$x\left(t\right)={t^{3}}.$$

(21)

Remark 5. The values of the parameters were chosen in such a way that the condition of Theorems 1 and 2 would be satisfied.

Due to the fact that the exact solution (21) of the Cauchy problem (20) is known, we can calculate the errors of numerical methods from the following relations:

$$\xi_{F}^{EX}=\mathop{\max}\limits_{j}\left({\left|{x\left({{t_{j}}}\right)-x_{j}^{F}}\right|}\right),\xi_{PC}^{EX}=\mathop{\max}\limits_{j}\left({\left|{x\left({{t_{j}}}\right)-x_{j}^{PC}}\right|}\right),$$

(22)

where $\xi_{F}^{EX},\xi_{PC}^{EX}$   — errors of ENFDS and ABM methods, orders of computational accuracy   — $p_{F}$ and $p_{PC}$, which calculated by the formula (19).

From Table 1 shows that with an increase in the nodes of the computational grid, the errors of numerical methods decrease. It can be seen that the ABM method converges faster than the JANCRS method. Based on (9), the ENFDS approximation is of the first order. According to the relation (17) for Example 1, the order of the ABM method.

From Table 1, we also see that with an increase in the nodes of the computational grid, the computational accuracy increases and tends to unity, and $p_{PC}$   — almost immediately reached the value 1.9. On the one hand, this confirms the validity of Lemma 3 and Theorem 3, and on the other hand, the faster convergence and high accuracy of the ABM method in comparison with the ENFDS method.

As a comparison with the results from Table 1, we calculate the errors and orders of magnitude of the computational accuracy of numerical methods taking into account the Runge rule (Table 2).

From Table 2 that the double recalculation method based on the Runge rule is in good agreement with the results from Table 1.

In the classical case, when the ABM method has the second order, and the ENFDS is still the same first order of approximation. Numerical analysis of errors for Example 1 is given in Table 3.

From Table 3 that the orders of computational accuracy and tend to the orders of approximation 1 and 2 for the corresponding numerical methods.

Figure 1 illustrates the effectiveness of the ABM method in comparison with the ENFDS. It can be seen that, due to the fast convergence and higher accuracy, the ABM method gives the best result.

However, the ABM method does not always work better than explicit methods like ENFDS. Let us show this with the following example.

Example 2. Consider an analogue of the harmonic oscillator (AHO). To do this, in the Cauchy problem (3) choose:

$$\omega\left(t\right)=\omega_{0}^{\beta},f\left(t\right)=\frac{{\mu{t^{\delta}}}}{{\Gamma\left({1+\delta}\right)}},\lambda=0.$$

Then we come to the following problem:

$$\partial_{0t}^{\beta}x\left(\tau\right)+\omega_{0}^{\beta}x\left(t\right)=\frac{{\mu{t^{\delta}}}}{{\Gamma\left({1+\delta}\right)}},x\left(0\right)={\alpha_{1}},\dot{x}\left(0\right)={\alpha_{2}}.$$

(23)

The solution to the Cauchy problem (23) is the function:

$$x\left(t\right)={\alpha_{1}}{E_{\beta,1}}\left({-{{\left({{\omega_{0}}t}\right)}^{\beta}}}\right)+{\alpha_{2}}t{E_{\beta,2}}\left({-{{\left({{\omega_{0}}t}\right)}^{\beta}}}\right)+\mu{t^{\beta+g}}{E_{\beta,\,\beta+\delta+1}}\left({-{{\left({{\omega_{0}}t}\right)}^{\beta}}}\right),\delta>-\beta,$$

(24)

where ${E_{\alpha,\beta}}\left(z\right)=\sum\limits_{k=0}^{\infty}{\dfrac{{{z^{k}}}}{{\Gamma\left({\alpha k+\beta}\right)}}}$ — two-parameter Mittag-Leffler function [14].

Let’s consider the numerical methods ABM and ENFDS for solving the Cauchy problem (23) and compare them with the exact solution (24). First, consider the usual harmonic oscillator ($\beta=2$).

The parameter values in the problem (24) are chosen as follows: $\delta=0.3,\mu=0.5,{\omega_{0}}=2,{\alpha_{1}}={\alpha_{2}}=1,t\in\left[{0,1}\right]$.

We see in Figure 2 that the ABM method works in much the same way as the ENFDS method. The order of computational accuracy is shown in the following Table 4. In the case of AHO, the error analysis is given in Table 5.

We see that the order of the computational accuracy of the methods tends to unity, i.e. the accuracy of the ABM method is similar to that of the ENFDS method, although the convergence rate is higher.

In this work, a numerical analysis of the Cauchy problem for a model class of fractional linear oscillators was carried out. The methods of ENFDS and ABM were considered as numerical methods. It has been shown on test examples that the ABM method converges faster and is more accurate (not always) than the ENFDS method.

The study of FLO with variable memory is also of interest. In this case, derivatives of fractional variable orders appear in the model equation of the Cauchy problem (3). Some aspects of the numerical analysis of such a generalized Cauchy problem (3) were considered by the author in the work [17].