Mathematical Models of Oscillators with Memory

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.


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 [6-9, 11-14, 18-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 Wrighttype 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 selfoscillator 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.

Formulation of the problem
Consider the following model integrodifferential equation for the function x t ð Þ ∈ C 3 0; T ð Þ, where T > 0: where € x t ð Þ ¼ d 2 x=dt 2 , _ x t ð Þ ¼ dx=dt, λ > 0-given constant. Eq. (1) describes a wide class of hereditary, depending on the form of the righthand side (function f x t ð Þ; t ð Þ) of linear or nonlinear oscillators.
(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].
Þon the right-hand side of Eq. (1) satisfies a Lipschitz condition with respect to a variable x t ð Þ: L-Lipschitz constant. Eq. (1) describes a broad class of hereditary nonlinear oscillators, depending on the form of the function f x t ð Þ; t ð Þon its right-hand side, and the parameter λ has the meaning of the coefficient of friction.
Note that in Definition 1, the memory functions K 1 t À η ð Þand K 2 t À η ð Þ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 where γ t ð Þ, β t ð Þ ∈ C 0; T ½ , Γ t ð Þ-Euler gamma function. We give the following definitions. Definition 3. Derivatives of fractional variables of orders β t ð Þ and γ t ð Þ Gerasimov-Caputo type: we call the following operators of fractional differentiation: We note that in the case when the functions β t ð Þ and γ t ð Þ 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 β ¼ 2 and γ ¼ 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: For Eq. (5), the initial conditions in the local formulation are valid: where α 0 and α 1 -given constants. As a result, we arrive at the Cauchy problems (5) and (6), which we will investigate.

Explicit finite-difference scheme
We construct an explicit finite-difference scheme for the Cauchy problems (5) and (6). We divide the time interval 0; T ½ into N equal parts with a step τ ¼ T=N. We introduce the grid function x t k ð Þ ¼ x k , where t k ¼ kτ, k ¼ 1, …, N À 1. The derivatives of the variables of fractional orders in Eq. (5) are approximated according to the relations in [36,37]: then the formulas will refer to the formula (7) Here, to shorten the record, β t k ð Þ ¼ β k , γ t k ð Þ ¼ γ k . Taking into account relation (7), the Cauchy problems (5) and (6) in the difference formulation will have the form Here, to shorten the record, f k ¼ f x k ; t k ð Þ. We write the problem (8) explicitly: 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: Proof. The first property follows the definition of weight coefficients: The second property is proven in the following way. We introduce two functions: where x > 0. These functions are decreasing. Really derived from these functions Therefore, the second property holds. The third property follows also the properties of the gamma function. The lemma is proven.
Then, we have the following lemma.
x η ð Þ satisfy the following estimates: where С 1 and С 2 -constants that are independent of the parameter τ.
Proof. Using the first property of Lemma 1 and Definition 3, we obtain Similarly, we can obtain the second estimate from Eq. (10).
The lemma is proven.
Investigation. According to Lemma 2, it can be shown that the explicit finitedifference scheme (9) has an error ε ¼ O τ ð Þ. 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: 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 Theorem 1. An explicit finite-difference scheme (9) converges with the first order x k À x k j j¼ O τ ð Þ if the following condition is satisfied: Proof. Let X k ¼ x 0 ; …; x NÀ2 ð Þ T be the exact solution of system (8) and the error vector e kþ1 ¼ X kþ1 À X kþ1 , e 0 ¼ 0:. Then, system (8), with allowance for Lemma 2, can be written as follows: where We note that for any k, the inequality holds L k j j≤ L. Consider the norm for the According to Lemma 1, we note that the inequality holds 2A 1 A 1 þB 1 ≥ 0. Suppose that the condition is satisfied A iÀ1 ≥ B iÀ1 , then the second diagonal element satisfies the inequality 1 ≤ 2A 1 A 1 þB 1 ≤ 2 in the matrix M, and the remaining diagonal elements are equal to the inequality 0 ≤ A iÀ1 þB iÀ1 ≤ 1; these elements with i ! 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 λ.
Therefore, the sum of the elements of the second row in the matrix M satisfies the condition 1 ≤ 2A 1 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 ≤ M k k ∞ ≤ 3. Note that for the values of the parameter λ ≫ 1 the norm M k k ∞ ! 1, however, the condition number μ M ð Þ ≫ 1 is violated and the diagonal transformation is violated; therefore, it is necessary to decrease the step τ.
Further, for any constant С > 0 independent of τ, and the error rate, the following estimate holds: We introduce the notation in Eq. (17): Then, we obtain the following estimate: ≤ s k s kÀ1 s kÀ2 s kÀ3 e kÀ3 k k ∞ þ s À Á þ s s k s kÀ1 þ s k þ 1 ð Þ ¼ s k s kÀ1 s kÀ2 s kÀ3 e kÀ3 k k ∞ þ s s k s kÀ1 s kÀ2 þ s k s kÀ1 þ s k þ 1 ð Þ ≤ … ≤ s k s kÀ1 Á … Á s kÀr e kÀr k k þ s s k s kÀ1 Á … Á s kÀrþ1 þ … þ s k þ 1 ð Þ : From the second initial condition (6) it follows: e 1 k k≤ e 0 k k and С 0 ¼ Q k p¼1 s p . Now according to our assumption A iÀ1 ≥ B iÀ1 , which leads us to the relation Condition (19) begins to work at such values λ, when many of conditionalities arise μ М ð Þ ≫ 1, and for sufficiently small values λ, it suffices that the step satisfies the inequality τ ≤ 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 Y k À X k j j≤ C Y 0 À X 0 j jfor any k, where С > 0 does not depend on the step τ.
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 ¼ Me k þ F e, k . Here, as it was said in. F e, k ¼ 1 Þ¼ΔF k e k According to Theorem 1, we have the following estimate: Therefore, the following estimate holds k k ∞ ≤ s k s kÀ1 e kÀ1 k k ∞ ≤ s k s kÀ1 s kÀ2 e kÀ2 k k ∞ ≤ … ≤ s k s kÀ1 Á … Á s kÀr e kÀr k k: 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 λ, many of conditionalities μ М ð Þ≫ 1 arise, and therefore it is necessary to decrease the step τ; according to Eq. (19), for small values λ, the estimate is valid τ ≤ 1. Then, the system is stable if condition (14) is satisfied. The theorem is proven.

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: ξ ¼ max i x i À x 2i j j -absolute error between the numerical solution x i in step τ and the numerical solution x 2i in step τ=2. Then, the order of computational accuracy p can be estimated by the formula We note that in the case when the fractional parameters in the scheme (9) do not change and have the following values of β k ¼ 2 and γ 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 Þ¼δ sin φt ð Þ þ tx t ð Þ: Then, Eq. (5) describes a linear hereditary Airy oscillator, which was considered in the author's papers [21,39] and has the following form.
We choose the initial condition (6) for simplicity by homogeneous: 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: For the explicit finite-difference scheme (20), we choose the following values of the control parameters: T ¼ 1. λ ¼ 1, δ ¼ 5, φ ¼ 10, ω ¼ 10 and β t ð Þ ¼ 1:8 À 0:03 sin ωt ð Þ, γ t ð Þ ¼ 0:8 À 0:05 cos ωt ð Þ. 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).
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 δ ¼ 0 ( Figure 2).
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 λ and also the "memory" effect, which gives an additional term in the ratio for the total energy of the oscillatory system ( Figure 3). This fact is confirmed by the results of [2]. Consider the following example of a nonlinear hereditary oscillator.
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.  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 ξ 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 δ ¼ 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 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 λ ¼ 0 in the oscillatory system, we obtain an oscillogram and a phase trajectory as in Figure 6.
Example 3. Suppose that in Eq. (1) the right-hand side has the form f x t ð Þ; t ð Þ¼bt þ c ∑ 7 n¼1 a n sin nx t ð Þ ð ÞÀω β t ð Þ x t ð Þ, where b is the spring travel speed; c is the surface adhesion energy; ω is the frequency of free oscillations; and a n ¼ 2n Ð 1 0 cos πnτ ð Þdτ cosh 2 πτ ð Þ 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].

Conclusion
A mathematical model characterizing a wide class of hereditary oscillators is proposed and studied. The model is a differential Cauchy problem with derivatives