Dubovsky’s Class of Mathematical Models for Describing Economic Cycles with Heredity Effects

Abstract

The article is devoted to the study of economic cycles and crises, which are studied within the framework of the theory of N.D. Kondratiev long waves (K-waves). The object of the study is the fractional mathematical models of S. V. Dubovsky, consisting of two nonlinear differential equations of fractional order and describing the dynamics of the efficiency of new technologies and the efficiency of capital productivity, taking into account constant and variable heredity. Fractional mathematical models also take into account the dependence of the rate of accumulation on capital productivity and the influx of external investment and new technological solutions. The effects of heredity lead to a delayed effect of the response of the system in question to the impact. The property of heredity in mathematical models is taken into account using fractional derivatives of constant and variable orders in the sense of Gerasimov–Caputo. The fractional mathematical models of S. V. Dubovsky are further studied numerically using the Adams–Bashforth–Moulton algorithm. Using a numerical algorithm, oscillograms and phase trajectories were constructed for various values and model parameters. It is shown that the fractional mathematical models of S. V. Dubovsky may have limit cycles, which are not always stable.

1. Introduction

Numerous research results show that N.D. Kondratiev long waves (K-waves) play an important role in forecasting economic crises and cycles [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. K-waves can be described using dynamic (oscillatory) systems, which mainly characterize nonlinear oscillatory processes [11,14,18,27]. Such processes are studied to describe the chaotic and regular dynamics of nonlinear oscillatory systems. In practice, identifying chaotic regimes, for example, in clinical medicine or in the study of turbulence in the atmosphere, is one of the main tasks. Regular modes of nonlinear oscillations are of particular interest when studying economic cycles and crises, various chemical reactions, biological systems, etc. [35,36]. Nonlinear oscillatory processes within the framework of the K-wave theory can also have heredity (memory) effects, which can be described from a mathematical point of view using derivatives of fractional orders [37,38,39,40,41].

The review work of Tarasov V.E. [42] provides 260 literature sources on the application of fractional derivatives to various economic processes. This article contains a brief overview of the stages in the history of the application of fractional calculus in modern mathematical economics and economic theory. The first stage of the application of heredity in economics is associated with the works published in 1966 and 1980 by Clive W. J. Granger [43,44], who received the Nobel Prize in Economics in 2003. Further, the history of the application of fractional calculus in economics can be divided into the following five stages of development (approaches): ARFIMA (models FARMA, GARMA, etc.) [43,44,45,46,47]; fractional Brownian motion [48]; econophysics [49,50]; deterministic chaos [51]; mathematical economics [52,53]. The modern stage (mathematical economics) is intended to include new economic concepts and concepts into modern economic theory that allow us to take into account the presence of memory in economic processes. Further, in the review article, some comments are given on possible further directions for the development of fractional mathematical economics.

In this article, within the framework of mathematical economics, we will study K-waves using a generalization of the classical mathematical model of S.V. Dubovsky, proposed in article [14] for the case of memory effects. Note that K-waves are well described in monographs [15,16,17,20,21,22,24]. However, an analysis of the literature on the research topic showed that the application of fractional calculus to the dynamic system of S.V. Dubovsky has not yet been carried out to describe K-waves. In addition, even the classical dynamical system of S.V. Dubovsky without taking into account heredity (with integer derivatives), proposed in article [14], did not take into account the influx of external investments and new technologies, as well as the rate of accumulation as a function of capital productivity. Therefore, the proposed mathematical models with derivatives of fractional constant and variable orders in the sense of Gerasimov–Caputo are completely new objects of study in this article.

The main objective of the study is to establish the existence of closed-phase trajectories in the proposed fractional mathematical models. Next, in the presence of closed trajectories using computer modeling, we study how the values of various parameters of the model and their various functional dependencies influence the shape of these cycles, as well as study cycles for stability. Other dynamic regimes arise that also deserve research.

The research plan in the article is structured as follows: The Section 1 provides an overview of the current state of research on the theory of K-waves, revealing some problems that have not been studied within the framework of the classical model of S.V. Dubovsky and which are solved in this article. The Section 2 examines the classical model of S.V. Dubovsky: the basic equations of the model are derived, and two modifications of the model are given, taking into account the rate of accumulation depending on capital productivity and external inflow of investments and new technologies. It is shown using Bendixson’s theorem that the proposed model can have a closed-phase trajectory. Oscillograms and phase trajectories were constructed using the Adams–Bashforth–Moulton numerical algorithm and an interpretation of the research results was given. The Section 3 provides a generalization of the classical model of S.V. Dubovsky in the case of taking into account constant heredity: the definition of the fractional derivative in the sense of Gerasimov–Caputo is given; the Cauchy problem is posed, which is studied by the Adams–Bashforth–Moulton numerical method; the order of convergence of the method is studied using the Runge rule; oscillograms and phase trajectories are constructed for various values of the model parameters and an interpretation of the research results is given. The Section 4 provides a generalization of the fractional model of S.V. Dubovsky in the case of taking into account variable heredity: the definition of a derivative of a fractional variable order of the Gerasimov–Caputo type is introduced; the Cauchy problem is posed, which is studied using the Adams–Bashforth–Moulton numerical method; the order of convergence of the method is studied using the Runge rule; oscillograms and phase trajectories are constructed for various values of the model parameters; an interpretation of the research results is given. The Section 5 provides conclusions based on the results of the study and also indicates their possible further development.

2. Classical Mathematical Model of S.V. Dubovsky and Some Modifications for Describing K-Waves

In this section, the classical mathematical model of S.V. Dubovsky was studied to describe N.D. Kondratiev long waves (K-waves). This model describes the dynamics of free fluctuations in the efficiency of new technologies and the efficiency of capital productivity. From the point of view of mathematics, it is a system of nonlinear ordinary differential equations of the first order (Cauchy problem).

The purpose of the research is to visualize the results of the solution using numerical modeling, taking into account the modification of the mathematical model of S.V. Dubovsky, which consists of taking into account the dependence of the accumulation rate on capital productivity and external inflow of investments and new technological solutions.

2.1. Derivation of the Basic Equations of the Classical Mathematical Model by S.V. Dubovsky

Derivation of the basic equations of the classical mathematical model by S.V. Dubovsky

(1)${\displaystyle \frac{\dot{Y}\left(t\right)}{Y\left(t\right)}}=n{\displaystyle \frac{\dot{K}\left(t\right)}{K\left(t\right)}}+\left(1-n\right)\left(l+{\displaystyle \frac{\dot{U}\left(t\right)}{U\left(t\right)}}\right),$

(2)${\displaystyle \frac{\dot{K}\left(t\right)}{K\left(t\right)}}={\displaystyle \frac{nY\left(t\right)}{K\left(t\right)}}-\mu ,$

(3)${\displaystyle \frac{\dot{U}\left(t\right)}{U\left(t\right)}}={\displaystyle \frac{nY\left(t\right)}{K\left(t\right)}}\left({\displaystyle \frac{u\left(t\right)}{U\left(t\right)}}-1\right),$

where $Y\left(t\right)$ is the gross domestic product, $K\left(t\right)$ is capital, $U\left(t\right)$ is the average technological level of the economy, n is the accumulation rate, $\mu$ is the capital depreciation rate, l is the employment growth rate and $u\left(t\right)$ is the latest technological level for newly introduced fixed production assets.

On the right side of Equation (1), in order to simplify the model, all types of income from production and/or other capital (working capital, banking, etc.) not used for investment are omitted. Equation (1) includes an indicator of quantitative measurement of the results of scientific and technological progress according to Harrod, which increases labor productivity.

Differential Equation (2) is more interesting from the point of view of revealing the fundamentals of the production function and is the standard equation for the dynamics of the overall production function.

A more detailed derivation of Equation (3) to describe the dynamics of the average technological level (the standard equation used in economics) is presented in a number of scientific works by S.V. Dubovsky [14,25].

S.V. Dubovsky, based on the model of N.D. Kondratiev and standard equations reflecting economic dynamics, in particular, equations for the rate of employment growth and the rate of depreciation of production assets, the latter depending on the rate of renewal of the latter, expressed the general interaction and interdependence in the form of the following formulas:

(4)$\mu \left(t\right)={\mu }_0{\displaystyle \frac{K_0}{I_0}}{\displaystyle \frac{I\left(t\right)}{K\left(t\right)}},l\left(t\right)=l_0{\displaystyle \frac{K_0}{I_0}}{\displaystyle \frac{I\left(t\right)}{K\left(t\right)}},I\left(t\right)=nY\left(t\right),$

here the parameters with zero indices take their constant value at the initial time.

Relations (4) are transformed into an equation for new variables: capital productivity efficiency: $y\left(t\right)={\displaystyle \frac{Y\left(t\right)}{K\left(t\right)}}$; efficiency of new technologies: $x\left(t\right)={\displaystyle \frac{u\left(t\right)}{U\left(t\right)}}$.

(5)$\dot{y}\left(t\right)=y^2\left(t\right)\left(1-n\right)n\left(x\left(t\right)-x^{\ast }\right),x^{\ast }=2-{\displaystyle \frac{l_0+{\mu }_0}{n_0{y}_0}}.$

Based on the Formula (5) S.V. Dubovsky derived a differential equation to describe the dynamics of the efficiency of new technologies, taking into account the assumption that the growth of the latest technological level is proportional to the latest technological level and its relative efficiency, as well as financial indicators, the rate of accumulation and return on assets.

The above assumption of S.V. Dubovsky is shown in the following equation:

(6)$\dot{U}\left(t\right)=U\left(t\right)\left({\displaystyle \frac{u\left(t\right)}{U\left(t\right)}}-1\right)n\left(a+by\left(t\right)\right),$

where a and b are constant parameters.

Equation (6) divided by $u\left(t\right)$ and taking into account relation (4) allowed S.V. Dubovsky write down the differential equation for the variable $x\left(t\right)$ in the form of Equation (4):

(7)$\dot{x}\left(t\right)=-\lambda x\left(t\right)\left(x\left(t\right)-1\right)n\left(y\left(t\right)-y^{\ast }\right),\lambda =1-b,y^{\ast }={\displaystyle \frac{a}{1-b}}.$

As a result, S.V. Dubovsky managed to reduce the basic model (1)–(3) to a system of two nonlinear differential Equations (5) and (7) with respect to two variables—the efficiency of new technologies $x\left(t\right)$ and capital productivity $y\left(t\right)$. Here, the first term, as a rule, correlates with economic profit, and the second with the depreciation of production assets and the costs of capital owners.

Remark 1.

It is necessary to emphasize that system (5) and (7) has an equilibrium stationary solution: $x\left(t\right)=x^{\ast },y\left(t\right)=y^{\ast }$, which agrees well with real statistics.

The introduction of new variables, and variations relative to the point $x^{\ast },y^{\ast }$ in the form $\delta x=x-x^{\ast },\delta y=y-y^{\ast }$. S.V. Dubovsky obtained linear differential equations in variations from Equations (5) and (7):

(8)${\displaystyle \frac{d\delta y}{dt}}={\left(y^{\ast }\right)}^2\left(1-n\right)n\delta x,{\displaystyle \frac{d\delta x}{dt}}=-\lambda \left(x^{\ast }-1\right)n{x}^{\ast }\delta y.$

Having multiplied Equation (8) by suitable coefficients, S.V. Dubovsky obtained the following relation (first integral):

(9)${\displaystyle \frac{{\left(\delta y\right)}^2}{{\left(y^{\ast }\right)}^2\left(1-n\right)}}+{\displaystyle \frac{{\left(\delta x\right)}^2}{\lambda x^{\ast }\left(x^{\ast }-1\right)}}=C^2.$

Equation (9) obtained by S.V. Dubovsky describes a family of ellipses for different values of the derivative of the constant C.

2.2. Statement of the Problem, Solution Method and Its Properties

Based on Equations (5) and (7) in article [14] S.V. Dubovsky proposed the following Cauchy problem

(10)$\left\{\begin{array}{c}\dot{x}\left(t\right)=-\lambda nx\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right),\hfill \\ \dot{y}\left(t\right)=n\left(1-n\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right),\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

where $x\left(t\right),y\left(t\right)\in C^1[0,T]$ are functions of the efficiency of new technologies and capital productivity, $t\in \left[0,T\right]$ process time, $T>0$ simulation time, $0<n<1$ accumulation rate, and $\lambda >0$ constant coefficient, which is determined from statistics; $x^{\ast }$ and $y^{\ast }$ are constant values that correspond to the rest point of the system, a and b are given constants that determine the initial conditions, $\dot{x}\left(t\right)=dx/dt,\dot{y}\left(t\right)=dy/dt$.

Definition 1.

We will call the system of nonlinear differential Equation (10), which represents the Cauchy problem, the Classical Mathematical Model of S.V. Dubovsky (CMMD).

Remark 2.

Note that the classical mathematical model of S.V. Dubovsky (10) describes free nonlinear fluctuations in the efficiency of new technologies and capital productivity.

Let us show that the dynamical system (10) can have closed-phase trajectories. To conduct this, we introduce the following definition:

Definition 2.

A simply connected region Ω is a region in which all closed paths are homotopic to zero, i.e., they can be continuously contracted to a point. This means that the region does not contain holes or self-intersections.

Let $f_1\left(x,y\right)=-\lambda nx\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right)$ and $f_2\left(x,y\right)=n\left(1-n\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right)$—right sides (10), which are continuous functions, $f=\left(f_1,f_2\right)$ is a vector field defined in the simply connected domain $\mathsf{\Omega }\subset R^2$. Then, the following theorem is true:

Theorem 1

(Bendixson’s theorem [54]). If the divergence of the vector field f does not change sign and is nonzero throughout the entire region Ω, then the phase trajectories of the autonomous system of differential equations $\dot{\mathbf{x}}=f\left(\mathbf{x}\right)$ do not have closed trajectories that lie entirely in Ω.

Theorem 2.

An autonomous system of differential Equation (10) can have closed-phase trajectories.

Proof. 

Using Theorem 1, we can show that the divergence of the vector field f changes its sign and is not identically zero in a simply connected domain.

${\displaystyle \frac{\partial f_1}{\partial x}}=-2\lambda n\left(x-{\displaystyle \frac{1}{2}}\right)\left(y-y^{\ast }\right),{\displaystyle \frac{\partial f_2}{\partial y}}=-2n\left(1-n\right)y\left(x-x^{\ast }\right),$

Then

(11)$\mathrm{div}f={\displaystyle \frac{\partial f_1}{\partial x}}+{\displaystyle \frac{\partial f_2}{\partial y}}=-2n\left(\lambda \left(x-{\displaystyle \frac{1}{2}}\right)\left(y-y^{\ast }\right)-\left(1-n\right)y\left(x-x^{\ast }\right)\right).$

From the economic meaning of the problem parameters ($x\left(t\right)>0,y\left(t\right)>0,\lambda >0$, $0<n<1$) we can conclude that the change in the sign of the function (11) in a simply connected domain $\mathsf{\Omega }\subset R^2$ can occur. Indeed, if you choose

(12)$0.5<x<x^{\ast },0<y<{\displaystyle \frac{\lambda y^{\ast }}{\lambda -1+n}}.$

Then, the divergence of the vector field f will change its sign in the simply connected region. If condition (12) is violated, then the divergence retains its sign. Using the visualization tools of the Maple computer mathematics environment, we will construct an example of a simply connected region in which changes in the signs of the parameter values are performed, that is, the conditions of Theorem 2 are satisfied (Figure 1): $n=0.2,\lambda =2.25,x^{\ast }=1.3,y^{\ast }=0.5$.

Based on the above, we can conclude that the system (10) may have closed phase trajectories. □

The theorem indicates that the classical mathematical model of S.V. Dubovsky (10) can be suitable for studying cyclical regimes that can arise in the theory of economic cycles and crises. Let us construct oscillograms and phase trajectories for the classical model of S.V. Dubovsky for different values of the initial conditions a and parameter $\lambda$ (Figure 2). For this purpose, we will use the one-step Adams–Bashforth–Moulton method (ABM) of the second order of convergence (trapezoidal rule) [55], which was implemented in the MATLAB computer algebra system.

From Figure 2 we see that an increase in the value of the parameter $\lambda$, as well as the values of a ($b=0.5$) leads to an increase in the efficiency of new technologies and capital productivity. This, in turn, leads to an increase in the orbit of the phase trajectory, which is closed (the equilibrium point of the system (10) $\left(x^{\ast },y^{\ast }\right)$ is the center). The periods of phase trajectories according to [14]: T = 50.1, 54.9, 60.6 years; their direction of travel is counterclockwise.

The interpretation of the economic cycle is given according to the phase portrait, by analogy with article [14]. The phase portrait can be represented as four quadrants with the center at the point $\left(1.3,0.5\right)$. In the first quadrant (northeast coordinate angle), due to increased economic activity, the supply of highly effective innovations is depleted, and the system moves into the second quadrant (counterclockwise movement). In the second quadrant, due to the reduced efficiency of innovations, economic activity decreases, and the system moves into the third quadrant. In the third quadrant, due to low economic activity, a stock of highly effective innovations accumulates, and the system moves into the fourth quadrant. In the fourth quadrant, due to the high efficiency of innovations, economic activity begins to grow, and the system moves to the first quadrant, and the cycle repeats.

Remark 3.

It should be noted that the simulated cycles in Figure 2 also have asymmetry. This property of cycles is associated with different times of entry into one or another quadrant due to the fact that the cycles have elliptical orbits.

Remark 4.

Note that according to A.A. Akaev [26], results obtained by S.V. Dubovsky in [14] were somewhat overestimated, and their error reached 50%.

Based on the above, there is a need to clarify the classical model of S.V. Dubovsky (10). Such clarification may be based on the dependence of the rate of accumulation on capital productivity, the introduction of external factors such as investment inflows and new technologies or management decisions, as well as the presence of heredity (memory) effects.

2.3. The Rate of Accumulation as a Function of Capital Productivity

In [14] S.V. Dubovsky pointed out that the rate of accumulation in the system (10) can be a function of capital productivity of the form:

(13)$n\left(y\left(t\right)\right)=l_1\left(l_2-1\right)/y\left(t\right),$

where $l_1>0,l_2>0$ are coefficients that are determined from statistical data as parameters in the regression equation $I\left(t\right)=l_{1}Y\left(t\right)-l_{2}K\left(t\right)$ [56]. However, in his work [14] S.V. Dubovsky did not investigate dependence (13).

Remark 5.

In the articles of the authors [32,33] dependence (13) was studied; oscillograms and phase trajectories were obtained for various values of the model parameters. It is shown that the period of cyclic fluctuations depends on the rate of accumulation.

Remark 6.

It should be noted that in a later work, ref. [25] S.V. Dubovsky suggested using another dependency:

(14) $n\left(y\left(t\right)\right)=l_1-l_2/y\left(t\right),$

where $l_1>0,l_2>0$ are constants that are also determined from statistics [56].

However, the authors’ articles [32,33] did not study the dependence (14). Therefore, in this article we will eliminate this drawback; we will conduct a study of the classical model of S.V. Dubovsky (10) taking into account the dependence (14), and also provide an economic interpretation of the research results.

Let us make a number of comments on the application of functional dependencies (13) and (14) in the classical mathematical model of S.V. Dubovsky.

Remark 7.

Due to the fact that the accumulation rate must be positive, then for dependence (13) the condition $l_2>1$ is satisfied, and for dependence (14) the condition $l_1>l_2/y\left(t\right)$.

Remark 8.

Note that according to dependence (13), with an increase in capital productivity, the accumulation rate decreases. A decrease in the accumulation rate depending on capital productivity means that by increasing the efficiency of using capital productivity, the organization can reduce the amount of its own funds that need to be accumulated to ensure the stable functioning and development of the business. This allows you to optimize your capital structure and improve your return on investment. Reducing the accrual rate can be especially important for organizations with limited financial resources, as it allows for more efficient use of available funds.

Remark 9.

You can also notice that dependence (14) indicates an increase in the rate of accumulation from capital productivity. Increasing the savings rate based on capital productivity means that the higher the level of efficiency in the use of capital productivity, the more funds a company can accumulate for future investments or financial purposes. Increasing capital productivity can lead to increased profits, which, in turn, allows the organization to increase the amount of savings. It is an important aspect of an organization’s financial strategy because it shows how effectively available resources are used to achieve financial goals.

Let us study dependencies (13) and (14). Taking into account dependence (13) or (14), the classical model of S.V. Dubovsky (10) takes the form $\left(y\left(t\right)\ne 0\right)$:

(15)$\left\{\begin{array}{c}\dot{x}\left(t\right)=-\lambda n\left(y\left(t\right)\right)x\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right),\hfill \\ \dot{y}\left(t\right)=n\left(y\left(t\right)\right)\left(1-n\left(y\left(t\right)\right)\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right),\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

The Cauchy problem (15), like the Cauchy problem (10), can have closed phase trajectories. This can be shown using Bendixson’s theorem.

Figure 3 shows examples of simply connected regions in which the divergence of the vector field changes its sign: Figure 3a is the accumulation rate function (13), and Figure 3b is the accumulation rate function (14).

We will seek a solution to problem (15) using the numerical one-step Adams–Bashforth–Moulton method of the second order of convergence (trapezoidal method). The phase trajectories obtained according to the numerical method for model (15) are shown in Figure 4 taking into account dependencies (13).

Figure 4 shows the calculation of oscillograms and phase trajectories for the dependence (13) of the accumulation rate on capital productivity, taking into account different values of $l_1,l_2$ and a. It can be noted that the rate of accumulation increases the asymmetry of the cycle. It can also be noted that an increase in the values of the parameters $l_1,l_2$ affects the frequency and, consequently, the oscillation period. Let the time axis be the number of years. Then, for example, we see that for the following set of values: $l_1=0.8,l_2=1.3,a=1.45$ period fluctuations are $T\approx 40$ years. Indeed, for the efficiency of new technologies $x\left(t\right)$: $T=118.121-77.441=40.68$ years, and for the efficiency of capital productivity $y\left(t\right)$: $T=86.361-47.041=$39.32 years. This period value corresponds to that indicated by A.A. Akaev in article [26], in which the author claims that the duration of N.D. cycles Kondratieva was on average 40 years old in the 20th century.

Figure 5 shows the calculation of oscillograms and phase trajectories for dependence (14) for various values of $l_1,l_2$ and a.

We see that the phase trajectories are closed. The oscillograms in Figure 5 differ from the oscillograms in Figure 4; however, the phase trajectories have similar orbits. Here too, an increase in the values of $l_1,l_2$ and a leads to greater cycle asymmetry, which is also reflected in the oscillation period. The oscillation period, for example, for a cycle with parameters $l_1=1.3,l_2=0.8,a=1.45$ is ≈56 years. Such values of oscillation periods were obtained by S.V. Dubovsky in his work [14]; however, these values were somewhat overestimated compared to the real economic cycle. However, despite this, all of the above indicates that the period of a particular cycle can be estimated based on dependencies (13) and (14), i.e., with the correct choice of parameter values $l_1$ and $l_2$. Let us remember that the values of these parameters are taken from statistics.

In the case when the selected parameter values give an incorrect value for the cycle period, it is necessary to refine the model. For example, take into account some properties of the system under consideration. As such, external influences on the system (10) can be taken into account, which determines the influx of external investments and new technologies or management decisions.

2.4. Influx of External Investments and New Technologies

The influx of external investments and new technologies can be taken into account in the mathematical model (15) by introducing the functions $f_1\left(t\right)$ and $f_2\left(t\right)$ into its right-hand sides

(16)$\left\{\begin{array}{c}\dot{x}\left(t\right)=-\lambda n\left(y\left(t\right)\right)x\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right)+f_1\left(t\right),\hfill \\ \dot{y}\left(t\right)=n\left(y\left(t\right)\right)\left(1-n\left(y\left(t\right)\right)\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right)+f_2\left(t\right),\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

Remark 10.

The mathematical model (16) describes forced nonlinear fluctuations in the efficiency of new technologies and capital productivity.

Remark 11.

In our work we will consider the harmonic functions $f_1\left(t\right)={\delta }_{1}cos\left({\omega }_{1}t\right)$ and $f_2\left(t\right)={\delta }_{2}cos\left({\omega }_{2}t\right)$ in model (16), where ${\delta }_1$ and ${\delta }_2$ are amplitudes, and ${\omega }_1$ and ${\omega }_2$ are frequencies. The choice of such dependencies is due to the fact that the influx of external investments and new technologies is periodic. In addition, such dependences lead to the occurrence of resonance effects in the oscillatory system, which is of great practical importance. In this work, we will demonstrate these properties using computer simulation.

Using the Adams–Bashforth–Moulton numerical algorithm of the second order of convergence for solving the Cauchy problem (16), we can construct oscillograms and phase trajectories for various values of its parameters, taking into account dependence (13).

Figure 6 shows an example of constructing oscillograms and a phase trajectory obtained with the following parameter values: $T=500, \tau =0.01, \lambda =2.25, l_1=0.6$, $l_2=1.1$, $a=1.35, b=0.5$, ${\delta }_1={\delta }_2=0.1, {\omega }_1={\omega }_2=6$. Note that the phase trajectory is closed; in addition, it experiences oscillations due to external harmonic functions.

Figure 7 shows another example with amplitude values ${\delta }_1={\delta }_2=0.5$, other parameters remained unchanged. Here we see an increase in the amplitude of natural oscillations on the oscillograms. The phase trajectory is not closed and has the appearance of an unwinding spiral. Just like in Figure 7, the phase trajectory experiences oscillations due to external harmonic functions.

For dependence (14), both cyclic oscillations (Figure 8) and resonance effects (Figure 9) can also occur.

In Figure 9 we see that the phase trajectory is not closed, it unwinds due to the fact that the amplitude of oscillations increases under resonance conditions. In addition, the trajectory itself experiences vibrations due to external harmonic influence.

In Figure 6, Figure 7, Figure 8 and Figure 9 it seems that the phase trajectory is a three-dimensional figure. This is due to the fact that the boundary of the phase trajectory also experiences oscillations and time t here becomes the third degree of freedom.

It should be noted that the functions $f_1\left(t\right)={\delta }_{1}cos\left({\omega }_{1}t\right)$ and $f_2\left(t\right)={\delta }_{2}cos\left({\omega }_{2}t\right)$ can produce more complex cycles (Figure 10).

These complex cycles are similar to Lissajous figures, the shape of which is determined by the ratio of the frequencies ${\omega }_1$ and ${\omega }_2$. It is known that if the frequency ratio ${\omega }_1/{\omega }_2$ is a rational number, then the trajectory will be closed, and if it is irrational, then it will not be. In any case, the study of such figures deserves attention.

In the works of the authors [30,31,32,33,34], a generalization of the CMMD was proposed taking into account the heredity effect [57]. This property indicates that the system remembers the impact exerted on it and, from a mathematical point of view, it can be described using integro-differential equations or using a fractional derivative [37,38,39,40,41].

It is shown in [31] that the introduction of fractional derivatives into system (10) leads to damped oscillations, and therefore, in order for cyclic oscillations to exist, it is necessary to have an external source of harmonic form $f_1\left(t\right)={\delta }_{1}cos\left({\omega }_{1}t\right)$ and $f_2\left(t\right)={\delta }_{2}cos\left({\omega }_{2}t\right)$. This makes it possible for the phase trajectory to reach a limit cycle after some time. Let us consider these modifications of the CMMD (16).

3. Fractional Mathematical Model S.V. Dubovsky with Constant Heredity

In this section, we will consider a generalization of the CMMD (16) in the case of constant heredity using the apparatus of fractional calculus. Let us consider the formulation of the problem in terms of Gerasimov–Caputo fractional derivatives. Next, a solution technique based on the numerical ABM algorithm will be given. The errors of the method are investigated, and an assessment of the computational accuracy using Runge’s rule is given. A visualization of the simulation results and their interpretation, as well as a software implementation of the numerical algorithm, are presented.

3.1. Some Definitions from the Theory of Fractional Calculus

It should be noted that more detailed information about fractional calculus can be found in well-known monographs [37,38,39,40,41].

Definition 3.

Euler’s gamma function is defined as follows:

(17) $\mathsf{\Gamma }\left(x\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-t}{t}^{x-1}dt,x\in \mathbb{C}:\Re \left(x\right)>0.$

Definition 4.

Let $x\left(t\right)\in C[0,T],T>0$. Integral

(18) $I_{0t}^{q}x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(q\right)}}\underset{0}{\overset{t}{\int }}{\displaystyle \frac{x\left(\tau \right)d\tau }{{\left(t-\tau \right)}^{1-q}}},$

where $q>0$, $\mathsf{\Gamma }\left(·\right)$ is the Euler gamma function (17), called the Riemann–Liouville fractional integral.

The fractional Riemann–Liouville integral (18) has the following properties:

(19)$1)I_{0t}^{0}x\left(t\right)=x\left(t\right); 2)I_{0t}^{\alpha }{I}_{0t}^{\beta }x\left(t\right)=I_{0t}^{\beta }{I}_{0t}^{\alpha }=I_{0t}^{\alpha +\beta }x\left(t\right).$

Definition 5.

The Gerasimov–Caputo fractional derivative of the constant order $0<q<1$ of the function $x\left(t\right)\in C^1[0,T],T>0$ has the following form:

(20) ${\partial }_{0t}^{q}x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(1-q\right)}}\underset{0}{\overset{t}{\int }}{\displaystyle \frac{\dot{x}\left(\tau \right)d\tau }{{\left(t-\tau \right)}^q}},$

where $\mathsf{\Gamma }\left(·\right)$ is the Euler gamma function (17).

The Riemann–Liouville fractional integral (18) and the Gerasimov–Caputo fractional derivative (20) have the semigroup property:

(21)${\partial }_{0t}^q{I}_{0t}^{q}x\left(t\right)=x\left(t\right).$

Remark 12.

It should be noted that the operator (20) began to be called Caputo’s fractional derivative and denoted as ${}^{C}D_{0t}^{\alpha }x\left(t\right)$, thanks to articles [58,59]. This mathematical construction was presented by the Italian mathematician M. Caputo in 1967 in work [60]. However, earlier the Soviet mechanic A.N. Gerasimov in 1948, in his work [61], devoted to plasticity problems within the framework of hereditary mechanics [62], introduced a partial fractional derivative of the order $0<\alpha <1$:

$D_{-\infty }^{\alpha }u\left(x,t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(1-\alpha \right)}}\underset{-\infty }{\overset{t}{\int }}{\displaystyle \frac{u_{\tau }\left(x,\tau \right)d\tau }{{\left(t-\tau \right)}^{\alpha }}}.$

Therefore, we will call the fractional derivative operator (20) the Gerasimov–Caputo fractional derivative.

Remark 13.

Note that the listed Formulas (18) and (20) are valid for functions of many variables. In this case, we will work with partial derivatives of fractional order.

Remark 14.

There are other definitions of fractional derivatives and integrals. Some are modifications of the Riemann–Liouville or Gerasimov–Caputo operators. Others are defined differently. The question of choosing one or another fractional operator to solve an applied problem is open. Fractional operators are usually chosen for simplicity in mathematical transformations and interpretation of simulation results. For example, for an equation with the Gerasimov–Caputo operator, traditional initial and boundary conditions are set, which is important for physical applications. The Riesz operator is convenient for the Fourier integral transform.

3.2. Problem Statement and Solution Method

Let us consider a generalization of the CMMD (16), which according to (20) can be represented as the following Cauchy problem:

(22)$\left\{\begin{array}{c}{\partial }_{0t}^{\alpha }x\left(t\right)=-\lambda n\left(y\right)x\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right)+f_1\left(t\right),0<\alpha <1,\hfill \\ {\partial }_{0t}^{\beta }y\left(t\right)=n\left(y\right)\left(1-n\left(y\right)\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right)+f_2\left(t\right),0<\beta <1,\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

Remark 15.

Note that in the limiting case when $\alpha =\beta =1$ we move from the Cauchy problem (22) to the classical Cauchy problem (16).

Remark 16.

We will call the Cauchy problem (22) the fractional mathematical model of S.V. Dubovsky (FMMD).

We will not dwell here in detail on the derivation of the model equations in (22). We will indicate only some approaches. The first approach is to use the semigroup property (21) [63]. For example, if the Cauchy problem is written in the general form as follows:

(23)$\left\{\begin{array}{c}\dot{x}\left(t\right)=I_{0t}^{\alpha }{f}_1\left(x,y,t\right),0<\alpha <1,\hfill \\ \dot{y}\left(t\right)=I_{0t}^{\beta }{f}_2\left(x,y,t\right),0<\beta <1,\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

Then, act on the first Equation (23) on the left with the operator ${\partial }_{0t}^{\alpha }$, and for the second equation on ${\partial }_{0t}^{\beta }$, we arrive at the problem Cauchy form (22) with specific functions $f_1\left(x,y,t\right)$ and $f_2\left(x,y,t\right)$.

The second approach is based on the formal replacement of integer derivatives from the right side of the system (16) with derivatives of fractional orders [64]:

(24)${\displaystyle \frac{d}{dt}}\to {\displaystyle \frac{1}{{\theta }_{\alpha }^{1-\alpha }}}{\partial }_{0t}^{\alpha },{\displaystyle \frac{d}{dt}}\to {\displaystyle \frac{1}{{\theta }_{\beta }^{1-\beta }}}{\partial }_{0t}^{\beta },$

where ${\theta }_{\alpha },{\theta }_{\beta }$ are parameters that have time dimensions; they are necessary to comply with the dimensions in the model equations. You can always introduce the characteristic time of the process, perform dimensionlessness and, in principle, arrive at system (22).

There is an approach associated with hereditary (hereditary) equations, in the terminology of Vito Volterra [57]. For example, we can write the original Cauchy problem in the integrodifferential formulation:

(25)$\left\{\begin{array}{c}\underset{0}{\overset{t}{\int }}{K}_1\left(t-\tau \right)x\left(\tau \right)d\tau =-\lambda n\left(y\right)x\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right)+f_1\left(t\right),\hfill \\ \underset{0}{\overset{t}{\int }}{K}_2\left(t-\tau \right)y\left(\tau \right)d\tau =n\left(y\right)\left(1-n\left(y\right)\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right)+f_2\left(t\right),\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

where $K_1\left(t-\tau \right)$ and $K_2\left(t-\tau \right)$ are difference kernels or heredity functions. Further, according to Volterra’s hereditary principles, the memory functions in (25) should decay as $\tau \to -\infty$. So, based on this, we can select memory functions as follows:

(26)$K_1\left(t-\tau \right)={\displaystyle \frac{{\left(t-\tau \right)}^{-\alpha }}{\mathsf{\Gamma }\left(1-\alpha \right)}},K_2\left(t-\tau \right)={\displaystyle \frac{{\left(t-\tau \right)}^{-\beta }}{\mathsf{\Gamma }\left(1-\beta \right)}}.$

Substituting relation (26) into the system (25) taking into account (20), we can obtain system (22).

Remark 17.

Note that in [65] it is noted that if the memory function is chosen in the form $K\left(t-\tau \right)=\omega \delta \left(t-\tau \right)$, where ω is a positive constant, $\delta \left(·\right)$ is the Dirac function, then we obtain the absence of memory by analogy with the Markov process. The case when in a dynamic system the memory function has the form $K\left(t-\tau \right)=t-H\left(t-\tau \right)$, where $H\left(·\right)$ is the function Heaviside, leads to systems with complete memory availability. The power-law nature of memory functions in the form (26) leads to dynamic systems with “partial loss of memory”.

Due to the nonlinearity of the Cauchy problem (22), it is necessary to use numerical methods. We will use the one-step Adams–Bashforth–Moulton method, adapted for solving fractional order differential equations. For the Cauchy problem (22), a generalization of the numerical ABM method is proposed [66,67,68]. To conduct this, we introduce a uniform finite-difference grid with a constant step $\tau =T/N$, as well as grid functions $x_{k+1},y_{k+1}$, $x_{k+1}^p,y_{k+1}^p$, $n_k=n\left(y_k\right)$, ${\alpha }_1=\alpha ,{\alpha }_2=\beta$, for $k=0,...,N-1$. A modification of the ABM method takes the following form. For the predictor, we obtain the calculation formulas:

(27)$\left\{\begin{array}{c}{x}_{k+1}^p=x_0+{\displaystyle \frac{{\tau }_1^{\alpha }}{\mathsf{\Gamma }\left({\alpha }_1+1\right)}}\sum _{j=0}^k{\theta }_{j,k+1}^1\left(-\lambda n_j{x}_j\left(x_j-1\right)\left(y_j-y^{\ast }\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right),\hfill \\ y_{k+1}^p=y_0+{\displaystyle \frac{{\tau }_2^{\alpha }}{\mathsf{\Gamma }\left({\alpha }_2+1\right)}}\sum _{j=0}^k{\theta }_{j,k+1}^2\left(n_j\left(1-n_j\right)y_j^2\left(x_j-x^{\ast }\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right),\hfill \\ {\theta }_{j,k+1}^i={\left(k-j+1\right)}^{{\alpha }_i}-{\left(k-j\right)}^{{\alpha }_i},i=1,2.\hfill \end{array}\right.$

Then, using the Adams–Moulton formula for the corrector, we obtain the following:

(28)$\left\{\begin{array}{c}{x}_{k+1}=x_0+{\displaystyle \frac{{\tau }_1^{\alpha }}{\mathsf{\Gamma }\left({\alpha }_1+2\right)}}\left(\begin{array}{c}\left(-\lambda n_{k+1}{x}_{k+1}^p\left(x_{k+1}^p-1\right)\left(y_{k+1}^p-y^{\ast }\right)+{\delta }_{1}cos\left({\omega }_1\tau \left(k+1\right)\right)\right)+\hfill \\ +\sum _{j=0}^k{\rho }_{j,k+1}^1\left(-\lambda n_j{x}_j\left(x_j-1\right)\left(y_j-y^{\ast }\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right)\hfill \end{array}\right),\hfill \\ y_{k+1}=y_0+{\displaystyle \frac{{\tau }_2^{\alpha }}{\mathsf{\Gamma }\left({\alpha }_2+2\right)}}\left(\begin{array}{c}\left(n_{k+1}\left(1-n_{k+1}\right){\left(y_{k+1}^p\right)}^2\left(x_{k+1}^p-x^{\ast }\right)+{\delta }_{2}cos\left({\omega }_2\tau \left(k+1\right)\right)\right)+\hfill \\ +\sum _{j=0}^k{\rho }_{j,k+1}^2\left(n_j\left(1-n_j\right)y_j^2\left(x_j-x^{\ast }\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right)\hfill \end{array}\right),\hfill \end{array}\right.$

where the weighting coefficients in (28) are determined by the following formula:

${\rho }_{j,k+1}^i=\left\{\begin{array}{c}{k}^{{\alpha }_i+1}-\left(n-{\alpha }_i\right){\left(k+1\right)}^{{\alpha }_i},j=0,\hfill \\ {\left(k-j+2\right)}^{{\alpha }_i+1}+{\left(k-j\right)}^{{\alpha }_i+1}-2{\left(k-j+1\right)}^{{\alpha }_i+1},1\le j\le k,\hfill \\ 1,j=k+1,\hfill \\ i=1,2.\hfill \end{array}\right.$

Remark 18.

Formulas (27) and (28) with values ${\alpha }_1={\alpha }_2=1$ transform into formulas for the classical ABM method.

3.3. Error Analysis of the Adams–Bashforth–Moulton Method

This section presents lemmas and the error theorem for the modified ABM method (27) and (28), more detailed proofs of which can be found in [66].

Lemma 1.

If the function $z\left(t\right)\in C^1\left[0,T\right]$, then the following estimate holds [66]:

(29) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}z\left(t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}z\left(t_j\right)\right|\le C{∥z^{\prime }∥}_{\infty }{t}_{k+1}^{\alpha }\tau ,$

where $\parallel {\parallel }_{\infty }$, Chebyshev norm, $C=1/\alpha$.

Lemma 2.

If the function $z\left(t\right)\in C^2\left[0,T\right]$, then there is a constant $C_{\alpha }$ that depends only from α such that [66]:

(30) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}z\left(t\right)dt-\sum _{j=0}^k{\rho }_{j,k+1}z\left(t_j\right)\right|\le C_{\alpha }{∥z^{″}∥}_{\infty }{t}_{k+1}^{\alpha }{\tau }^2,$

Based on Lemma 1 and Lemma 2, we formulate the following theorem.

Theorem 3.

If ${\partial }_{0t}^{{\alpha }_i}{x}_i\left(t\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)$ [67] and

(31) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{{\alpha }_i-1}{\partial }_{0t}^{{\alpha }_i}{x}_i\left(t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}^i{\partial }_{0t}^{{\alpha }_i}{x}_i\left(t_j\right)\right|\le C_1{t}_{k+1}^{{\alpha }_i}\tau ,$

(32) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{{\alpha }_i-1}{\partial }_{0t}^{{\alpha }_i}{x}_i\left(t\right)dt-\sum _{j=0}^k{\rho }_{j,k+1}^i{\partial }_{0t}^{{\alpha }_i}{x}_i\left(t_j\right)\right|\le C_2{t}_{k+1}^{{\alpha }_i}{\tau }^2,$

where $C_1$ and $C_2$ are constants, then the Adams–Bashforth–Moulton numerical method converges with the order

(33) $\underset{1\le j\le k}{max}\left|x_i\left(t_j\right)-x_{i,j}\right|=O\left({\tau }^q\right),q=1+\underset{i}{min}{\alpha }_i.$

Remark 19.

Note that in the case when ${\alpha }_1={\alpha }_2=1$ we have the second order of convergence according to estimate (33), which corresponds to the order of convergence of the classical ABM method.

Next, using test examples, we will evaluate the computational accuracy of the ABM method (27) and (28) using Runge’s rule (double recalculation method) [69]. For these purposes, we introduce the following calculation formulas:

(34)${\xi }_x^i={\displaystyle \frac{max\left|x_i-x_{2i}\right|}{2^q-1}},{\xi }_y^i={\displaystyle \frac{max\left|y_i-y_{2i}\right|}{2^q-1}},$

where $q=1+min\left({\alpha }_1,{\alpha }_2\right)$ is the theoretical order of accuracy of the ABM method according to estimate (32), ${\xi }_x^i$, ${\xi }_y^i$—error of the ABM method, $x_i,y_i$—numerical solution at step ${\tau }_i$, $x_{2i},y_{2i}$—numerical solution at step ${\tau }_i/2$. The computational accuracy p will then be determined by the following formulas:

(35)$p_x={log}_2\left({\displaystyle \frac{{\xi }_x^i}{{\xi }_x^{i+1}}}\right),p_y={log}_2\left({\displaystyle \frac{{\xi }_y^i}{{\xi }_y^{i+1}}}\right),$

where ${\tau }_i,{\tau }_{i+1}={\tau }_i/2$ computational grid steps, ${\xi }_x^{i+1}$, ${\xi }_y^{i+1}$ solution errors at step ${\tau }_{i+1}$.

Example 1

(Case of CMMD). Let us consider a test example in the case ${\alpha }_1={\alpha }_2=1$. We choose the parameter values for the Cauchy problem (22) as follows: ${\delta }_1={\delta }_2=1, {\omega }_1={\omega }_2=0.5$, $n=0.2, \lambda =1.5, x\left(0\right)=5,$ $y\left(0\right)=4, x^{\ast }=1.35, y^{\ast }=0.5, t\in \left[0.1\right]$. The results of the numerical analysis are shown in Table 1.

From Table 1 it is clear that with an increase in the calculated grid nodes, the computational accuracies $p_x,p_y$ tend to $q=2$, the order of accuracy of the classical ABM method.

Example 2

(Case of FMMD). Let us consider the case ${\alpha }_1=0.9,{\alpha }_2=0.8$, and take the remaining parameters from the previous example. The simulation results are shown in Table 2.

From Table 2 it is clear that with an increase in the calculated grid nodes, the computational accuracies $p_x,p_y$ tend to $q=1+min\left(0.9,0.8\right)=1.8$—the order of accuracy of the ABM method (27) and (28).

3.4. Simulation Results

Let us look at some examples of simulation results within the FMMD model and their visualization. The construction of oscillograms and phase trajectories was carried out using the HFMD 1.0 computer program written in Phyton [33].

Example 3

(Case of CMMD). Let us consider the case when in the numerical ABM algorithm (27), (28) and the parameter values are as follows: ${\alpha }_1={\alpha }_2=1$, ${\delta }_1={\delta }_2=0$ and n is a constant. This leads us to the model (10), which was considered by S.V. Dubovsky in his work [14]. Let us choose the values of the remaining parameters: $n=0.2,x^{\ast }=1.3,y^{\ast }=0.5,T=100$, $\tau =0.2,b=0.5$. Next, using the numerical ABM algorithm (27), (28), we construct oscillograms and phase trajectories for various values of $\lambda =\left[2.25,2.5,2.75\right]$ and $a=\left[1.35,1.4,1.45\right]$.

We see that in Figure 11 the results coincided with the previously obtained results in Figure 2. This once again confirms that the fractional mathematical model of S.V. Dubovsky is a generalization of the classical mathematical model [14].

Example 4

(Case of FMMD). Let ${\delta }_1={\delta }_2=0$ and $n=0.2,\lambda =2.25,x^{\ast }=1.3,y^{\ast }=0.5$, $T=500,\tau =0.04$, $a=1.35,b=0.5$. In this case, the oscillatory system experiences its own vibrations. Let us construct oscillograms and phase trajectories for various values of ${\alpha }_1$ with a fixed value of ${\alpha }_2=1$.

In Figure 12 we see that as the values of the parameter ${\alpha }_1$ decrease; the oscillations dampen (dissipation) and on the phase plane the rest point $\left(x^{\ast },y^{\ast }\right)$ is a stable focus. Let us give another example. Let us fix the value ${\alpha }_1$ and change the values ${\alpha }_2=\left[0.8,0.9,1\right]$.

As in Figure 12 and Figure 13 we see damped oscillations when changing the values of the parameter ${\alpha }_2$, the rest point $\left(x^{\ast },y^{\ast }\right)$ is also a stable focus (phase trajectories, twisting spirals).

Based on the damped oscillatory mode, we can conclude that the presence of memory in the system (22), provided there is no external influence, does not lead to a closed phase trajectory, i.e., there is no self-oscillating mode in the system. Therefore, the existence of a closed trajectory is possible only under external influence. Let us now consider an example taken from [32].

Example 5.

Fractional mathematical model S.V. Dubovsky with parameter values [29]: $n=0.2,\lambda =2.25,$ $x^{\ast }=1.3,y^{\ast }=0.5,T=500,\tau =0.04$, ${\alpha }_1={\alpha }_2=0.1$, ${\delta }_1={\omega }_1=0$, ${\delta }_2=0.5,{\omega }_2=1.5$, $a=1.35,b=0.5$.

It can be noted that the phase trajectory in Figure 14 goes to the limit cycle. Let us construct phase trajectories depending on different values of the initial conditions.

Figure 15 shows a family of phase trajectories obtained for different initial values of $\left(a,b\right)$. These phase trajectories show that the limit cycle presented in Figure 14 is not stable.

Definition 6.

A limit cycle is called stable if all phase trajectories starting in the ε-neighborhood asymptotically approach the limit cycle as $t\to \infty$. Limit cycles in which close phase trajectories approach them indefinitely are called attractors. Limit cycles in which close phase trajectories repel from them are called repellers.

Indeed, if we take the initial values inside this limit cycle, for example, the value $\left(1.4,0.45\right)$, then the phase trajectory goes to another limit cycle with a larger orbit, the red phase trajectory in Figure 15. Next, we can see that with a further increase in the initial value along the x-axis and a decrease in the initial value along the ordinate axis, we obtain another limit cycle again with a larger orbit, the blue phase trajectory in Figure 15.

Example 6.

Fractional mathematical model S.V. Dubovsky with comparable parameter values: ${\alpha }_1={\alpha }_2=0.4,{\alpha }_1={\alpha }_2=0.6,{\alpha }_1={\alpha }_2=0.8$. We will select the values of the remaining parameters from Example 5.

Figure 16 shows a visualization of the simulation results for various parameter values: ${\alpha }_1={\alpha }_2=0.4,{\alpha }_1={\alpha }_2=0.6,{\alpha }_1={\alpha }_2=0.8$. We see that all phase trajectories reach limit cycles. As the values of the parameters ${\alpha }_1$ and ${\alpha }_2$ decrease, the limit cycles have an increasingly larger orbit. In addition, we see that as the values of the parameters ${\alpha }_1$ and ${\alpha }_2$ decrease, the time to reach the limit cycle also increases. This indicates that a memory effect occurs in the system (22), i.e., it takes time to stabilize the oscillatory modes.

Let us consider the following more general case (Figure 17). Let us choose the following AZ parameter values: ${\delta }_1=1,{\omega }_1=3,{\delta }_2=0.5,{\omega }_2=1.5$, and leave the remaining parameters unchanged.

The presence of external influence in the form of harmonic functions, taking into account heredity, significantly complicates the nature of oscillations and the shape of phase trajectories (Figure 17). The shape of the phase trajectories is more similar to Lissajous figures, which are found in technology and depend on the frequency ratio ${\omega }_1/{\omega }_2$. However, due to the memory effect, these figures are quite strongly deformed, for example, (see Figure 10). We also see that the trend from previous examples remains the same: as the values of ${\alpha }_1$ and ${\alpha }_2$ decrease, the time to reach the limit cycle increases, as well as its orbit.

Let us now consider the fractional mathematical model of S.V. Dubovsky, when the accumulation rate n is a function of capital productivity $y\left(t\right)$ and changes according to law (13), where $l_1=0.6$ and $l_2=1.1$. Let us construct phase trajectories for various values of ${\alpha }_1,{\alpha }_2$ and ${\delta }_1=3,{\delta }_2=1,{\omega }_1=1,{\omega }_2=1$. The values of the remaining parameters will be left unchanged.

We see that when external influence is turned on with parameter values ${\delta }_1=3$, ${\delta }_2=1,{\omega }_1=1,{\omega }_2=1$ the phase trajectories reach the limit cycle. Note that the shape of limit cycles differs from those previously obtained; this is due to the influence of the accumulation rate as a function of capital productivity. Moreover, in this inversely proportional relationship (13), the greater the capital productivity, the lower the accumulation rate. Here, you can also notice that the efficiency of new technologies $x\left(t\right)$ and the efficiency of fixed assets $y\left(t\right)$ behave the same way: they either increase together or decrease together.

Figure 19 shows the phase trajectories for various values of the parameters ${\alpha }_1,{\alpha }_2$, as well as $l_1$ and $l_2$. We will take the remaining parameter values from the previous example. It can be seen that in this case, the general trend continues, as in the previous case. However, we see that the phase trajectories are different here: they are more deformed and the orbits are more enlarged. In Figure 18 and Figure 19, from an economic point of view, you can notice the following: first, the system needs to reach the limit cycle for some time; at the limit cycle there is a simultaneous increase in the efficiency of new technologies and the efficiency of fixed assets; after some time, there is a simultaneous decline in these indicators, which intensifies as a result of an increase in the accumulation rate; the latter means that the growth of capital productivity is slowed down and its decline is accelerating.

Figure 20 shows another example of constructing a phase trajectory for the values of the parameters ${\alpha }_1,{\alpha }_2$ with fixed values of $l_1$ and $l_2$ for dependence (14). Here, we see that the phase trajectories reach a limit cycle similar to Figure 17. Such limit cycles deserve special attention and economic interpretation.

4. Fractional Mathematical Model S.V. Dubovsky with Variable Heredity

This section proposes the next stage of modernization of FMMD—taking into account variable heredity. This means that the economic system in question has memory effects that change over time. From a mathematical point of view, the change in memory over time can be described using derivatives of fractional order variables.

4.1. Definition of Fractional Derivative of Variable Order

Fractional derivatives of variable orders are studied within the framework of the theory of fractional calculus. There is a well-known review article [70], in which the authors covered in more detail various definitions of the derivative of a fractional variable order and its applications. As one of the applications, we can cite article [71], where a description of the financial system was given using a derivative of a fractional variable order.

Definition 7.

Fractional derivative of variable order of Gerasimov–Caputo type for $x\left(t\right)\in C^1[0,T]$, $T>0$:

(36) ${\partial }_{0t}^{q\left(t\right)}x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(1-q\left(t\right)\right)}}\underset{0}{\overset{t}{\int }}{\displaystyle \frac{\dot{x}\left(\tau \right)d\tau }{{\left(t-\tau \right)}^{q\left(t\right)}}},$

where $0<q\left(t\right)<1$ is a known function from the class $C\left[0,T\right]$.

Remark 20.

In the case when in definition (36) the function $q\left(t\right)$ is a constant, then we come to the definition of the Gerasimov–Caputo fractional derivative (20).

Remark 21.

Similarly, one can introduce the definition of a fractional integral of variable order of the Riemann–Liouville type, as for constant order (18). However, it should be noted here that the semigroup property between the variable-order fractional derivative operators and the variable-order fractional integral is violated [72,73].

4.2. Problem Statement and Solution Method

Based on Remark 21 and taking into account Definition 7, we can write the following Cauchy problem:

(37)$\left\{\begin{array}{c}{\partial }_{0t}^{\alpha \left(t\right)}x\left(t\right)=-\lambda n\left(y\right)x\left(t\right)\left(x\left(t\right)-1\right)\left(y\left(t\right)-y^{\ast }\right)+f_1\left(t\right),\hfill \\ {\partial }_{0t}^{\beta \left(t\right)}y\left(t\right)=n\left(y\right)\left(1-n\left(y\right)\right)y^2\left(t\right)\left(x\left(t\right)-x^{\ast }\right)+f_2\left(t\right),\hfill \\ x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$

Remark 22.

Note that in the case when the functions $\alpha \left(t\right)$ and $\beta \left(t\right)$ are constants, then the Cauchy problem (37) becomes the Cauchy problem (22).

Remark 23.

We will call the Cauchy problem (37) the fractional mathematical model of S.V. Dubov variable order (FMMD VO).

Methods for numerical analysis of equations with derivatives of fractional order variables are quite well developed, some of them can be studied in articles [74,75,76,77].

For the Cauchy problem (37), a generalization of the numerical ABM method (27), (28) is proposed [73]. To conduct this, we introduce a uniform finite-difference grid with a constant step $\tau =T/N$, as well as grid functions: $x_{k+1},y_{k+1}$, $x_{k+1}^p,y_{k+1}^p$, ${\alpha }_{1,k+1}=\alpha \left(t_{k+1}\right)$, ${\alpha }_{2,k+1}=\beta \left(t_{k+1}\right)$, $n\left(y_k\right)=n_k$ for $k=0,...,N-1$. A modification of the ABM method takes the following form:

(38)$\left\{\begin{array}{c}{x}_{k+1}^p=x_0+{\displaystyle \frac{{\tau }^{{\alpha }_{1,k+1}}}{\mathsf{\Gamma }\left({\alpha }_{1,k+1}+1\right)}}\sum _{j=0}^k{\theta }_{j,k+1}^1\left(-\lambda n_j{x}_j\left(x_j-1\right)\left(y_j-y^{\ast }\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right),\hfill \\ y_{k+1}^p=y_0+{\displaystyle \frac{{\tau }^{{\alpha }_{2,k+1}}}{\mathsf{\Gamma }\left({\alpha }_{2,k+1}+1\right)}}\sum _{j=0}^k{\theta }_{j,k+1}^2\left(n_j\left(1-n_j\right)y_j^2\left(x_j-x^{\ast }\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right),\hfill \\ {\theta }_{j,k+1}^i={\left(k-j+1\right)}^{{\alpha }_{i,k+1}}-{\left(k-j\right)}^{{\alpha }_{i,k+1}},i=1,2.\hfill \end{array}\right.$

Then, using the Adams–Moulton formula for the corrector, we obtain the following:

(39)$\left\{\begin{array}{c}{x}_{k+1}=x_0+{\displaystyle \frac{{\tau }^{{\alpha }_{1,k+1}}}{\mathsf{\Gamma }\left({\alpha }_{1,k+1}+2\right)}}\left(\begin{array}{c}\left(-\lambda n_{k+1}{x}_{k+1}^p\left(x_{k+1}^p-1\right)\left(y_{k+1}^p-y^{\ast }\right)+{\delta }_{1}cos\left({\omega }_1\tau \left(k+1\right)\right)\right)+\hfill \\ +\sum _{j=0}^k{\rho }_{j,k+1}^1\left(-\lambda n_j{x}_j\left(x_j-1\right)\left(y_j-y^{\ast }\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right)\hfill \end{array}\right),\hfill \\ y_{k+1}=y_0+{\displaystyle \frac{{\tau }^{{\alpha }_{2,k+1}}}{\mathsf{\Gamma }\left({\alpha }_{2,k+1}+2\right)}}\left(\begin{array}{c}\left(n_{k+1}\left(1-n_{k+1}\right){\left(y_{k+1}^p\right)}^2\left(x_{k+1}^p-x^{\ast }\right)+{\delta }_{2}cos\left({\omega }_2\tau \left(k+1\right)\right)\right)+\hfill \\ +\sum _{j=0}^k{\rho }_{j,k+1}^2\left(n_j\left(1-n_j\right)y_j^2\left(x_j-x^{\ast }\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right)\hfill \end{array}\right),\hfill \end{array}\right.$

where the weighting coefficients in (39) are determined by the following formula:

${\rho }_{j,k+1}^i=\left\{\begin{array}{c}{k}^{{\alpha }_{i,k+1}+1}-\left(n-{\alpha }_{i,k+1}\right){\left(k+1\right)}^{{\alpha }_{i,k+1}},j=0,\hfill \\ {\left(k-j+2\right)}^{{\alpha }_{i,k+1}+1}+{\left(k-j\right)}^{{\alpha }_{i,k+1}+1}-2{\left(k-j+1\right)}^{{\alpha }_{i,k+1}+1},1\le j\le k,\hfill \\ 1,j=k+1,\hfill \\ i=1,2.\hfill \end{array}\right.$

4.3. Error Analysis of the Modified Adams–Bashforth–Moulton Method

In this section, we present lemmas and a theorem, which are a generalization of the results from Section 3.3; proofs of lemmas and theorems are carried out similarly to the method proposed in Section 3.3 [71].

Lemma 3.

If the function $z\left(t\right)\in C^1\left[0,T\right]$, then the following estimate holds:

(40) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{{\alpha }_{k+1}-1}z\left(t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}z\left(t_j\right)\right|\le C{∥z^{\prime }∥}_{\infty }{t}_{k+1}^{{\alpha }_{k+1}}\tau ,$

where $\parallel {\parallel }_{\infty }$—Chebyshev norm, $C=1/{\alpha }_{k+1}$.

Lemma 4.

If the function $z\left(t\right)\in C^2\left[0,T\right]$, then there is a constant $C_{{\alpha }_{k+1}}$, which depends only on ${\alpha }_{k+1}$ such that

(41) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{{\alpha }_{k+1}-1}z\left(t\right)dt-\sum _{j=0}^k{\rho }_{j,k+1}z\left(t_j\right)\right|\le C_{{\alpha }_{k+1}}{∥z^{″}∥}_{\infty }{t}_{k+1}^{{\alpha }_{k+1}}{\tau }^2,$

Theorem 4.

If ${\partial }_{0t}^{{\alpha }_{i,k+1}}{x}_i\left(t\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)$ and

(42) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{{\alpha }_{i,k+1}-1}{\partial }_{0t}^{{\alpha }_{i,k+1}}{x}_i\left(t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}{\partial }_{0t}^{{\alpha }_{i,k+1}}{x}_i\left(t_j\right)\right|\le C_1{t}_{k+1}^{{\alpha }_{i,k+1}}\tau ,$

(43) $\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{{\alpha }_{i,k+1}-1}{\partial }_{0t}^{{\alpha }_{i,k+1}}{x}_i\left(t\right)dt-\sum _{j=0}^k{\rho }_{j,k+1}{\partial }_{0t}^{{\alpha }_{i,k+1}}{x}_i\left(t_j\right)\right|\le C_2{t}_{k+1}^{{\alpha }_{i,k+1}}{\tau }^2,$

then the ABM method (38) and (39) converges with the order

(44) $\underset{1\le j\le k}{max}\left|x_i\left(t_j\right)-x_{i,j}\right|=O\left({\tau }^q\right),q=1+\underset{i}{min}{\alpha }_{i,m},{\alpha }_{i,m}=\underset{\tau \in \left[0,T\right]}{min}{\alpha }_i\left(\tau \right).$

Let us evaluate the computational accuracy of the ABM method (38), (39) using Runge’s rule (34) and (35), which is compared with the theoretical one. A good agreement was obtained. In the test examples, we used linearly monotonically decreasing functions ${\alpha }_1\left(t\right)$ and ${\alpha }_2\left(t\right)$ of the form

(45)${\alpha }_1\left(t\right)={\alpha }_{10}-{\displaystyle \frac{M_{1}t}{T}},{\alpha }_2\left(t\right)={\alpha }_{20}-{\displaystyle \frac{M_{2}t}{T}}.$

where ${\alpha }_{10},{\alpha }_{20},M_1,M_2$ are given constants whose values satisfy the inequalities $0<{\alpha }_1\left(t\right),{\alpha }_2\left(t\right)<1$.

Example 7.

$n=0.2,{\alpha }_{10}=1,{\alpha }_{20}=0.8,t\in \left[0.1\right],x_0=2,y_0=1,\lambda =2.25,{\phi }_1=1,{\phi }_2=0.5,$ ${\delta }_1=0.8,{\delta }_2=1.3,M_1=0.1,M_2=0.3$. The theoretical order of convergence of the ABM method: $q=1.500468750$.

From the Table 3 we see that the computational accuracy approaches the theoretical q in three iterations according to the ABM method (38) and (39).

Example 8.

$l_1=0.6,l_2=1.1\left(14\right),{\alpha }_{10}=0.7,{\alpha }_{20}=0.6,$ ${\delta }_1=0.8,{\delta }_2=1.3,M_1=0.1$, $M_2=0.3$, $t\in \left[0.1\right],x_0=2,y_0=1,\lambda =2.25,{\phi }_1=1,{\phi }_2=0.5$. The theoretical order of convergence of the ABM method: $q=1.300468750$.

From the Table 4 we see that the computational accuracy approaches the theoretical q in five iterations according to the ABM method.

According to Examples 7 and 8, we can conclude that the estimate of the computational accuracy tends to be the theoretical one within a finite number of iterations of the ABM method. The latter indicates the adequacy of the estimate (44).

4.4. Simulation Results

Let us visualize the simulation results for different parameter values. To conduct this, we will use procedures written in MATLAB, which are included in the FDDSVO computer program (Appendix A).

Example 9.

Let us consider the case when the orders ${\alpha }_1\left(t\right)$ and ${\alpha }_2\left(t\right)$ of fractional derivatives linearly decrease on $\left[0.1\right]$ according to Formula (45) [33]. Parameter values: ${\delta }_1={\omega }_1=0$, ${\delta }_2=0.5$, ${\omega }_2=1.5$, $n=0.2$, $\lambda =2.25$, $x^{\ast }=1.3,y^{\ast }=0.5$, $a=1.35,b=0.5$. For mathematical modeling, we choose the following: $T=500,\tau =0.4$.

To find a numerical solution to system (37), a function was written in MATLAB (Appendix A.1) [33] which outputs oscillograms and phase trajectories based on the numerical solution of system (37) using the ABM method (38) and (39). The application of the function for different values of parameters ${\alpha }_{10}=\left(0.8,0.6,0.4\right)$ and ${\alpha }_{20}=\left(0.7,0.5,0.3\right)$, which are included in Formula (45).

Figure 21 shows oscillograms for Example 9. Here, we can observe that over time the oscillations reach a constant level, and the rest point $\left(x^{\ast },y^{\ast }\right)$ for the corresponding closed-phase trajectory is called the limit cycle. It can be noted that as the values of ${\alpha }_1\left(t\right)$ and ${\alpha }_2\left(t\right)$ decrease, the exit time increases to the limit cycle, and the orbit of the limit cycle itself becomes larger.

In Figure 22 we see that the phase trajectories reach limit cycles, which are located quite close to each other. Let us consider the case when the parameter values ${\alpha }_{10}$, ${\alpha }_{20}$, $M_1$, $M_2$ are fixed in Formula (45), and the values of the parameters ${\omega }_2$ and $\lambda$ vary.

In Figure 23 we see that increasing the values of the parameter $\lambda$ and the frequency of external influence ${\omega }_2$ leads to a decrease in the orbits of limit cycles. Let us consider the case when the parameters ${\delta }_1$ and ${\omega }_1$ change (Figure 24).

Here, we also see that the phase trajectories become closed over time. Let us consider the case for Example 9, when the initial conditions of the system (37) change.

In Figure 25 we see that the phase trajectories reach a stable limit cycle.

In Example 9 the accumulation rate coefficient n was constant. Let us consider an example when the accumulation rate is a function of capital productivity $y\left(t\right)$ according to relation (14).

Example 10.

Let us consider the case when relations (14) and (45) are satisfied. We will take the parameters for modeling from the previous Example 9. The function that implements the numerical algorithm for solving the fractional system (37) taking into account (14) and (45) can be represented in MATLAB as the function ABMDubovskiyFracLineny.m (see Appendix A.2). Let us consider the case where in relation (14) the values $l_1=0.6$, $l_2=1.1$, and the remaining parameters are taken from the previous Example 9. The values of the parameters $M_1$ and $M_2$ in Formula (45) change.

We see in Figure 26 that the limit cycles have a complex shape and are strongly deformed as the values of $M_1$ and $M_2$ increase. An interesting case is when the values of the parameters $l_1$ and $l_2$ change. The type of functions is presented below.

Figure 27 shows oscillograms and phase trajectories for Example 10 with parameter values $l_1=\left(0.6,0.4,0.2\right)$ and $l_2=\left(1.1,1.5,1.9\right)$. Here, we also see the exit of the phase trajectory to the limit cycle.

To conclude this example, we present the case where the initial conditions a and b change (Figure 28).

5. Conclusions

The following conclusions can be drawn from the conducted research. Mathematical models have been developed taking into account constant and variable heredity, which generalize the previously known classical model of S.V. Dubovsky to describe N.D. Kondratiev long waves within the framework of the interaction between the effectiveness of new technologies and the efficiency of capital productivity. These models are based on a system of nonlinear ordinary differential equations with derivatives of fractional constant and variable orders, which are understood in the Gerasimov–Caputo sense with corresponding local initial conditions. In addition, the models took into account the rate of accumulation as a function of capital productivity and external inflow of investment and new technologies. Using the Bendixson criterion, it is shown that the classical model of S.V. Dubovsky can have closed trajectories.

Numerical algorithms have been developed to search for solutions to the proposed mathematical models based on the Adams–Bashforth–Moulton method from the family of predictor-correctors. For such methods, an error analysis was carried out, and the results were formulated in the form of lemmas and theorems. Next, using test examples and Runge’s rule, estimates of the accuracy of calculations were obtained, which gave good agreement with the theoretical results.

Computer experiments were carried out using the proposed numerical algorithms. Oscillograms and phase trajectories were constructed for various values of the model parameters. It is shown that in the absence of external influence, the fractional mathematical model of S.V. Dubovsky oscillations are only damped, i.e., there are no self-oscillations. The presence of external influence makes the existence of a limit cycle possible, which may not always be stable. The limit cycle can be greatly deformed if the orders of fractional derivatives are functions of time, and the shape of the limit cycle can also be influenced by the function of the rate of accumulation from capital productivity. Complex oscillatory regimes deserve separate study and, accordingly, economic interpretation.

It should be noted that further development of the mathematical models proposed in this article may consist of taking into account the dependence of the orders of fractional derivatives not only on time but also on the solution function by analogy with article [78]. Another direction of research is related to the qualitative analysis of dynamic regimes and the construction of their maps by analogy with article [79]. It is also necessary to apply the developed models to the study of experimental data on the crises and cycles of various economies of countries in order to obtain estimates of the periods of economic cycles for their possible forecasting in the future.

Author Contributions

Conceptualization, R.P. and D.M.; methodology, R.P.; software, D.M.; validation, R.P., Z.R. and D.M.; formal analysis, R.P. and Z.R.; investigation, D.M.; resources, D.M.; writing—original draft preparation, D.M. and R.P.; writing—review and editing, D.M., Z.R. and R.P.; visualization, D.M.; supervision, R.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

CMMD	Classical Mathematical Model of S.V. Dubovsky
ABM	Adams–Bashforth–Moulton Method
FMMD	Fractional Mathematical Model of S.V. Dubovsky
FMMDVO	Fractional Mathematical Model of S.V. Dubov Variable Order

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Example of a simply connected region for [Forumla omitted. See PDF.].

Figure 2. Oscillograms and phase trajectories, constructed for various values of [Forumla omitted. See PDF.] and a.

Figure 3. Examples of a simply connected region: (a) for dependence (13); (b) for dependence (14).

Figure 4. Examples of oscillograms and phase trajectories for dependence (13). The points on the oscillogram are the coordinates between the two maximum peaks.

Figure 5. Examples of oscillograms and phase trajectories for dependence (14). The points on the oscillogram are the coordinates between the two maximum peaks.

Figure 6. Oscillograms and phase trajectory of system (16) for dependence (13) taking into account the following parameter values: [Forumla omitted. See PDF.].

Figure 7. Oscillograms and phase trajectory of system (16) for dependence (13) taking into account the following parameter values: [Forumla omitted. See PDF.].

Figure 8. Oscillograms and phase trajectory of system (16) for dependence (14) taking into account the values of the following parameters: [Forumla omitted. See PDF.].

Figure 9. Oscillograms and phase trajectory of system (16) for dependence (14) taking into account the values of the following parameters: [Forumla omitted. See PDF.].

Figure 10. Oscillograms and phase trajectory of system (16) for dependence (14) taking into account the values of the following parameters: [Forumla omitted. See PDF.].

Figure 11. Oscillograms and phase trajectories, constructed for various values of [Forumla omitted. See PDF.] and a.

Figure 12. Oscillograms and phase trajectories, constructed for various values [Forumla omitted. See PDF.].

Figure 13. Oscillograms and phase trajectories, constructed for various values [Forumla omitted. See PDF.].

Figure 14. Oscillogram for [Forumla omitted. See PDF.].

Figure 15. Oscillograms and phase trajectories constructed for different values of the initial conditions [Forumla omitted. See PDF.].

Figure 16. Oscillograms and phase trajectories plotted for different values of [Forumla omitted. See PDF.].

Figure 17. Oscillograms and phase trajectories plotted for different values of [Forumla omitted. See PDF.].

Figure 18. Phase trajectories obtained for different values of the parameters [Forumla omitted. See PDF.] with fixed values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 19. Phase trajectories obtained for different values of the parameters [Forumla omitted. See PDF.] with fixed values [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 20. Phase trajectories obtained for different values of the parameters [Forumla omitted. See PDF.] with fixed values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for dependency (14).

Figure 21. Simulation results in Example 9 for different values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 22. Simulation results in Example 9 for different values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 23. Simulation results in Example 9 for various values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 24. Simulation results in Example 9 for various values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 25. Simulation results in Example 9 for various values of a and b.

Figure 26. Simulation results in Example 10 for various values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 27. Simulation results in Example 10 for various values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 28. Simulation results in Example 10 for various values of a and b.

Table 1

Case of CMMD.

N	$\mathit{\tau }$	${\mathit{\xi }}_{\mathit{x}}$	${\mathit{\xi }}_{\mathit{y}}$	${\mathit{p}}_{\mathit{x}}$	${\mathit{p}}_{\mathit{y}}$
10	1/10	0.0397415350	0.0799740787	–	–
20	1/20	0.0060318953	0.0187709253	2.72	2.09
40	1/40	0.0010575820	0.0045870660	2.51	2.03
80	1/80	0.0002336673	0.0011345727	2.18	2.02
160	1/160	0.0000552493	0.0002821520	2.08	2.01
320	1/320	0.000013453	0.0000703530	2.04	2.00

Table 2

Case of FMMD.

N	$\mathit{\tau }$	${\mathit{\xi }}_{\mathit{x}}$	${\mathit{\xi }}_{\mathit{y}}$	${\mathit{p}}_{\mathit{x}}$	${\mathit{p}}_{\mathit{y}}$
10	1/10	0.1171054420	0.2158809510	–	–
20	1/20	0.0235086005	0.0586392164	2.31	1.88
40	1/40	0.0041498810	0.0173332737	2.50	1.75
80	1/80	0.0009959756	0.0051032771	2.05	1.76
160	1/160	0.0002742665	0.0014889669	1.86	1.77
320	1/320	0.0000779949	0.0004308895	1.81	1.79

Table 3

Case of VO FMMD.

N	$\mathit{\tau }$	${\mathit{\xi }}_{\mathit{x}}$	${\mathit{\xi }}_{\mathit{y}}$	${\mathit{p}}_{\mathit{x}}$	${\mathit{p}}_{\mathit{y}}$
10	1/10	0.4740440$·{10}^{-3}$	0.1036998$·{10}^{-2}$	–	–
20	1/20	0.9818199$·{10}^{-4}$	0.4061758$·{10}^{-3}$	1.57	1.41
40	1/40	0.3289388$·{10}^{-4}$	0.1522497$·{10}^{-3}$	1.54	1.45

Table 4

Case of VO FMMD.

N	$\mathit{\tau }$	${\mathit{\xi }}_{\mathit{x}}$	${\mathit{\xi }}_{\mathit{y}}$	${\mathit{p}}_{\mathit{x}}$	${\mathit{p}}_{\mathit{y}}$
10	1/10	0.5457955$·{10}^{-3}$	0.4676925$·{10}^{-3}$	–	–
20	1/20	0.2583299$·{10}^{-3}$	0.2446991$·{10}^{-3}$	1.18	1.11
40	1/40	0.1137234$·{10}^{-3}$	0.1127146$·{10}^{-3}$	1.26	1.19
80	1/80	0.4747948$·{10}^{-4}$	0.4909588$·{10}^{-4}$	1.29	1.22
160	1/160	0.1937985$·{10}^{-4}$	0.2094110$·{10}^{-4}$	1.30	1.29

Appendix A. Listing of FDDSVO Program Procedures

Appendix A.1. ABMDubovskiyFracLine.m

Description; [x,y]=ABMDubovskiyFracLine(T,tau,alpha,delta,phi,par,Init) T - simulation time; tau - computational mesh discretization step; alpha - [alpha10 alpha20 M1 M2] parameters for formulas (50); delta - [delta1 delta 2] amplitude values on the right sides of the system (42); phi - [phi1 phi 2] frequency values on the right sides of the system (42); par - [n lambda] savings ratio values and lambda; Init - [a b], defining the initial conditions of the problem (42). The program code that implements the function is given below: function [x,y]=ABMDubovskiyFracLine(T,tau,alpha,delta,phi,par,Init) alpha10=alpha(1); alpha20=alpha(2);M1=alpha(3);M2=alpha(4); delta1=delta(1); delta2=delta(2); phi1=phi(1); phi2=phi(2); n=par(1); lambda=par(2); N=round(T/tau); a=zeros(N,N); b=zeros(N,N); x=zeros(1,N); y=zeros(1,N); xp=zeros(1,N); yp=zeros(1,N); x1=1.3; y1=0.5; xp(1)=Init(1); yp(1)=Init(2); x(1)=xp(1); y(1)=yp(1);alpha1=zeros(1,N);alpha2=zeros(1,N); for i=0:N alpha1(i+1)=alpha10-M1.*tau.*i/T; alpha2(i+1)=alpha20-M2.*tau.*i/T; a(1,i+1)=i.^(alpha1(i+1)+1)-(i-alpha1(i+1))*(i+1).^(alpha1(i+1)); b(1,i+1)=i.^(alpha2(i+1)+1)-(i-alpha2(i+1))*(i+1).^(alpha2(i+1)); end for i=1:N a(i+1,i+1)=1; b(i+1,i+1)=1; for j=2:i a(j,i+1)=(i-j+2)^(alpha1(i)+1)+(i-j)^(alpha1(i)+1)-2*(i-j+1)^(alpha1(i)+1); b(j,i+1)=(i-j+2)^(alpha2(i)+1)+(i-j)^(alpha2(i)+1)-2*(i-j+1)^(alpha2(i)+1); end end for k=1:N-1 p=1:k; xp(k+1)=xp(1)+tau^(alpha1(k))/gamma(alpha1(k)+1)*... sum(((k-p+1).^alpha1(k)-(k-p).^alpha1(k)).*... (delta1*cos(phi1*p*tau)-lambda*n*xp(p).*(xp(p)-1).*(yp(p)-y1))); yp(k+1)=yp(1)+tau^(alpha2(k))/gamma(alpha2(k)+1)*... sum(((k-p+1).^alpha2(k)-(k-p).^alpha2(k)).*... (delta2*cos(phi2*p*tau)+n*(1-n)*yp(p).^2.*(xp(p)-x1))); x(k+1)=x(1)+tau^(alpha1(k))/gamma(alpha1(k)+2)*... (delta1*cos(phi1*(k+1)*tau)-lambda*n*xp(k+1).*(xp(k+1)-1).*... (yp(k+1)-y1)+sum((delta1*cos(phi1*p*tau)-... lambda*n*x(p).*(x(p)-1).*(y(p)-y1))*a(p,k+1))); y(k+1)=y(1)+tau^(alpha2(k))/gamma(alpha2(k)+2)*... (delta2*cos(phi2*(k+1)*tau)+n*(1-n)*yp(k+1).^2.*(xp(k+1)-x1)+... sum((delta2*cos(phi2*p*tau)+n*(1-n)*y(p).^2.*(x(p)-x1))*b(p,k+1))); end end Example of performing a procedure and visualizing the results. [x1,y1]=ABMDubovskiyFracLine(500,0.04,[0.8 0.7 0.01 0.03],[0 0.5],[0 1.5], [0.2 2.25],[1.35 0.5]); [x2,y2]=ABMDubovskiyFracLine(500,0.04,[0.6 0.5 0.01 0.03],[0 0.5],[0 1.5], [0.2 2.25],[1.35 0.5]); [x3,y3]=ABMDubovskiyFracLine(500,0.04,[0.4 0.3 0.01 0.03],[0 0.5],[0 1.5], [0.2 2.25],[1.35 0.5]); t=0.01:0.04:500; figure(1) subplot(2,2,1); plot(t,x1,’k’,t,x2,’b’,t,x3,’r’) xlabel(’t’); ylabel(’x’); title(’Oscillogram x(t)’) subplot(2,2,2); plot(t,y1,’k’,t,y2,’b’,t,y3,’r’) xlabel(’t’); ylabel(’y’); title(’Oscillogram y(t)’) subplot(2,2,[3,4]); plot(x1,y1,’k’, x2,y2,’b’,x3,y3,’r’) xlabel(’x’); ylabel(’y’); title(’Phase trajectory’)

Appendix A.2. ABMDubovskiyFracLineny.m

Description; [x,y]=ABMDubovskiyFracLineny(T,tau,alpha,delta,phi,par,Init) T - simulation time; tau - computational mesh discretization step; alpha - [alpha10 alpha20 M1 M2] parameters for formulas (50); delta - [delta1 delta 2] amplitude values on the right sides of the system (42); phi - [phi1 phi 2] frequency values on the right sides of the system (42); par - [l1 l2 lambda] values l1 and l2 from (14), as well as the accumulation coefficient and lambda; Init - [a b], defining the initial conditions of the problem (42). The program code that implements the function is given below: function [x,y]=ABMDubovskiyFracLineny(T,tau,alpha,delta,phi,par,Init) N=round(T/tau); alpha10=alpha(1); alpha20=alpha(2);M1=alpha(3); M2=alpha(4); delta1=delta(1);delta2=delta(2); phi1=phi(1);phi2=phi(2); l1=par(1); l2=par(2);lambda=par(3);tau=T/N; a=zeros(N,N);b=zeros(N,N);xp=zeros(1,N); yp=zeros(1,N);x=zeros(1,N);y=zeros(1,N);alpha1=zeros(1,N);alpha2=zeros(1,N); x1=1.3; y1=0.5; xp(1)=Init(1); yp(1)=Init(2); x(1)=xp(1); y(1)=yp(1); for i=0:N alpha1(i+1)=alpha10-M1.*tau.*i/T; alpha2(i+1)=alpha20-M2.*tau.*i/T; a(1,i+1)=i^(alpha1(i+1)+1)-(i-alpha1(i+1))*(i+1)^(alpha1(i+1)); b(1,i+1)=i^(alpha2(i+1)+1)-(i-alpha2(i+1))*(i+1)^(alpha2(i+1)); end for n=1:N a(i+1,i+1)=1; b(i+1,i+1)=1; for j=2:i a(j,i+1)=(i-j+2)^(alpha1(i)+1)+(i-j)^(alpha1(i)+1)-2*(i-j+1)^(alpha1(i)+1); b(j,i+1)=(i-j+2)^(alpha2(i)+1)+(i-j)^(alpha2(i)+1)-2*(i-j+1)^(alpha2(i)+1); end end for k=1:N-1 p=1:k; xp(k+1)=xp(1)+tau^(alpha1(k))/gamma(alpha1(k)+1)*... sum(((k-p+1).^alpha1(k)-(k-p).^alpha1(k)).*(delta1*cos(phi1*p*tau)-... lambda*(l1*(l2-1)./yp(p)).*xp(p).*(xp(p)-1).*(yp(p)-y1))); yp(k+1)=yp(1)+tau^(alpha2(k))/gamma(alpha2(k)+1)*... sum(((k-p+1).^alpha2(k)-(k-p).^alpha2(k)).*(delta2*cos(phi2*p*tau)+... (l1*(l2-1)./yp(p)).*(1-(l1*(l2-1)./yp(p))).*yp(p).^2.*(xp(p)-x1))); x(k+1)=x(1)+tau^(alpha1(k))/gamma(alpha1(k)+2)*... (delta1*cos(phi1*(k+1)*tau)-lambda*(l1*(l2-1)./yp(k+1)).*xp(k+1).*... 5(xp(k+1)-1).*(yp(k+1)-y1)+... sum((delta1*cos(phi1*p*tau)-lambda*(l1*(l2-1)./y(p)).*x(p).*(x(p)-1).*... (y(p)-y1))*a(p,k+1))); y(k+1)=y(1)+tau^(alpha2(k))/gamma(alpha2(k)+2)*... (delta2*cos(phi2*(k+1)*tau)+(l1*(l2-1)./yp(k+1)).*(1-(l1*(l2-1)./yp(k+1))).*... yp(k+1).^2.*(xp(k+1)-x1)+... sum((delta2*cos(phi2*p*tau)+(l1*(l2-1)./y(p)).*(1-(l1*(l2-1)./y(p))).*... y(p).^2.*(x(p)-x1))*b(p,k+1))); end end

References

1. Clarke, H. Physical Economy. A Preliminary Inquiry into the Physical Laws Governing the Periods of Famines and Panics. Railw. Regist.; 1847; 5, pp. 155-169.

2. Juglar, C. Des Crises Commercials et de leur Retour Periodique en France, en Angleterre et aux EU; Guillaumin: Paris, France, 1862; 258p.

3. Jevons, W.S. The solar-commercial cycle. Nature; 1882; 26, pp. 226-228. [DOI: https://dx.doi.org/10.1038/026226a0]

4. Tugan-Baranovskiy, M. Industrial Crises. Essay on the Social History of England; O. N. Popova: Saint Petersburg, Russia, 1900; 354p.

5. Kitchin, J. Cycles and Trends in Economic Factors. Rev. Econ. Stat.; 1923; 5, pp. 10-16. [DOI: https://dx.doi.org/10.2307/1927031]

6. Kondratyev, N.D.; Oparin, D.N. Large Cycles of Market Conditions; Institute of Economics: Moscow, Russia, 1928; 287p.

7. Solow, R.M. Technical Change and the Aggregate Production Function. Rev. Econ. Stat.; 1957; 39, pp. 312-320. [DOI: https://dx.doi.org/10.2307/1926047]

8. Rostow, W.W. Technical The Stages of Economic Growth. Econ. Hist. Rev. New Ser.; 1959; 12, pp. 1-16. [DOI: https://dx.doi.org/10.1111/j.1468-0289.1959.tb01829.x]

9. Mensch, G. Stalemate in Technology. Innovations Overcome the Depression; Ballinger Pub Co.: Cambridge, UK, 1979; 241p.

10. Gattei, G. Every 25 years? Strike waves and long economic cycles. Proceedings of the Colloquium “The Long Waves in the Economic Conjuncture—The Present State of the International Debates”; Brussels, Belgium, 12–14 January 1989; pp. 12-14.

11. Menshikov, S.M.; Klimenko, L.A. Long Waves in Economics; International Relations: Moscow, Russia, 1989; pp. 145-149.

12. Amos, O.M., Jr. Growth pole cycles: A synthesis of growth pole and long wave theories. Rev. Reg. Stud.; 1990; 20, pp. 37-48.

13. Tylecotte, A. Long Wave in the World Economy: The Present Crisis in Historical Perspective; Psychology Press: London, UK, New York, NY, USA, 1993; 338p.

14. Dubovsky, S.V. Modeling object Kondratiev cycle. Math. Model.; 1995; 7, pp. 65-74.

15. Alexander, M.A. The Kondratiev Cycle: A Generational Interpretation; IUniverse: Bloomington, IN, USA, 2002; 314p.

16. Kondratiev, N.D.; Yakovets, Y.V.; Abalkin, L.I. Large Cycles of Conjuncture and the Theory of Foresight: Selected Works; Economics: Moscow, Russia, 2002; 766p.

17. Hirooka, M. Innovation Dynamism and Economic Growth. A Nonlinear Perspective; Edward Elgar: Cheltenham, Northampton, UK, 2006.

18. Akaev, A.A. Derivation of the general macroeconomic dynamics equation describing the joint interaction of long-term growth and business cycles. Dokl. Math.; 2007; 76, pp. 879-881. [DOI: https://dx.doi.org/10.1134/S1064562407060191]

19. Akaev, A.A. Influence of business cycles on long-term economic growth. Dokl. Math.; 2008; 78, pp. 621-625. [DOI: https://dx.doi.org/10.1134/S106456240804039X]

20. Muller, K.H. Farewell to Long Waves: Substituting Cyclical Approaches in Innovation and Technology Research with a RISC-Framework; University of Ljubljana: Ljubljana, Slovenia, 2008; 54p.

21. Korotaev, A.V.; Tsirel, S.V. Kondratieff Waves in Global Economic Dynamics. System Monitoring. Global and Regional Development; Librocom: Moscow, Russia, 2009; 229p.

22. Dementyev, V.E. Long Waves of Economic Development and Financial Bubbles; CEMI RAS: Moscow, Russia, 2009.

23. Akaev, A.A. Qualitative analysis of business cycles’ influence on economic growth. Econ. Math. Methods; 2009; 45, pp. 78-137.

24. Gladkikh, I.P. Long Waves in the Post-Industrial Economy: Theoretical Foundations and Features; Lambert Academic: Saarbrücken, Germany, 2012.

25. Dubovsky, S.V. Modeling Kondratiev cycles and forecasting crises. Kondratiev Waves. Aspects and Prospects; Akaev, A.A. Uchitel’: Volgograd, Russia, 2012; pp. 179-188.

26. Akaev, A.A. Mathematical foundations of the Schumpeter-Kondratiev innovation-cyclical theory of economic development. Kondratiev Waves. Aspects and Prospects; Akaev, A.A. Uchitel’: Volgograd, Russia, 2012; pp. 314-341.

27. Makarov, D.V.; Parovik, R.I. Modeling of the economic cycles using the theory of fractional calculus. J. Internet Bank. Commer.; 2016; 21, S21.

28. Kleinknecht, A. Innovation Patterns in Crisis and Prosperity: Schumpeter’s Long Cycle Reconsidered; Springer: Germany, 2016; 235p.

29. Kondratyev, N.D. Large Cycles of Market Conditions. Selected Works; Yurayt Publishing House: Moscow, Russia, 2021; 490p.

30. Makarov, D.; Parovik, R. Numerical modeling of Kondratyev’s long waves taking into account heredity. Aip Conf. Proc.; 2021; 2365, 020007.

31. Makarov, D.V.; Parovik, R.I. A computer program for the numerical analysis of economic cycles within the framework of the Dubovsky generalized model. Aip Conf. Proc.; 2022; 2467, 060015.

32. Makarov, D.V.; Parovik, R.I. Fractional mathematical model S.V. Dubovsky and economic cycles. Probl. Comput. Appl. Math.; 2023; 52, pp. 8-25.

33. Makarov, D.V.; Parovik, R.I. Fractional mathematical model S.V. Dubovsky with the effect of variable heredity. Probl. Comput. Appl. Math.; 2023; 53, pp. 5-22.

34. Makarov, D.V. The classical mathematical model of S.V. Dubovsky and some of its modifications for describing K-waves. Vestn. Kraunc.-Fiz.-Mat. Nauk.; 2024; 46, pp. 52-69.

35. Landa, P.S. Regular and Chaotic Oscillations; Springer: Berlin, Germany, 2012; 397p.

36. Lichtenberg, A.J.; Lieberman, M.A. Regular and Chaotic Dynamics; Springer: Berlin, Germany, 2013; 692p.

37. Oldham, K.B.; Spanier, J. The Fractional Calculus. Theory and Applications of diffErentiation and Integration to Arbitrary Order; Academic Press: London, UK, 1974; 240p.

38. Miller, K.S.; Ross, B. An Introduction to the Fractional Calculus and Fractional diffErential Equations; A Wiley-Interscience publication: New York, NY, USA, 1993; 384p.

39. Nakhushev, A.M. Fractional Calculus and its Application; Fizmatlit: Moscow, Russia, 2003; 272p.

40. Pskhu, A.V. Fractional Partial Differential Equations; Nauka: Moscow, Russia, 2005; 199p.

41. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006; 523p.

42. Tarasov, V.E. On history of mathematical economics: Application of fractional calculus. Mathematics; 2019; 7, 509. [DOI: https://dx.doi.org/10.3390/math7060509]

43. Granger, C.W.J. The typical spectral shape of an economic variable. Econometrica; 1966; 34, pp. 150-161. [DOI: https://dx.doi.org/10.2307/1909859]

44. Granger, C.W.J.; Joyeux, R. An introduction to long memory time series models and fractional differencing. J. Time Ser. Anal.; 1980; 1, pp. 15-39. [DOI: https://dx.doi.org/10.1111/j.1467-9892.1980.tb00297.x]

45. Hosking, J.R.M. Equivalent Forms of the Multivariate Portmanteau Statistic. J. R. Stat. Soc.; 1981; B 43, pp. 261-262. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1981.tb01179.x]

46. Gray, H.L.; Zhang, N.; Woodward, W.A. On generalized fractional processes. J. Time Ser. Anal.; 1989; 10, pp. 233-257. Erratum in J. Time Ser. Anal. 1994, 15, 561–562 [DOI: https://dx.doi.org/10.1111/j.1467-9892.1989.tb00026.x]

47. Gil-Alana, L.A. Testing stochastic cycles in macroeconomic time series. J. Time Ser. Anal.; 2001; 22, pp. 411-430. [DOI: https://dx.doi.org/10.1111/1467-9892.00233]

48. Rogers, L.C.G. Arbitrage with fractional Brownian motion. Math. Financ.; 1997; 7, pp. 95-105. [DOI: https://dx.doi.org/10.1111/1467-9965.00025]

49. Scalas, E.; Gorenflo, R.; Mainardi, F. Fractional calculus and continuous-time finance. Phys. A Stat. Mech. Its Appl.; 2000; 284, pp. 376-384. [DOI: https://dx.doi.org/10.1016/S0378-4371(00)00255-7]

50. Mainardi, F.; Raberto, M.; Gorenflo, R.; Scalas, E. Fractional calculus and continuous-time finance II: The waiting-time distribution. Phys. A Stat. Mech. Appl.; 2000; 287, pp. 468-481. [DOI: https://dx.doi.org/10.1016/S0378-4371(00)00386-1]

51. Chen, W.C. Nonlinear dynamics and chaos in a fractional-order financial system. Chaos Solitons Fract.; 2008; 36, pp. 1305-1314. [DOI: https://dx.doi.org/10.1016/j.chaos.2006.07.051]

52. Tarasov, V.E.; Tarasova, V.V. Criterion of existence of power-law memory for economic processes. Entropy; 2018; 20, 414. [DOI: https://dx.doi.org/10.3390/e20060414] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33265504]

53. Pakhira, R.; Ghosh, U.; Sarkar, S. Study of memory effect in an inventory model with quadratic type demand rate and salvage value. Appl. Math. Sci.; 2019; 13, pp. 209-223. [DOI: https://dx.doi.org/10.12988/ams.2019.9111]

54. Bendixson, I. Sur les courbes définies par des équations différentielles. Acta Math.; 1901; 24, pp. 1-88. [DOI: https://dx.doi.org/10.1007/BF02403068]

55. Iserles, A. A First Course in the Numerical Analysis of Differential Equations; Cambridge University Press: Cambridge, UK, 2009; 459p.

56. Stoleru, L. Economic Equilibrium and Growth; North-Holland Publishing Company: Amsterdam, The Netherlands, 1975; 300p.

57. Volterra, V. Functional Theory, Integral and Integro-Differential Equations; Dover Publications: New York, NY, USA, 2005; 288p.

58. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: Amsterdam, The Netherlands, 1998; 340p.

59. Mainardi, F. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models; World Scientific: Singapore, 2010; 368p.

60. Gerasimov, A.N. Generalization of the laws of linear deformation and their application to problems of internal friction. Acad. Sci. Ssr. Appl. Math. Mech.; 1948; 44, pp. 62-78.

61. Caputo, M. Linear models of dissipation whose Q is almost frequency independent—II. Geophys. J. Int.; 1967; 13, pp. 529-539. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1967.tb02303.x]

62. Rabotnov, Y.N. Elements of Hereditary Mechanics of Solids; MIR Publishers: Moscow, Russia, 1980; 387p.

63. Pskhu, A.V.; Rekhviashvili, S.S. Analysis of Forced Oscillations of a Fractional Oscillator. Tech. Phys. Lett.; 2018; 44, pp. 1218-1221. [DOI: https://dx.doi.org/10.1134/S1063785019010164]

64. Aguilar, J.F.G.; Hernández, M.M. Space-time fractional diffusion-advection equation with Caputo derivative. Abstr. Appl. Anal.; 2014; 2014, 283019. [DOI: https://dx.doi.org/10.1155/2014/283019]

65. Butenkov, S.A. Mathematical models of processes on fractal structures with given properties based on granulation methods. News South. Fed. Univ. Tech. Sci.; 2011; 121, pp. 199-209.

66. Diethelm, K.; Ford, N.J.; Freed, A.D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn.; 2002; 29, pp. 3-22. [DOI: https://dx.doi.org/10.1023/A:1016592219341]

67. Yang, C.; Liu, F. A computationally effective predictor-corrector method for simulating fractional order dynamical control system. ANZIAM J.; 2005; 47, pp. 168-184. [DOI: https://dx.doi.org/10.21914/anziamj.v47i0.1037]

68. Garrappa, R. Numerical solution of fractional differential equations: A survey and a software tutorial. Mathematics; 2018; 6, 016. [DOI: https://dx.doi.org/10.3390/math6020016]

69. Gavrilyuk, I. Exact and Truncated Difference Schemes for Boundary Value ODEs; Springer Science & Business Media: Basel, Switzerland, 2011; 247p.

70. Patnaik, S.; Hollkamp, J.P.; Semperlotti, F. Applications of variable-order fractional operators: A review. Proc. R. Soc. A; 2020; 476, 20190498. [DOI: https://dx.doi.org/10.1098/rspa.2019.0498] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32201475]

71. Ma, S.; Xu, Y.; Yue, W. Numerical Solutions of a Variable-Order Fractional Financial System. J. Appl. Math.; 2012; 2012, 417942. [DOI: https://dx.doi.org/10.1155/2012/417942]

72. Samko, S.G.; Ross, B. Integration and differentiation of a variable fractional order. Integral Transform. Spec. Funct.; 1993; 1, pp. 277-300. [DOI: https://dx.doi.org/10.1080/10652469308819027]

73. Samko, S.G. Fractional integration and differentiation of variable order. Anal. Math.; 1995; 21, pp. 213-236. [DOI: https://dx.doi.org/10.1007/BF01911126]

74. Sun, H.; Chang, A.; Zhang, Y.; Chen, W. A review on variable-order fractional differential equations: Mathematical foundations, physical models, numerical methods and applications. Fract. Calc. Appl. Anal.; 2019; 22, pp. 27-59. [DOI: https://dx.doi.org/10.1515/fca-2019-0003]

75. Zhou, Y.; Suzuki, J.L.; Zhang, C.; Zayernouri, M. Implicit-explicit time integration of nonlinear fractional differential equations. Appl. Numer. Math.; 2020; 156, pp. 555-583. [DOI: https://dx.doi.org/10.1016/j.apnum.2020.04.006]

76. Parovik, R.I. On a Finite-Difference Scheme for an Hereditary Oscillatory Equation. J. Math. Sci.; 2021; 253, pp. 547-557. [DOI: https://dx.doi.org/10.1007/s10958-021-05252-2]

77. Yousuf, M.; Sarwar, S. Highly Efficient Numerical Algorithm for Nonlinear Space Variable-Order Fractional Reaction–Diffusion Models. Fractal Fract.; 2023; 7, 688. [DOI: https://dx.doi.org/10.3390/fractalfract7090688]

78. Kanwal, A.; Boulaaras, S.; Shafqat, R.; Taufeeq, B.; Ur Rahman, M. Explicit scheme for solving variable-order time-fractional initial boundary value problems. Sci. Rep; 2024; 14, 5396. [DOI: https://dx.doi.org/10.1038/s41598-024-55943-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38443513]

79. Parovik, R.I.; Yakovleva, T.P. Construction of maps for dynamic modes and bifurcation diagrams in nonlinear dynamics using the Maple computer mathematics software package. J. Phys. Conf. Ser.; 2022; 2373, 52022. [DOI: https://dx.doi.org/10.1088/1742-6596/2373/5/052022]