Mathematical Model of Fractional Dufﬁng Oscillator with Variable Memory

: The article investigates a mathematical model of the Dufﬁng oscillator with a variable fractional order derivative of the Riemann–Liouville type. The study of the model is carried out using a numerical scheme based on the approximation of the fractional derivative of the Riemann–Liouville type by a discrete analog—the fractional derivative of Grunwald–Letnikov. The adequacy of the numerical scheme is veriﬁed using speciﬁc examples. Using a numerical algorithm, oscillograms and phase trajectories are constructed depending on the values of the model parameters. Chaotic regimes of the Dufﬁng fractional oscillator are investigated using the Wolf–Bennetin algorithm. The forced oscillations of the Dufﬁng fractional oscillator are investigated using the harmonic balance method. Analytical formulas for the amplitude-frequency, phase-frequency characteristics, and also the quality factor are obtained. It is shown that the fractional Dufﬁng oscillator possesses different modes: regular, chaotic, multi-periodic. The relationship between the order of the fractional derivative and the quality factor of the oscillatory system is established. We group expressions by integrals from sine cosine differentials cosine


Introduction
At present, one of the scientific directions has received wide development-hereditary dynamics, the founder of which is the Italian mathematician Vito Volterra [1]. Hereditarity is a property of a dynamic system in which its current state depends on previous states. Hereditary dynamical systems find their application in hereditary mechanics [2], biology [3], economics [4] and other fields of knowledge.
From the point of view of mathematics, a hereditary dynamical system can be described using integro-differential equations with difference kernels, which are called memory functions [1]. If the memory functions are chosen as exponential, then we turn to equations with derivatives of fractional derivatives, which are studied in the framework of fractional calculus [5][6][7]. Fractional derivatives have many definitions and unique properties, but all of them to some extent describe the memory effect that characterizes information about previous states of the system. This effect determines additional degrees of freedom-orders of fractional derivatives. The orders of fractional derivatives can be constant or functions of time, which in the latter case significantly complicates the study of a dynamical system.
It should be noted that in article [8] the Duffing oscillator with a constant order Riemann-Liouville derivative was considered. Chaotic and regular regimes are investigated using a numerical scheme and various tests. Similar studies were carried out in article [9], but with the Caputo derivative of fractional constant order. In article [10], on the basis of Melnikov's method, the necessary condition for a chaotic regime of a Duffing oscillator with a fractional order derivative under harmonic excitation is investigated. In article [11], forced oscillations of the Duffing oscillator with a constant-order Caputo derivative were investigated. The amplitude-frequency, phase-frequency characteristics, as well as the bistable mode are investigated.
Article [12] proposes an analytical method for solving the Duffing equation with fractional Caputo derivatives based on a combination of the perturbation and Laplace methods. Paper [13] uses He's fractional derivative to describe the fractional Duffing equation, which is determined using a variational iterative algorithm. To solve such an equation, a complex fractional transformation was used. In article [14], a solution of the Duffing oscillator with fractional derivatives was obtained using the method of homotopy analysis.
This article proposes a generalization of the Duffing oscillator, a nonlinear oscillatory system with chaotic dynamics. This generalization is based on the introduction of a variable order derivative of the Riemann-Liouville type to describe the energy dissipation of an oscillatory system. This article analyzes the chaotic and regular modes, frequency characteristics of the Duffing oscillator with variable memory, which is a generalization of the research from article [8].
The work is organized as follows. Section 2 gives some definitions from the theory of fractional disappearance; Section 3 presents a statement of the Cauchy problem for the Duffing equation with a fractional variable-order derivative of the Riemann-Liouville type in the dissipative term. It also shows a numerical scheme obtained on the basis of the discrete fractional Grunwald-Letnikov derivative. It is shown on a specific example that the circuit has the first order of accuracy. Section 4 contains an algorithm for calculating the maximum Lyapunov exponents. It is shown with a specific example that chaotic regimes are possible under certain conditions. In Section 5, we studied forced oscillations of the Duffing fractional oscillator. Using the harmonic balance method, formulas are derived for calculating the amplitude-frequency and phase-frequency characteristics, as well as the quality factor. Surfaces of the corresponding characteristics are constructed. The conclusion is given in Section 6.

Preliminary Material on Fractional Calculus
Definition 1. The fractional derivative of variable order 0 < q (t) < 1 of the Riemann-Liouville type has the form: where Γ (1 − q (t)) is Euler's gamma function.

Definition 2.
The fractional derivative of variable order 0 < q (t) < 1 of the Grunwald-Letnikov type has the form [15]: where q (t) j -binomial coefficient.
The weighting factor in Formula (2) can be determined as follows:

Problem Statement and Numerical Solution
Consider the following Cauchy problem with offset [16]: where x (t)-displacement, solution function, x 0 and y 0 -given constants, λ is the coefficient of friction, δ and ω are the amplitude and frequency of the external periodic impact, 0 < q (t) < 1, t ∈ [0, T], T > 0-simulation time.

Remark 2.
The variable order q (t) of the fractional derivative determines the intensity of energy dissipation in the oscillatory system. If this order is constant and equal to one, then the Cauchy problem (3) becomes the Cauchy problem for the classical Duffing oscillator [4] .
Let x (t) ∈ C 3 [0, T] to achieve the required smoothness. Divide the time interval [0, T] into N equal parts with a constant step h = T/N. The solution function x (t) will go to the grid function x (t k ) = x k , where t k+1 = (k + 1) h, k = 1, . . . , N − 1.

Remark 3. The following relation is valid
The Cauchy problem (3) may be represented as a system: Let us rewrite the Cauchy problem (5) in the difference setting taking into account the relation (4) Remark 4. Note that relation (6) is an explicit Euler approximation, which is conditionally stable, and taking into account relation (4), the scheme has the first order of accuracy.
Let us show the last remark with a specific example.
Example 1. Consider the case: Consider the change of the absolute error ε and computational accuracy α = ln (ε) / ln (h/2) of scheme (6) when changing the step h. To calculate the absolute error, we use the Runge rule: ε = max

Lyapunov Exponents
In the study of nonlinear systems, one of the important tasks is to determine the type of oscillations-periodic, quasi-periodic, random, chaotic. It is especially difficult to distinguish quasiperiodic oscillations from chaotic and random ones, since quasiperiodic oscillations often have a very complex shape that is visually indistinguishable from "random" [17].
A feature of chaotic oscillations is their high sensitivity to small changes in initial conditions. Therefore, one of the most reliable ways to detect chaos is to determine the speed of trajectory run-up, which is estimated using a spectrum of Lyapunov exponents [18].
The geometric meaning of the Lyapunov exponents is that two solutions whose initial values are located in a certain neighborhood of radius will diverge in time T into an n-dimensional ellipsoid along n main semi-axes, and at time t the radii will be determined by the values where Λ i -Lyapunov exponents.
If two trajectories remain close over time, then their change obeys the law where J is Jacobi matrix, u j = (x j , y j ) T is column vector, f 1 , f 2 -right sides of the system (5). The Formula (7) defines the perturbed equations.

Remark 5.
The presence of at least one positive Lyapunov exponent in the spectrum means the presence of a chaotic regime (asymptotic instability) of the phase trajectory under consideration. A negative reading indicates a limit cycle, while a zero reading indicates that analysis is not enough to draw a conclusion about stability and instability according to Lyapunov. The vectors of the initial conditions for Equation (7) must be orthogonal.
Consider the Benettin-Wolf algorithm for calculating Lyapunov exponents and implement it in the MAPLE environment [19].
1. Select the starting point of the vector x 0 and together with it we will track the K disturbed trajectories.
Let K = 3. 2. We solve the original equation numerically together with three sets of perturbed equations or equations in variations (7). As the initial vectors for equations in variations, you must select a set of vectorsx 0 0 ,ỹ 0 0 ,z 0 0 , that are orthogonal and normalized by one. 3. After time T, the trajectory will move to the point vector x 1 , the perturbation vectorsx 1 ,ỹ 1 ,z 1 , which we renormalize using the Gram-Schmidt method.
x 0 1 =x 1 x 1 , 1. Next, we continue counting from the point x 1 and the perturbation vectorsx 0 1 ,ỹ 0 1 ,z 0 1 . After a regular time interval T, we get a new set of perturbation vectorsx 2 ,ỹ 2 ,z 2 , which is subjected to orthogonalization and renormalization.
2. Steps 2-4 are repeated M times, and the sums of the calculations are calculated: , which include perturbation vectors before renormalization, but after orthogonalization. 3. The estimation of the first three Lyapunov exponents is calculated using the formula: For a more simplified calculation, take K = 2. As a result, we get the following equations: where xx k , yy k , XX k , YY k solutions of perturbed Equations (9) and (10) (Figure 2b-d). The figures show that for positive values of Lyapunov exponents, chaotic modes are observed, and the lower the value λ, the more closely the trajectory approaches that of the classical Duffing oscillator. For negative values of Lyapunov exponents, the phase trajectories are limit cycles, and the greater the modulus of the negative Lyapunov exponent, the more stable the limit cycle will be.

Forced Oscillations
As shown by the results of articles [20,21], the orders of fractional derivatives in oscillatory systems with memory are associated with the quality factor and are responsible for the intensity of energy dissipation. Therefore, let us investigate the forced oscillations of the Duffing fractional oscillator (3).
Consider the dependence of the amplitude A and phase ϕ of steady-state oscillations for operator (1) on the frequency of the external force. To calculate the frequency response and frequency response, we use the harmonic linearization method [21][22][23].
Consider the Cauchy problem (3). We look for a solution in the form: We group expressions by integrals from sine and cosine and differentials from integrals from sine and cosine [21].
[cos(ϕ) sin(ωt) + sin(ϕ) cos(ωt)] t 0 v −q(t) cos(ωv)dv [− cos(ϕ) cos(ωt) + sin(ϕ) sin(ωt)] t 0 v −q(t) sin(ωv)dv Thus, the Riemann-Liouville operator for solution (11) has the form: In the case of steady oscillations at large times t (t → ∞), the integrals (12) can be written in the form [23]: Substituting (13) and (14) in (12) we'll get Substituting (11) and (15) in Equation (3), given that: x(t) = A cos(ωt + ϕ), x (t) = −Aω sin(ωt + ϕ), x (t) = −Aω 2 cos(ωt + ϕ), We get the equation [23]: where There In practice it is enough to take n = 50. Equation (16) is a classical linear oscillator for which the formulas for the Amplitude-frequency and Phase-frequency characteristics (AFC and PFC) are derived using the harmonic linearization method: The Q-factor for a fractional Duffing oscillator can be determined from Equation (16) using the formula [21,23]: where s and p are defined by the Equalities (17) and (18).  Figure 4 shows that when the q parameter is reduced, the quality-factor increases (21). The maximum amplitude corresponds to the maximum Q-factor, and when the frequency decreases, the Q-factor decreases. However, most of all, Q-factor depends on the parameter q(t). Figure 5 shows AFC and PFC versus amplitude and phase versus frequency and the exponent of the derivative q(t). In Figure 5a, we see separate areas from low frequencies to high frequencies.

Conclusions
The article proposes a mathematical model of a fractional Duffing oscillator with a variable fractional order derivative of the Riemann-Liouville type in the dissipative term. Using the discrete Grunwald-Letnikov derivative, a first-order numerical scheme is constructed. Using Runge's method, it is shown that with an increase in the number of nodes in the computational grid, the order of computational accuracy tends to unity. A numerical algorithm is proposed for calculating the spectra of the maximum Lyapunov exponents. It is shown that the spectrum of the maximum Lyapunov exponents can contain positive values, which indicates the existence of chaotic regimes. Knowing the conditions for the existence of chaotic regimes, we can investigate regular regimes, for example, regimes of steady oscillations under the action of an external harmonic force. The study of forced vibrations using analytical formulas for the frequency characteristics and quality factor of the oscillatory system showed that the variable order of the fractional derivative determines the intensity of energy dissipation in the system. With an increase in the order, the quality factor of the oscillatory system increases, and with a decrease, it increases. This fact suggests that oscillatory modes can be controlled using the order of the fractional derivative. This fact is very important for the development of applied engineering research. For example, in machine learning and intelligent communications [24].
It is necessary to study the fractional Duffing oscillator with a variable fractional derivative of the Riemann-Liouville type, taking into account different definitions and types of its function [15]. It is of interest to study forced oscillations, resonance curves as in the Figure 5a. Furthermore, also the presence of bistable behavior [11].
It is of interest to study the fractional Duffing oscillator with higher order nonlinearities, for example, by analogy with the article [25].  Acknowledgments: I would like to express my gratitude to Arsen Pskhu for valuable advice, which served to better understand the research results.

Conflicts of Interest:
The authors declare no conflict of interest.