Mathematical Models of Oscillators with Memory

Abstract

The chapter proposes a mathematical model for a wide class of hereditary oscillators, which is a Cauchy problem in the local formulation. As an initial model equation, an integrodifferential equation of Voltaire type was introduced, which was reduced by means of a special choice of difference kernels to a differential equation with nonlocal derivatives of fractional-order variables. An explicit finite-difference scheme is proposed, and questions of its stability and convergence are investigated. A computer study of the proposed numerical algorithm on various test examples of the hereditary oscillators Airy, Duffing, and others was carried out. Oscillograms and phase trajectories are plotted and constructed.

1. Introduction

In the paper of the Italian mathematician Vito Volterra [1], the notion of heredity (memory), a property of a dynamical system characterized by nonlocality in time, is introduced, which consists in the dependence of its current state on a finite number of previous states. In another paper [2], Volterra investigated the hereditary oscillator—a vibration system with memory, which was written in the form of an integrodifferential equation with a difference kernel, a function of memory. Further, for such an oscillator, Volterra derived the law of total energy, in which an additional term appeared, responsible for the dissipation of energy in the vibrational system. This fact was confirmed in subsequent works.

In papers [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], fractal oscillators were considered, which represent the class of hereditary oscillators with a power-law function of memory. The peculiarity of such oscillators is that their mathematical description can be reduced to differential equations with nonlocal derivatives of fractional constant orders, which are investigated within the framework of the theory of fractional calculus [22, 23, 24].

In papers [6, 7, 8, 9, 11, 12, 13, 14, 18, 19, 20, 21], models of fractal linear oscillators were investigated in the sense of the Gerasimov-Caputo derivative and in papers [4, 10, 16]—in the sense of the Riemann-Liouville derivative. Analytical solutions of model equations in terms of a special function of Mittag-Leffler-type and generalized Wright-type function, oscillograms, and phase trajectories are constructed. It is shown that in the regime of free oscillations, the presence of memory effects in the system leads to attenuation of oscillations as a result of energy dissipation, and with allowance for external periodic action, it is possible to stabilize the amplitude of the oscillations, with the phase trajectories reaching the limit cycle and also the resonance effect.

In papers of the author [25, 26], the fractal parametric resonance (the fractal Mathieu oscillator) was investigated, and the Strutt-Ince diagrams of parametric resonance existence areas were constructed. It is shown that these regions strongly depend on the orders of the value fractional derivatives entering into the initial equation.

In a monograph by the Slovak mathematician Ivo Petras [10], the fractal nonlinear oscillator models whose differential equations contained fractional derivatives in the sense of Riemann-Liouville were considered and analyzed using numerical methods and considered the stability of the rest point of oscillatory systems. However, the stability and convergence of numerical methods have not been considered.

A further continuation of the investigation of hereditary oscillators is associated with the introduction of the derivatives of fractional variable orders in the model equations. This is due to the fact that the orders of fractional derivatives are related to the properties of the medium in which this or that process takes place and changes with time under the influence of external influence. Therefore, papers [27, 28, 29, 30] proposed that the models of fractal nonlinear oscillators (Duffing, Van der Pol, Van der Pol-Duffing, FitzHugh-Nagumo) were proposed and investigated using explicit finite-difference schemes, whose equations contain both the derivatives of the constants and variable fractional orders of the Gerasimov-Caputo and Riemann-Liouville types. With the help of computer experiments, the convergence of finite-difference schemes was shown, and estimates of the computational accuracy of the method were obtained; oscillograms and phase trajectories were constructed. However, the questions of stability and convergence were not formulated in the form of corresponding theorems.

In [31, 32], a new class of fractal oscillators was proposed and investigated; active fractal oscillators (AFOs)—nonlinear oscillators with external influences, which include the fractional Riemann-Liouville integral—were investigated. Such oscillators are constructed on the basis of the scheme of the radioelectronic аutogenerator with a fractional feedback circuit. Authors use the method of equivalent linearization to investigate AFOs and come to the conclusion that the self-oscillator is isochronous.

From the analysis of the above publications on the study of hereditary oscillator, we can conclude that the main tool for their study is numerical methods, for example, finite-difference schemes. In most cases, the authors leave without considering the questions of stability and convergence of finite-difference schemes and, even if they touch, then without formulating the corresponding theorems and proofs. Therefore, the goal of the present paper is to construct a finite-difference scheme for a wide class of hereditary (fractal) linear and nonlinear oscillators, prove its stability and convergence, formulate results in the form of corresponding theorems, and study finite-difference schemes on specific test examples.

2. Formulation of the problem

Consider the following model integrodifferential equation for the function $x\left(t\right)\in C^3\left(0,T\right)$ , where $T>0$ :

$$\underset{0}{\overset{t}{\int }}{K}_1\left(t-\eta \right)\ddot{x}\left(\eta \right)\mathit{d\eta }+\lambda \underset{0}{\overset{t}{\int }}{K}_2\left(t-\eta \right)\dot{x}\left(\eta \right)\mathit{d\eta }=f\left(x\left(t\right),t\right),$$ E1

where $\ddot{x}\left(t\right)=d^{2}x/{\mathit{dt}}^2,\dot{x}\left(t\right)=\mathit{dx}/\mathit{dt}$ , $\lambda >0$ —given constant.

Eq. (1) describes a wide class of hereditary, depending on the form of the right-hand side (function $f\left(x\left(t\right),t\right)$ ) of linear or nonlinear oscillators.

Definition 1. Functions $K_1\left(t-\eta \right)$ and $K_2\left(t-\eta \right)$ —difference kernels in Eq. (1)—will be called memory functions, since they define the notion of heredity (memory), which was introduced in the work of the Italian mathematician Vito Volterra [2].

Definition 2. A nonlinear function $f\left(x\left(t\right),t\right)$ on the right-hand side of Eq. (1) satisfies a Lipschitz condition with respect to a variable $x\left(t\right)$ :

$$\left|f\left(x_1\left(t\right),t\right)-f\left(x_2\left(t\right),t\right)\right|\le L\left|x_1\left(t\right)-x_2\left(t\right)\right|,$$ E2

$L$ —Lipschitz constant.

Eq. (1) describes a broad class of hereditary nonlinear oscillators, depending on the form of the function $f\left(x\left(t\right),t\right)$ on its right-hand side, and the parameter $\lambda$ has the meaning of the coefficient of friction.

Note that in Definition 1, the memory functions $K_1\left(t-\eta \right)$ and $K_2\left(t-\eta \right)$ can be chosen arbitrarily, depending on the conditions of the particular problem. We will choose these functions power law, since power laws are often found in various fields of knowledge [33]. We choose the memory functions $K_1\left(t-\eta \right)$ and $K_2\left(t-\eta \right)$ in the form

$$K_1\left(t-\eta \right)=\frac{{\left(t-\eta \right)}^{1-\beta \left(t\right)}}{\mathrm{\Gamma }\left(2-\beta \left(t\right)\right)},K_2\left(t-\eta \right)=\frac{{\left(t-\eta \right)}^{-\gamma \left(t\right)}}{\mathrm{\Gamma }\left(2-\gamma \left(t\right)\right)},1<\beta \left(t\right)<2,0<\gamma \left(t\right)<1,$$ E3

where $\gamma \left(t\right),\beta \left(t\right)\in C\left[0,T\right]$ , $\mathrm{\Gamma }\left(t\right)$ —Euler gamma function.

We give the following definitions.

Definition 3. Derivatives of fractional variables of orders $\beta \left(t\right)$ and $\gamma \left(t\right)$ Gerasimov-Caputo type: we call the following operators of fractional differentiation:

$${\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)=\frac{1}{\mathrm{\Gamma }\left(2-\beta \left(t\right)\right)}\underset{0}{\overset{t}{\int }}\frac{\ddot{x}\left(\eta \right)\mathit{d\eta }}{{\left(t-\eta \right)}^{\beta \left(t\right)-1}},{\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)=\frac{1}{\mathrm{\Gamma }\left(1-\gamma \left(t\right)\right)}\underset{0}{\overset{t}{\int }}\frac{\dot{x}\left(\eta \right)\mathit{d\eta }}{{\left(t-\eta \right)}^{\gamma \left(t\right)}}.$$ E4

We note that in the case when the functions $\beta \left(t\right)$ and $\gamma \left(t\right)$ in the relations (4) and (5) are constants, we arrive at the definitions of the fractional derivative in the sense of Gerasimov-Caputo [34, 35], and in the case when these constants $\beta =2$ and $\gamma =1$ , the operators of fractional differentiation (4) become classical derivatives of the second and first orders.

Taking into account Definition 3, the model Eq. (1) can be rewritten in a more compact form:

$${\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)+\lambda {\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)=f\left(x\left(t\right),t\right).$$ E5

For Eq. (5), the initial conditions in the local formulation are valid:

$$x\left(0\right)={\alpha }_0,\dot{x}\left(0\right)={\alpha }_1,$$ E6

where ${\alpha }_0$ and ${\alpha }_1$ —given constants. As a result, we arrive at the Cauchy problems (5) and (6), which we will investigate.

3. Explicit finite-difference scheme

We construct an explicit finite-difference scheme for the Cauchy problems (5) and (6). We divide the time interval $\left[0,T\right]$ into $N$ equal parts with a step $\tau =T/N$ . We introduce the grid function $x\left(t_k\right)=x_k$ , where $t_k=\mathit{k\tau },k=1,\dots ,N-1$ .

The derivatives of the variables of fractional orders in Eq. (5) are approximated according to the relations in [36, 37]:

$${\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)\approx A_k\sum _{j=0}^{k-1}{a}_{j,k}\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right),$$ E7

$${\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)\approx B_k\sum _{j=0}^{k-1}{b}_{j,k}\left(x_{k-j+1}-x_{k-j-1}\right),$$

then the formulas will refer to the formula (7)

$$a_{j,k}={\left(j+1\right)}^{2-{\beta }_k}-j^{2-{\beta }_k},b_{j,k}={\left(j+1\right)}^{1-{\gamma }_k}-j^{1-{\gamma }_k},$$

$$A_k=\frac{{\tau }^{-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)},B_k=\frac{{\mathit{\lambda \tau }}^{-{\gamma }_k}}{2\mathrm{\Gamma }\left(2-{\gamma }_k\right)},$$

Here, to shorten the record, $\beta \left(t_k\right)={\beta }_k,\gamma \left(t_k\right)={\gamma }_k$ .

Taking into account relation (7), the Cauchy problems (5) and (6) in the difference formulation will have the form

$$\begin{array}{c}{A}_k\sum _{j=0}^{k-1}{a}_{j,k}\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right)+B_k\sum _{j=0}^{k-1}{b}_{j,k}\left(x_{k-j+1}-x_{k-j-1}\right)=f_k,\\ x_0={\alpha }_0,x_1={\alpha }_1+{\mathit{\tau \alpha }}_0,\end{array}$$ E8

Here, to shorten the record, $f_k=f\left(x_k,t_k\right)$ . We write the problem (8) explicitly:

$$\begin{array}{c}{x}_{k+1}=\frac{1}{A_k+B_k}\left(2A_k{x}_k-\left(A_k-B_k\right)x_{k-1}-A_k\sum _{j=1}^{k-1}{a}_{j,k}\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right)\right)\\ -\frac{B_k}{A_k+B_k}\sum _{j=1}^{k-1}{b}_{j,k}\left(x_{k-j+1}-x_{k-j-1}\right)+f_k.\end{array}$$ E9

We note that the weight coefficients $a_{j,k}$ and $b_{j,k}$ have properties, which we formulate in the form of the following lemmas.

Lemma 1. For any $k$ weights, coefficients $a_{j,k}$ and $b_{j,k}$ , as well as coefficients $A_k$ and $B_k$ , have the following properties:

1. $\sum _{j=0}^{k-1}{a}_{j,k}=k^{2-{\beta }_k},\sum _{j=0}^{k-1}{b}_{j,k}=k^{1-{\gamma }_k},$

2. $1=a_{0,k}>a_{1,k}>\dots >0,1=b_{0,k}>b_{1,k}>\dots >0,$

3. $A_k\ge 0,B_k\ge 0.$

Proof. The first property follows the definition of weight coefficients:

$\begin{array}{l}\sum _{j=0}^{k-1}{a}_{j,k}=\sum _{j=0}^{k-1}\left[{\left(j+1\right)}^{2-{\beta }_k}-j^{2-{\beta }_k}\right]=1-0+2^{2-{\beta }_k}-1+3^{2-{\beta }_k}-2^{2-{\beta }_k}+.\dots +{\left(k-1\right)}^{2-{\beta }_k}+k^{2-{\beta }_k}\\ -{\left(k-1\right)}^{2-{\beta }_k}=k^{2-{\beta }_k}.\end{array}$

$\begin{array}{l}\sum _{j=0}^{k-1}{b}_{j,k}=\sum _{j=0}^{k-1}\left[{\left(j+1\right)}^{1-{\gamma }_k}-j^{1-{\gamma }_k}\right]=1-0+2^{1-{\gamma }_k}-1+3^{1-{\gamma }_k}-2^{1-{\gamma }_k}+.\dots +{\left(k-1\right)}^{1-{\gamma }_k}+k^{1-{\gamma }_k}\\ -{\left(k-1\right)}^{1-{\gamma }_k}=k^{1-{\gamma }_k}.\end{array}$

The second property is proven in the following way. We introduce two functions:

$$\phi \left(x\right)={\left(x+1\right)}^{2-{\beta }_k}-x^{2-{\beta }_k}\text{and}\eta \left(x\right)={\left(x+1\right)}^{1-{\gamma }_k}-x^{1-{\gamma }_k},$$

where $x>0$ . These functions are decreasing. Really derived from these functions

$${\phi }^{\prime }\left(x\right)=\left(2-{\beta }_k\right)\left[{\left(x+1\right)}^{1-{\beta }_k}-x^{1-{\beta }_k}\right]<0,{\eta }^{\prime }\left(x\right)=\left(1-{\gamma }_k\right)\left[{\left(x+1\right)}^{1-{\gamma }_k}-x^{1-{\gamma }_k}\right]<0.$$

Therefore, the second property holds. The third property follows also the properties of the gamma function. The lemma is proven.

Let ${\overline{\partial }}_{0t}^{\beta \left(t\right)}x\left(\eta \right)$ and ${\overline{\partial }}_{0t}^{\gamma \left(t\right)}x\left(\eta \right)$ —approximations of differential operators of Gerasimov-Caputo types ${\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)$ and ${\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)$ . Then, we have the following lemma.

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

$$\left|{\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)-{\overline{\partial }}_{0t}^{\beta \left(t\right)}x\left(\eta \right)\right|\le C_1\tau ,\left|{\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)-{\overline{\partial }}_{0t}^{\gamma \left(t\right)}x\left(\eta \right)\right|\le C_2\tau ,$$ E10

where ${С}_1$ and ${С}_2$ —constants that are independent of the parameter $\mathit{\tau }$ .

Proof. Using the first property of Lemma 1 and Definition 3, we obtain

$$\begin{array}{c}{\overline{\partial }}_{0t}^{\beta \left(t\right)}x\left(\eta \right)=\frac{{\tau }^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}\sum _{j=0}^{k-1}{a}_{j,k}\left[\ddot{x}\left(t-\mathit{j\tau }\right)+O\left({\tau }^2\right)\right]=\frac{{\tau }^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}\sum _{j=0}^{k-1}{a}_{j,k}\ddot{x}\left(t-\mathit{j\tau }\right)\\ +\frac{{\tau }^{2-{\beta }_k}{k}^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}O\left({\tau }^2\right)=\frac{{\tau }^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}\sum _{j=0}^{k-1}{a}_{j,k}\ddot{x}\left(t-\mathit{j\tau }\right)+\frac{t^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}O\left({\tau }^2\right)=\frac{{\tau }^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}\\ \sum _{j=0}^{k-1}{a}_{j,k}\ddot{x}\left(t-\mathit{j\tau }\right)+O\left({\tau }^2\right).\end{array}$$

$${\partial }_{0t}^{\beta }x\left(\eta \right)=\frac{1}{\mathrm{\Gamma }\left(2-{\beta }_k\right)}\sum _{j=0}^{k-1}\underset{\mathit{j\tau }}{\overset{\left(j+1\right)\tau }{\int }}{\xi }^{1-{\beta }_k}\ddot{x}\left(t-\xi \right)\mathit{d\xi }=\frac{1}{\mathrm{\Gamma }\left(2-{\beta }_k\right)}\sum _{j=0}^{k-1}{a}_{j,k}\ddot{x}\left(t-{\eta }_j\right),{\eta }_j\in \left[\mathit{j\tau },\left(j+1\right)\tau \right].\left|{\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)-{\overline{\partial }}_{0t}^{\beta \left(t\right)}x\left(\eta \right)\right|=\left|\frac{{\tau }^{2-\beta \left(t\right)}}{\mathrm{\Gamma }\left(3-\beta \left(t\right)\right)}\sum _{j=0}^{k-1}{a}_{j,k}\left[\ddot{x}\left(t-\mathit{j\tau }\right)-\ddot{x}\left(t-{\eta }_j\right)\right]+O\left({\tau }^2\right)\right|=\left|\frac{{\tau }^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}\sum _{j=0}^{k-1}{a}_{j,k}\cdot O\left(\tau \right)+O\left({\tau }^2\right)\right|=\left|\frac{{\tau }^{2-{\beta }_k}{k}^{2-{\beta }_k}}{\mathrm{\Gamma }\left(3-{\beta }_k\right)}O\left(\tau \right)+O\left({\tau }^2\right)\right|=O\left(\tau \right)+O\left({\tau }^2\right)=O\left(\tau \right).$$

Similarly, we can obtain the second estimate from Eq. (10).

$${\overline{\partial }}_{0t}^{\gamma \left(t\right)}x\left(\eta \right)=\frac{{\tau }^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\left[\dot{x}\left(t-\mathit{j\tau }\right)+O\left(\tau \right)\right]=\frac{{\tau }^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\dot{x}\left(t-\mathit{j\tau }\right)+\frac{{\tau }^{1-{\gamma }_k}{k}^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}$$

$$O\left(\tau \right)=\frac{{\tau }^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\dot{x}\left(t-\mathit{j\tau }\right)+\frac{t^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}O\left({\tau }^2\right)=\frac{{\tau }^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\dot{x}\left(t-\mathit{j\tau }\right)+O\left(\tau \right).$$

$${\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)=\frac{1}{\mathrm{\Gamma }\left(1-{\gamma }_k\right)}\sum _{j=0}^{k-1}\underset{\mathit{j\tau }}{\overset{\left(j+1\right)\tau }{\int }}{\xi }^{-{\gamma }_k}\dot{x}\left(t-\xi \right)\mathit{d\xi }=\frac{1}{\mathrm{\Gamma }\left(1-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\dot{x}\left(t-{\eta }_j\right),{\eta }_j\in \left[\mathit{j\tau },\left(j+1\right)\tau \right].$$

$$\left|{\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)-{\overline{\partial }}_{0t}^{\gamma \left(t\right)}x\left(\eta \right)\right|=\left|\frac{{\tau }^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\left[\dot{x}\left(t-\mathit{j\tau }\right)-\dot{x}\left(t-{\eta }_j\right)\right]+O\left(\tau \right)\right|=$$

$$=\left|\frac{{\tau }^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}\sum _{j=0}^{k-1}{b}_{j,k}\cdot O\left(\tau \right)+O\left(\tau \right)\right|=\left|\frac{{\tau }^{1-{\gamma }_k}{k}^{1-{\gamma }_k}}{\mathrm{\Gamma }\left(2-{\gamma }_k\right)}O\left(\tau \right)+O\left(\tau \right)\right|=O\left(\tau \right)+O\left(\tau \right)=O\left(\tau \right).$$

The lemma is proven.

Investigation. According to Lemma 2, it can be shown that the explicit finite-difference scheme (9) has an error $\epsilon =O\left(\tau \right)$ . This fact will be used in computer experiments in determining the computational accuracy of the numerical method.

Lemma 3. The sums in the finite-difference scheme (9) have the following representations:

$$\begin{array}{l}\sum _{j=1}^{k-1}{a}_{j,k}\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right)=a_{1,k}{x}_k+a_{k-1,k}{x}_0\\ +\left(a_{k-2,k}-2a_{k-1,k}\right)x_1+\left(a_{2,k}-2a_{1,k}\right)x_{k-1}+\sum _{j=2}^{k-2}\left(a_{j+1,k}-2a_{j,k}+a_{j-1,k}\right)x_{k-j},\\ \sum _{j=1}^{k-1}{b}_{j,k}\left(x_{k-j+1}-x_{k-j-1}\right)=b_{1,k}{x}_k-b_{k-1,k}{x}_0+b_{2,k}{x}_{k-1}-b_{k-2,k}{x}_1\\ +\sum _{j=2}^{k-2}\left(b_{j+1,k}-b_{j-1,k}\right)x_{k-j}.\end{array}$$ E11

Proof. The representation (11) follows the properties of the sum. Indeed, by opening the sums in Eq. (9) and grouping the terms properly, we arrive at the representation (11).

Using Lemma 3, the finite-difference scheme (9) can be rewritten as

$$x_1={\alpha }_0+{\mathit{\tau \alpha }}_1,$$

$$x_2=\frac{2A_1}{A_1+B_1}{x}_1-\frac{A_1-B_1}{A_1+B_1}{x}_0+\frac{f_1}{A_1+B_1},k=1,$$

$$\begin{array}{l}{x}_{k+1}=\frac{1}{A_k+B_k}\\ \left[\left(A_k\left(2-a_{1,k}\right)-B_k{b}_{1,k}\right)x_k-\left(A_k\left(a_{2,k}-2a_{1,k}+a_{0,k}\right)+B_k\left(b_{2,k}-b_{0,k}\right)\right)x_{k-1}+f_k\right]-\\ -\frac{\left(A_k\left(a_{k-2,k}-2a_{k-1,k}\right)-B_k{b}_{k-2,k}\right)x_1-\left(A_k{a}_{k-1,k}-B_k{b}_{k-1,k}\right)x_0}{A_k+B_k}\\ -\frac{1}{A_k+B_k}\sum _{j=2}^{k-2}\left[A_k\left(a_{j+1,k}-2a_{j,k}+a_{j-1,k}\right)+B_k\left(b_{j+1,k}-b_{j-1,k}\right)\right]x_{k-j},k=2,\dots ,n-1,\end{array}$$

Or in matrix form

$$X_{k+1}={\mathit{MX}}_k+F_k,$$ E12

$$X_{k+1}={\left(x_1,x_2,\dots ,x_{N-1}\right)}^T,X_k={\left(x_0,x_1,\dots ,x_{N-2}\right)}^T,F_k={\left(f_0,f_1,\dots ,f_{N-2}\right)}^T$$

where the matrix $M=\left(m_{\mathit{ij}}\right)$ , $i=1,\dots ,N-1,j=1,\dots ,N-1$ :

$$m_{\mathit{ij}}=\left\{\begin{array}{l}0,j\ge i+1,\\ \frac{A_{i-1}\left(2-a_{i-2,i-1}\right)-B_{i-1}{b}_{i-2,i-1}}{A_{i-1}+B_{i-1}},j=i=3,\dots ,N-1,\\ \frac{-A_{i-1}\left(a_{i-j+1,i-1}-2a_{i-j,i-1}+a_{i-j-1,i-1}\right)-B_{i-1}\left(b_{i-j+1,i-1}-b_{i-j-1,i-1}\right)}{A_{i-1}+B_{i-1}},j\le i-1,\end{array}\right)$$ E13

$$m_{1,1}=1,m_{2,2}=\frac{2A_1}{A_1+B_1},m_{i,1}=\frac{B_{i-1}{b}_{i-2,i-1}-A_{i-1}{a}_{i-2,i-1}}{A_{i-1}+B_{i-1}},i=2,\dots ,N-1,$$

$$m_{i,2}=\frac{A_{i-1}\left(2a_{i-2,i-1}-a_{i-3,i-1}\right)+B_{i-1}{b}_{i-3,i-1}}{A_{i-1}+B_{i-1}},i=3,\dots ,N-1.$$

Theorem 1. An explicit finite-difference scheme (9) converges with the first order $\left|{\overline{x}}_k-x_k\right|=O\left(\tau \right)$ if the following condition is satisfied:

$$\tau \le {\tau }_0=\text{min}\left(1,{\left(\frac{2\mathrm{\Gamma }\left(2-{\gamma }_{i-1}\right)}{\lambda \mathrm{\Gamma }\left(3-{\beta }_{i-1}\right)}\right)}^{\frac{1}{{\beta }_{i-1}-{\gamma }_{i-1}}}\right),i=2,\dots ,N-1.$$ E14

Proof. Let ${\overline{X}}_k={\left({\overline{x}}_0,\dots ,{\overline{x}}_{N-2}\right)}^T$ be the exact solution of system (8) and the error vector $e_{k+1}={\overline{X}}_{k+1}-X_{k+1},e_0=0.$ . Then, system (8), with allowance for Lemma 2, can be written as follows:

$$e_{k+1}={\mathit{Me}}_k+F_{e,k}+O\left(\tau \right),$$ E15

where

$$\begin{array}{c}{F}_{e,k}=\frac{1}{A_k+B_k}{\left(\left|f\left(x_1,t_k\right)-f\left({\overline{x}}_1,t_k\right)\right|,\dots ,\left|f\left(x_k,t_k\right)-f\left({\overline{x}}_k,t_k\right)\right|\right)}^T\\ \le \frac{1}{A_k+B_k}\left(L_1{e}_1,\dots ,L_k{e}_k\right)=\mathrm{\Delta }{F}_k{e}_k,\mathrm{\Delta }{F}_k=\frac{1}{A_k+B_k}\mathit{diag}{\left(L_1,\dots ,L_k\right)}^T.\end{array}$$

We note that for any $k$ , the inequality holds $\left|L_k\right|\le L$ . Consider the norm for the matrix $M$ : ${‖M‖}_{\infty }=\underset{i}{\text{max}}\left(\sum _{j=1}^{k-1}{m}_{\mathit{ij}}\right)$ , we obtain

$$\begin{array}{l}{‖M‖}_{\infty }=\underset{1\le i\le N-1}{\text{max}}\\ \left(\begin{array}{l}1+\left|\frac{B_1{b}_{0,1}-A_1{a}_{0,1}}{A_1+B_1}\right|+\left|\frac{2A_1}{A_1+B_1}\right|+\left|\frac{B_2{b}_{1,2}-A_2{a}_{1,2}}{A_2+B_2}\right|\\ +\left|\frac{A_2\left(2a_{1,2}-a_{0,2}\right)+B_2{b}_{0,2}}{A_2+B_2}\right|+\left|\frac{A_2\left(2-a_{1,2}\right)-B_2{b}_{1,2}}{A_2+B_2}\right|\\ +\left|\frac{A_3\left(2a_{1,3}-a_{2,3}-a_{0,3}\right)+B_3\left(b_{0,3}-b_{2,3}\right)}{A_3+B_3}\right|+\dots +\left|\frac{B_{i-1}{b}_{i-2,i-1}-A_{i-1}{a}_{i-2,i-1}}{A_{i-1}+B_{i-1}}\right|\\ +\left|\frac{A_{i-1}\left(2a_{i-2,i-1}-a_{i-3,i-1}\right)+B_{i-1}{b}_{i-3,i-1}}{A_{i-1}+B_{i-1}}\right|\\ +\left|\frac{A_{i-1}\left(2a_{i-3,i-1}-a_{i-2,i-1}-a_{i-4,i-1}\right)+B_{i-1}\left(b_{i-4,i-1}-b_{i-2,i-1}\right)}{A_{i-1}+B_{i-1}}\right|\\ +\left|\frac{A_{i-1}\left(2-a_{i-2,i-1}\right)-B_{i-1}{b}_{i-2,i-1}}{A_{i-1}+B_{i-1}}\right|+\dots \end{array}\right)\end{array}$$ E16

According to Lemma 1, we note that the inequality holds $\frac{2A_1}{A_1+B_1}\ge 0$ . Suppose that the condition is satisfied $A_{i-1}\ge B_{i-1}$ , then the second diagonal element satisfies the inequality $1\le \frac{2A_1}{A_1+B_1}\le 2$ in the matrix $M$ , and the remaining diagonal elements are equal to the inequality $0\le \frac{A_{i-1}\left(2-a_{i-2,i-1}\right)-B_{i-1}{b}_{i-2,i-1}}{A_{i-1}+B_{i-1}}\le 1$ ; these elements with $i\to N-1$ , in view of properties 2 and 3, tend to be zero.

The remaining elements of the matrix (16) also possess these properties. We also note that the matrix $M$ is a matrix with a diagonal predominance for small values $\lambda$ .

Therefore, the sum of the elements of the second row in the matrix $M$ satisfies the condition $1\le \frac{2A_1}{A_1+B_1}+\frac{A_1{a}_{0,1}-B_1{b}_{0,1}}{A_1+B_1}\le 3$ . Further, by virtue of properties 2 and 3 of Lemma 1, it is obvious that the sum of the remaining terms also satisfies these conditions. Therefore, the following estimate is valid for the norm: $1\le {‖M‖}_{\infty }\le 3$ .

Note that for the values of the parameter $\lambda \gg 1$ the norm ${‖M‖}_{\infty }\to 1$ , however, the condition number $\mu \left(M\right)\gg 1$ is violated and the diagonal transformation is violated; therefore, it is necessary to decrease the step $\tau$ .

Further, for any constant $С>0$ independent of $\tau$ , and the error rate, the following estimate holds:

$${‖e_{k+1}‖}_{\infty }\le {‖\mathrm{\Delta }{F}_k+M‖}_{\infty }{‖e_k‖}_{\infty }+\mathit{C\tau }\le \left(3+\frac{L}{A_k+B_k}\right){‖e_k‖}_{\infty }+\mathit{C\tau }.$$ E17

We introduce the notation in Eq. (17): $s_k=\left(3+\frac{L}{A_k+B_k}\right),s=\mathit{C\tau }$ . Then, we obtain the following estimate:

$${‖e_{k+1}‖}_{\infty }\le s_k{‖e_k‖}_{\infty }+s\le s_k\left(s_{k-1}{‖e_{k-1}‖}_{\infty }+s\right)+s=s_k{s}_{k-1}{‖e_{k-1}‖}_{\infty }+s\left(s_k+1\right)$$

$$\le s_k{s}_{k-1}\left(s_{k-2}{‖e_{k-2}‖}_{\infty }+s\right)+s\left(s_k+1\right)=s_k{s}_{k-1}{s}_{k-2}{‖e_{k-2}‖}_{\infty }+s\left(s_k{s}_{k-1}+s_k+1\right)$$

$$\le s_k{s}_{k-1}{s}_{k-2}\left(s_{k-3}{‖e_{k-3}‖}_{\infty }+s\right)+s\left(s_k{s}_{k-1}+s_k+1\right)=s_k{s}_{k-1}{s}_{k-2}{s}_{k-3}{‖e_{k-3}‖}_{\infty }$$ E18

$$\le s_k{s}_{k-1}{s}_{k-2}\left(s_{k-3}{‖e_{k-3}‖}_{\infty }+s\right)+s\left(s_k{s}_{k-1}+s_k+1\right)=s_k{s}_{k-1}{s}_{k-2}{s}_{k-3}{‖e_{k-3}‖}_{\infty }$$

$$\begin{array}{c}+s\left(s_k{s}_{k-1}{s}_{k-2}+s_k{s}_{k-1}+s_k+1\right)\le \dots \le s_k{s}_{k-1}\cdot \dots \cdot s_{k-r}‖e_{k-r}‖\\ +s\left(s_k{s}_{k-1}\cdot \dots \cdot s_{k-r+1}+\dots +s_k+1\right).\end{array}$$

Substituting into Eq. (18) $r=k-1$ , we obtain

$${‖e_{k+1}‖}_{\infty }\le s_k{s}_{k-1}\cdot \dots \cdot s_1‖e_1‖+s\left(s_k{s}_{k-1}\cdot \dots \cdot s_2+\dots +s_k+1\right)\le C_0‖e_0‖+O\left(\tau \right).$$

From the second initial condition (6) it follows: $‖e_1‖\le ‖e_0‖$ and ${С}_0=\prod _{p=1}^k{s}_p$ .

Now according to our assumption $A_{i-1}\ge B_{i-1}$ , which leads us to the relation

$$\tau \le {\left(\frac{2\mathrm{\Gamma }\left(2-{\gamma }_{i-1}\right)}{\lambda \mathrm{\Gamma }\left(3-{\beta }_{i-1}\right)}\right)}^{\frac{1}{{\beta }_{i-1}-{\gamma }_{i-1}}},i=2,\dots ,N-1.$$ E19

Condition (19) begins to work at such values $\lambda$ , when many of conditionalities arise $\mu \left(М\right)\gg 1$ , and for sufficiently small values $\lambda$ , it suffices that the step satisfies the inequality $\tau \le 1$ . Therefore, we arrive at the relation (14). The theorem is proven.

We note that in [38] the authors used the classical Lax theorem, which holds for local finite-difference schemes, to prove the convergence of the scheme. For nonlocal finite-difference schemes, the convergence must be proven independently.

We consider the stability of an explicit finite-difference scheme (4). Suppose that $X_k$ and two $Y_k$ are different solutions of the matrix Eq. (12) with initial conditions $X_0$ and $Y_0$ .

Theorem 2. An explicit finite-difference scheme (9) is conditionally stable if condition (14) is satisfied and the estimate holds $\left|Y_k-X_k\right|\le C\left|Y_0-X_0\right|$ for any $k$ , where $С>0$ does not depend on the step $\tau$ .

Proof. We introduce the notation: $e_{k+1}=Y_{k+1}-X_{k+1}.$ Then, Eq. (12) can be written in the form $e_{k+1}={\mathit{Me}}_k+F_{e,k}$ . Here, as it was said in.

$F_{e,k}=\frac{1}{A_k+B_k}{\left(\left|f\left(x_1,t_k\right)-f\left({\overline{x}}_1,t_k\right)\right|,\dots ,\left|f\left(x_k,t_k\right)-f\left({\overline{x}}_k,t_k\right)\right|\right)}^T\le \frac{1}{A_k+B_k}\left(L_1{e}_1,\dots ,L_k{e}_k\right)=\mathrm{\Delta }{F}_k{e}_k$

According to Theorem 1, we have the following estimate:

$$‖M+\mathrm{\Delta }{F}_k‖\le \left(3+\frac{L}{A_k+B_k}\right)=s_k.$$

Therefore, the following estimate holds

$${‖e_{k+1}‖}_{\infty }\le ‖M+\mathrm{\Delta }{F}_k‖{‖e_k‖}_{\infty }\le \left(3+\frac{L}{A_k+B_k}\right){‖e_k‖}_{\infty }$$

$$=s_k{‖e_k‖}_{\infty }\le s_k{s}_{k-1}{‖e_{k-1}‖}_{\infty }\le s_k{s}_{k-1}{s}_{k-2}{‖e_{k-2}‖}_{\infty }\le \dots \le s_k{s}_{k-1}\cdot \dots \cdot s_{k-r}‖e_{k-r}‖.$$

With $r=k-1$ , we obtain ${‖e_{k+1}‖}_{\infty }\le {С}_0‖e_1‖\le C_0‖e_0‖$ and ${С}_0=\prod _{p=1}^k{s}_p$ .

The last inequality follows the second condition of problem (6). Therefore, if $X_0$ there is a perturbation, then it does not lead to a large increase in the error of the numerical solution. However, for large values $\lambda$ , many of conditionalities $\mu \left(М\right)\gg 1$ arise, and therefore it is necessary to decrease the step $\tau$ ; according to Eq. (19), for small values $\lambda$ , the estimate is valid $\tau \le 1$ . Then, the system is stable if condition (14) is satisfied. The theorem is proven.

4. Results of modeling

Consider the work of the explicit finite-difference scheme (9) on specific examples. We show that the scheme (9) has the first order of accuracy. Since in the general case, the exact solution of the Cauchy problems (5) and (6) cannot be written in analytical form, we will use the double conversion method. For this, we introduce two parameters: $\xi =\underset{i}{\text{max}}\left|x_i-x_{2i}\right|$ —absolute error between the numerical solution $x_i$ in step $\tau$ and the numerical solution $x_{2i}$ in step $\tau /2$ . Then, the order of computational accuracy $p$ can be estimated by the formula

$$p={\text{log}}_2\left(\xi \right)/{\text{log}}_2\left(\tau /2\right).$$

We note that in the case when the fractional parameters in the scheme (9) do not change and have the following values of ${\beta }_k=2$ and ${\gamma }_k=1$ , we arrive at the classical local explicit finite-difference scheme with the second order of accuracy.

The numerical algorithm (9) was implemented in Maple software.

Example 1. Suppose that the right-hand side in Eq. (1) has the form

$$f\left(x\left(t\right),t\right)=\delta \text{sin}\left(\mathit{\phi t}\right)+\mathit{tx}\left(t\right).$$

Then, Eq. (5) describes a linear hereditary Airy oscillator, which was considered in the author’s papers [21, 39] and has the following form.

$${\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)+\lambda {\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)-\mathit{tx}\left(t\right)=\delta \text{cos}\left(\mathit{\phi t}\right).$$

We choose the initial condition (6) for simplicity by homogeneous:

$$x\left(0\right)=\dot{x}\left(0\right)=0.$$

We note that the Airy oscillator is used in optics in the simulation of Airy laser beams.

In this case, the explicit finite-difference scheme (9) has a more specific form:

$$x_0=x_1=0,$$

$$x_{k+1}=\frac{1}{A_k+B_k}\left(\left(2A_k-\mathit{k\tau }\right)x_k-\left(A_k-B_k\right)x_{k-1}\right)$$ E20

$$\begin{array}{c}-\frac{A_k}{A_k+B_k}\sum _{j=1}^{k-1}{a}_{j,k}\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right)\\ -\frac{B_k}{A_k+B_k}\sum _{j=1}^{k-1}{b}_{j,k}\left(x_{k-j+1}-x_{k-j-1}\right)+\delta \text{sin}\left(\mathit{\phi k\tau }\right).\end{array}$$

For the explicit finite-difference scheme (20), we choose the following values of the control parameters: $T=1$ . $\lambda =1$ , $\delta =5$ , $\phi =10$ , $\omega =10$ and $\beta \left(t\right)=1.8-0.03\text{sin}\left(\mathit{\omega t}\right),\gamma \left(t\right)=0.8-0.05\text{cos}\left(\mathit{\omega t}\right)$ . And during the simulation, we will change the number of nodes $N$ in the calculation grid.

Note that the values of the selected parameters for Example 1 satisfy the conditions of Theorems 1 and 2, which is indirectly confirmed by the results of modeling for different values $N$ of the nodes of the computational grid (Table 1).

$N$	$\mathit{\xi }$	$p$
640	0.0003331017	1.119146497
1280	0.0001745618	1.102636795
2560	0.0000906971	1.089811915

Table 1.

Results of numerical simulation.

From Table 1 we can notice that when the number of calculated nodes in the grid doubles in nodes $N$ , the maximum error in absolute value decreases twice, and the order of computational accuracy $p$ tends to unite.

This confirms that the explicit finite-difference scheme (9) and in particular the scheme (20) for Example 1 have the first order of accuracy, and since condition (14) is satisfied, then convergence with the same order.

In Figure 1 the oscillogram (Figure 1a) and the phase trajectory (Figure 1b) are shown for Example 1 at the parameter value $T=10,N=1000.$ It can be noted that with time the amplitude of the oscillations is established and as a result the phase trajectory reaches the limit cycle. Another situation arises in the case of free oscillations $\delta =0$ (Figure 2).

Figure 1.

The oscillogram (a) and the phase trajectory (b) for Example 1 with the parameter values $T=10,N=1000.$

Figure 2.

Oscillogram (a) and phase trajectory (b) for Example 1 with initial conditions $x\left(0\right)=0.1,\dot{x}\left(0\right)=0.2$ , and $\delta =0$ .

The amplitude of the oscillations decays (Figure 2a), and the phase trajectory twists into a spiral (Figure 2b). The dissipation of energy in this case occurs as a result of the presence of friction with a coefficient $\lambda$ and also the “memory” effect, which gives an additional term in the ratio for the total energy of the oscillatory system (Figure 3).

Figure 3.

The oscillogram (a) and the phase trajectory (b) for Example 1 with initial conditions $x\left(0\right)=0.1,\dot{x}\left(0\right)=0.2$ , and $\delta =0$ , $\lambda =0$ .

This fact is confirmed by the results of [2]. Consider the following example of a nonlinear hereditary oscillator.

Example 2. Let that in Eq. (1) the right-hand side has the form.

$$f\left(x\left(t\right),t\right)=\delta \text{sin}\left(\mathit{\phi t}\right)-\mathit{ax}\left(t\right)+{\mathit{bx}}^3\left(t\right),$$

and we choose the initial conditions (6) to be homogeneous:

$$x\left(0\right)=\dot{x}\left(0\right)=0.$$

In this case, Eq. (5) describes the Duffing fractional oscillator [18]:

$${\partial }_{0t}^{\beta \left(t\right)}x\left(\eta \right)+\lambda {\partial }_{0t}^{\gamma \left(t\right)}x\left(\eta \right)+{\mathit{bx}}^3\left(t\right)-\mathit{ax}\left(t\right)=\delta \text{sin}\left(\mathit{\phi t}\right).$$

The explicit finite-difference scheme (9) for this case has the form

$$x_0=x_1=0,$$

$$x_{k+1}=\frac{1}{A_k+B_k}\left(\left(2A_k+1\right)x_k-x_k^3-\left(A_k-B_k\right)x_{k-1}\right)$$ E21

$$\begin{array}{c}-\frac{A_k}{A_k+B_k}\sum _{j=1}^{k-1}{a}_{j,k}\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right)\\ -\frac{B_k}{A_k+B_k}\sum _{j=1}^{k-1}{b}_{j,k}\left(x_{k-j+1}-x_{k-j-1}\right)+\delta \text{sin}\left(\mathit{\phi k\tau }\right).\end{array}$$

For the explicit finite-difference scheme (21), we take the values of the control parameters as follows: $T=1$ , $\lambda =0.3$ , $\delta =2$ , and $\phi =\omega =1$ .

Remark. Note that this choice of control parameter values is ensured by the condition (14) for Theorems 1 and 2. The results of numerical simulation for Example 2 are given in Table 2.

$N$	$\mathit{\xi }$	$p$
640	0.0003619281	1.107545912
1280	0.0001896841	1.092050182
2560	0.0000991471	1.079382204

Table 2.

Results of numerical simulation.

Note from Table 2 that for Example 2, with an increase in the number of design nodes $N$ , the maximum error $\xi$ in absolute value decreases and the order of computational accuracy $p$ tends to unite. This indicates that the explicit finite-difference scheme (21) has the first order of accuracy.

Let’s perform numerical simulation according to the scheme (21) with the values of the following parameters, $T=100$ , $N=2000$ , and $\delta =50$ , and leave the remaining parameters unchanged. Let us construct an oscillogram and a phase trajectory (Figure 3).

The oscillogram (Figure 4a) has a constant amplitude of a more complex shape at its minima and maxima, which is reflected in the phase trajectory (Figure 4b). The phase trajectory enters a complex two-loop limit cycle. The presence of such loops, apparently, is associated with the effects of memory in the oscillatory system [40].

Figure 4.

The oscillogram (a) and the phase trajectory (b) for Example 2.

Figure 5 shows the case of free oscillations for Example 2. It is seen that the presence of friction and memory effects in the oscillatory system intensify energy dissipation, which leads to damping of the oscillations (Figure 5a) and a phase trajectory—a twisting spiral (Figure 5b). Indeed, if there is no friction $\lambda =0$ in the oscillatory system, we obtain an oscillogram and a phase trajectory as in Figure 6.

Figure 5.

Oscillogram (a) and phase trajectory (b) for Example 3 with initial conditions $x\left(0\right)=0.1,\dot{x}\left(0\right)=0.2$ , and $\delta =0$ .

Figure 6.

Oscillogram (a) and phase trajectory (b) for Example 2 with initial conditions $x\left(0\right)=0.1,\dot{x}\left(0\right)=0.2$ , and $\delta =0$ , $\lambda =0$ .

Example 3. Suppose that in Eq. (1) the right-hand side has the form

$$f\left(x\left(t\right),t\right)=\mathit{bt}+c\sum _{n=1}^7{a}_n\text{sin}\left(\mathit{nx}\left(t\right)\right)-{\omega }^{\beta \left(t\right)}x\left(t\right),$$ E22

where $b$ is the spring travel speed; $c$ is the surface adhesion energy; $\omega$ is the frequency of free oscillations; and $a_n=2n\underset{0}{\overset{1}{\int }}\frac{\text{cos}\left(\mathit{\pi n\tau }\right)\mathit{d\tau }}{{\text{cosh}}^2\left(\mathit{\pi \tau }\right)}$ is the coefficients of the expansion of the Fourier series.

Eq. (5) with the right-hand side of Eq. (22) describes the hereditary stick–slip effect [20]. The stick–slip effect is encountered in tribology problems, for example, when the movement of a load on a spring along a surface is investigated. Due to adhesion, the load adheres to the surface, and due to the tension of the spring, it breaks and slides along it, and its oscillations occur [41, 42]. The stick–slip effect can also be incorporated into the mechanical model of an earthquake in the subduction zone of lithospheres plates [43].

In [43] it was said that in order to obtain a reliable solution it suffices to take the first seven coefficients $a_n$ in the expansion of the function (22). The values of these coefficients are taken from [43] $a_1=0.436$ , $a_2=0.344$ , $a_3=0.164$ , $a_4=0.058$ , $a_5=0.021$ , $a_6=0.004$ , and $a_7=0.003$ . Values of control parameters are $\beta \left(t\right)=1.8-0.03\text{sin}\left(\mathit{\pi t}\right)$ $\gamma \left(t\right)=0.6-0.04\text{cos}\left(\mathit{\pi t}\right)$ , $N=3000$ $\delta =50$ , $\tau =0.05$ , $\lambda =0.3$ , $b=1$ , $\omega =1$ , and $x\left(0\right)=0,\dot{x}\left(0\right)=0.3$ .

Figure 7 shows the calculated displacement curves, displacement velocities, and phase trajectory. Figure 7a shows the oscillogram for Example 3. It can be seen that during the separation, the load experiences oscillations and the rate of such oscillations in the potential well attenuates rather slowly (Figure 7b). This effect is the eradication of the process. The phase trajectory in Figure 7c shows that the potential wells are stable focuses.

Figure 7.

Calculated curves obtained from formula (6) from [20] (curve 1) and formula (9) (curve 2): (a) oscillogram, (b) oscillator speed, and (c) phase trajectory.

5. Conclusion

A mathematical model characterizing a wide class of hereditary oscillators is proposed and studied. The model is a differential Cauchy problem with derivatives of fractional-order variables of the Gerasimov-Caputo types (5) and (6). Using the theory of finite-difference schemes, a nonlocal explicit finite-difference scheme (9) was constructed with the first order of accuracy. Questions of its stability and convergence, which are formulated in the form of corresponding theorems, were studied.

The main result of the paper can be formulated as follows: an explicit finite-difference scheme is conditionally stable and converges if criterion (14) is satisfied. With the help of computational examples, it was shown that the scheme (9) has the first order of accuracy. It is confirmed that in the case of free oscillations, the presence of friction and heredity increases dissipation of energy, which leads to attenuation of oscillations.

One of the continuations of the investigation of the Cauchy problems (5) and (6) is a generalization of it:

$${\partial }_{0t}^{\beta \left(x\left(t\right),t\right)}x\left(\eta \right)+\lambda \left(x\left(t\right),t\right){\partial }_{0t}^{\gamma \left(x\left(t\right),t\right)}x\left(\eta \right)=f\left(x\left(t\right),t\right),x\left(0\right)={\alpha }_0,\dot{x}\left(0\right)={\alpha }_1.$$

Another continuation of the research is related to the introduction of other memory functions $K_1\left(t-\tau \right),K_2\left(t-\tau \right)$ into the model Eq. (1), which leads to different model equations with different derivatives of fractional orders, and also the Cauchy problems (5) and (6) can be written in terms of the local fractional derivative [44–46].

The question of the stability of the rest points of dynamical systems described by the Cauchy problems (5) and (6) is also interesting, by analogy with the papers [47, 48].

Acknowledgments

The work was carried out according to the state task within the framework of research work Vitus Bering Kamchatka State University on the topic “Application of fractional calculus in the theory of oscillatory processes” No.AAAA-A17-117031050058-9 and with the support of the grant of the President of the Russian Federation MK-1152.2018.1.