Fractional Mathematical Models S.V. Dubovsky to Describe K-Waves in Economics

1. Introduction

Numerous research results show that long wavelength N.D. Kondratiev (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,35,36,37,38]. K-waves can be described using dynamic (oscillatory) systems, which mainly characterize nonlinear oscillatory processes [11,14,18,25,26,37]. Such processes are studied to describe the chaotic and regular dynamics of nonlinear oscillatory systems. In world 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. [39,40]. 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 [41,42,43,44,45].

In the review work of Tarasov V.E. [46] 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 [47,48], 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 [47,48]; fractional Brownian motion [49]; econophysics [50,51]; deterministic chaos [52]; mathematical economics [53,54]. 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,21,22,31,34]. However, the analysis of the literature on the research topic showed that the application of fractional calculus to the dynamic system of S.V. Dubovsky to describe K-waves has not yet been carried out. 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 external periodic influences - investments in the efficiency of fixed assets and the influx of 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 research in this article.

The research plan in the article is organized as follows. The first section provides an overview of the current state of research in 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 second section examines the classical model of S.V. Dubovsky: the basic equations of the model are derived, two modifications of the model are also given - taking into account the rate of accumulation depending on capital productivity and external inflow of investments and new technologies, it is shown that the model can have closed phase trajectories, oscillograms are constructed using the Adams-Bashforth-Moulton numerical algorithm and phase trajectories, interpretation of research results is given. The third section provides a generalization of the classical model of S.V. Dubovsky in the case of taking into account constant memory: the definition of the fractional derivative in the Gerasimov-Caputo sense 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 fourth section provides a generalization of the fractional model of S.V. Dubovsky in the case of taking variable memory into account: 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 different values of the model parameters, an interpretation of the research results is given. The fifth section provides conclusions based on the research results, 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. was studied. Dubovsky to describe long waves N.D. Kondratiev (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 in 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

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

(1)

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

(2)

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

(3)

where $Y\left(t\right)$ is 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, $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), more interesting from the point of view of revealing the fundamentals of the production function, 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:

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

(4)

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)=\frac{Y\left(t\right)}{K\left(t\right)}$; efficiency of new technologies: $x\left(t\right)=\frac{u\left(t\right)}{U\left(t\right)}$.

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

(5)

Based on 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 formally wrote this as the following equation:

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

(6)

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):

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

(7)

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^{*},y\left(t\right)=y^{*}$, which agrees well with real statistics.

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

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

(8)

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

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

(9)

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

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

(10)

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, $\lambda >0$ — constant coefficient, which is determined from statistics, $x^{*}$ and $y^{*}$ 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 equations (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 do 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^{*}\right)$ and $f_2\left(x,y\right)=n\left(1-n\right)y^2\left(t\right)\left(x\left(t\right)-x^{*}\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 $\Omega \subset R^2$. Then the following theorem is true.

Theorem 1.(Bendixson’s theorem [55]) 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 equations (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.

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

Then

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

(11)

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 $\Omega \subset R^2$ can occur. Indeed, if you choose

$$0.5<x<x^{*},0<y<\frac{\lambda y^{*}}{\lambda -1+n}.$$

(12)

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 the sign changes for the parameter values are performed (Figure 1): $n=0.2,\lambda =2.25,x^{*}=1.3,y^{*}=0.5$.

Based on the above, we can conclude that 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’s 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 2nd order of convergence (trapezoidal rule) [56], 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$) lead 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 system (10) $\left(x^{*},y^{*}\right)$ is the center) . The periods of phase trajectories according to [14]: T = 50.1, 54.9, 60.6 years, and their direction of travel is counterclockwise.

The interpretation of the economic cycle is given according to the phase portrait (Figure 2c), 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 — 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:

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

(13)

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)$ [57]. However, dependence (13) in his work [14] S.V. Dubovsky did not investigate.

Remark 5.

In article [14], dependence (13) was studied, oscillograms and phase trajectories were obtained for various parameter values, but their economic interpretation was not given.

Remark 6.

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

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

(14)

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

However, in [53] there was no study of dependence (14). Therefore, at this point we will eliminate this drawback — we will conduct a study of the classical model of S.V. Dubovsky (10) taking into account dependencies (13) and (14), and also give their economic interpretation.

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)$:

$$\left\{\begin{array}{cc}\hfill & \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^{*}\right),\hfill \\ \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^{*}\right),\hfill \\ \hfill & x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$$

(15)

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, in contrast to Figure 4c, the accumulation rate enhances 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 number of years. Then, for example, we see that for the 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

$$\left\{\begin{array}{cc}\hfill & \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^{*}\right)+f_1\left(t\right),\hfill \\ \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^{*}\right)+f_2\left(t\right),\hfill \\ \hfill & x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$$

(16)

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).

In 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 on the oscillograms an increase in the amplitude of natural oscillations. 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 [32,33,34,35,36,37,38], a generalization of the CMMD was proposed taking into account the heredity effect [58]. 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 [41,42,43,44,45].

It is shown in [35] 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’s 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 [41,42,43,44,45].

Definition 3.

Euler’s gamma function is defined:

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

(17)

Definition 4.

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

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

(18)

where $q>0$, $\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:

$$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).$$

(19)

Definition 5.

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

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

(20)

where $\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:

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

(21)

Remark 12.

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

$$D_{-\infty }^{\alpha }u\left(x,t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}\underset{-\infty }{\overset{t}{\int }}\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’s consider a generalization of the CMMD (16), which according to (20) can be represented as the following Cauchy problem:

$$\left\{\begin{array}{cc}\hfill & {\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^{*}\right)+f_1\left(t\right),0<\alpha <1,\hfill \\ \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^{*}\right)+f_2\left(t\right),0<\beta <1,\hfill \\ \hfill & x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$$

(22)

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) [64]. For example, if the Cauchy problem is written in general form:

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

(23)

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 }$, then 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 system (16) with derivatives of fractional orders [65]:

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

(24)

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 [58]. For example, we can write the original Cauchy problem in the integro-differential formulation:

$$\left\{\begin{array}{cc}\hfill & \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^{*}\right)+f_1\left(t\right),\hfill \\ \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^{*}\right)+f_2\left(t\right),\hfill \\ \hfill & x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$$

(25)

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:

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

(26)

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

Remark 17.

Note that in [66] 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 get 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), to solve it 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 [67,68,69]. To do 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 get the calculation formulas:

$$\left\{\begin{array}{cc}\hfill & x_{k+1}^p=x_0+\frac{{\tau }_1^{\alpha }}{\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^{*}\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right),\hfill \\ \hfill & y_{k+1}^p=y_0+\frac{{\tau }_2^{\alpha }}{\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^{*}\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right),\hfill \\ \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.$$

(27)

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

$$\left\{\begin{array}{cc}\hfill & x_{k+1}=x_0+\frac{{\tau }_1^{\alpha }}{\Gamma \left({\alpha }_1+2\right)}\left(\begin{array}{cc}\hfill & \left(-\lambda n_{k+1}{x}_{k+1}^p\left(x_{k+1}^p-1\right)\left(y_{k+1}^p-y^{*}\right)+{\delta }_{1}cos\left({\omega }_1\tau \left(k+1\right)\right)\right)+\hfill \\ \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^{*}\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right)\hfill \end{array}\right),\hfill \\ \hfill & y_{k+1}=y_0+\frac{{\tau }_2^{\alpha }}{\Gamma \left({\alpha }_2+2\right)}\left(\begin{array}{cc}\hfill & \left(n_{k+1}\left(1-n_{k+1}\right){\left(y_{k+1}^p\right)}^2\left(x_{k+1}^p-x^{*}\right)+{\delta }_{2}cos\left({\omega }_2\tau \left(k+1\right)\right)\right)+\hfill \\ \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^{*}\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right)\hfill \end{array}\right),\hfill \end{array}\right.$$

(28)

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

$${\rho }_{j,k+1}^i=\left\{\begin{array}{cc}& 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 & \hfill 1,,j=k+1,\\ & 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 provides a detailed analysis of the error of the modified ABM method (27), (28). The following auxiliary lemma holds for the ABM method.

Lemma 1.

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

$$\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 ,$$

(29)

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

Proof. 

When developing calculation formulas (27), (28), we used the trapezoidal quadrature formula. Let’s create an error in this representation.

$$\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)=\sum _{j=0}^k\underset{jh}{\overset{\left(j+1\right)h}{\int }}\left(t_{k+1}-t\right)\left(z\left(t\right)-z\left(t_j\right)\right)dt.$$

(30)

Next, we use the well-known mean value theorem (Lagrange’s theorem) [70]:

$$\left|f^{\prime }\left(c\right)\left(b-a\right)\right|\le {∥f^{\prime }∥}_{\infty }\left|a-b\right|,$$

where f is a continuous function on the interval $\left[a,b\right]$, and $c\in \left(a,b\right)$, ${∥ ∥}_{\infty }$ — Chebyshev norm. We get the following estimate:

$$\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 {∥z^{\prime }∥}_{\infty }\sum _{j=0}^k\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}\left(t-t_j\right)dt.$$

(31)

Let us calculate the integral on the right side of (31).

$$\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}\left(t-t_j\right)dt=\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}\left[\left(t-t_{k+1}\right)+\left(t_{k+1}-t_j\right)\right]dt=$$

$$=\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}\left(t-t_{k+1}\right)dt+\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}\left(t_{k+1}-t_j\right)dt=$$

$$=-\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha }dt+\left(t_{k+1}-t\right)\underset{t_j}{\overset{t_{j+1}}{\int }}{\left(t_{k+1}-t_j\right)}^{\alpha -1}dt=$$

$$={\left.\frac{{\left(t_{k+1}-t\right)}^{\alpha +1}}{\alpha +1}\right|}_{t_j}^{t_{j+1}}-\left(t_{k+1}-t_j\right){\left.\frac{{\left(t_{k+1}-t\right)}^{\alpha }}{\alpha }\right|}_{t_j}^{t_{j+1}}=$$

$$=\frac{1}{\alpha +1}\left[{\left(t_{k+1}-t_{j+1}\right)}^{\alpha +1}-{\left(t_{k+1}-t_j\right)}^{\alpha +1}\right]-\frac{1}{\alpha }\left[{\left(t_{k+1}-t_{j+1}\right)}^{\alpha }-{\left(t_{k+1}-t_j\right)}^{\alpha }\right]=$$

$$=\frac{1}{\alpha \left(\alpha +1\right)}\left[\alpha {\left(t_{k+1}-t_{j+1}\right)}^{\alpha +1}-\alpha {\left(t_{k+1}-t_j\right)}^{\alpha +1}-\left(\alpha +1\right){\left(t_{k+1}-t_{j+1}\right)}^{\alpha }-\left(\alpha +1\right)\left(t_{k+1}-t_j\right)\right]=$$

$$=\frac{1}{\alpha \left(\alpha +1\right)}\left[{\left(t_{k+1}-t_{j+1}\right)}^{\alpha +1}\left[\alpha \left(t_{k+1}-t_j\right)-\left(\alpha +1\right)\left(t_{k+1}-t_{j+1}\right)+{\left(t_{k+1}-t_j\right)}^{\alpha +1}\right]\right]=$$

$$=\frac{1}{\alpha \left(\alpha +1\right)}\left[{\left(t_{k+1}-t_{j+1}\right)}^{\alpha +1}\left[\alpha \left(t_{k+1}-t_j\right)-\left(\alpha +1\right)\left(t_{k+1}-t_{j+1}\right)+{\left(t_{k+1}-t_j\right)}^{\alpha +1}\right]\right]=$$

$$=\frac{{\tau }^{\alpha +1}}{\alpha \left(\alpha +1\right)}\left[{\left(k+1-j\right)}^{\alpha +1}-\left(\alpha +1\right){\left(k-j\right)}^{\alpha }-{\left(k-j\right)}^{\alpha +1}\right]=$$

$$=\frac{{\tau }^{\alpha +1}}{\alpha }\left[\frac{{\left(k+1-j\right)}^{\alpha +1}-\left(\alpha +1\right){\left(k-j\right)}^{\alpha }-{\left(k-j\right)}^{\alpha +1}}{\alpha +1}\right].$$

Then the right-hand side of estimate (31) will take the form, taking into account Theorem 97 from [71]:

$${∥z^{\prime }∥}_{\infty }\frac{{\tau }^{\alpha +1}}{\alpha }\sum _{j=0}^k\left[\frac{{\left(k+1-j\right)}^{\alpha +1}-\left(\alpha +1\right){\left(k-j\right)}^{\alpha }-{\left(k-j\right)}^{\alpha +1}}{\alpha +1}\right]=$$

$$={∥z^{\prime }∥}_{\infty }\frac{{\tau }^{\alpha +1}}{\alpha }\left[\frac{{\left(k+1\right)}^{\alpha +1}}{\alpha +1}-\sum _{j=0}^k{j}^{\alpha }\right]=$$

$$={∥z^{\prime }∥}_{\infty }\frac{{\tau }^{\alpha +1}}{\alpha }\left[\underset{0}{\overset{k+1}{\int }}{t}^{\alpha }-\sum _{j=0}^k{j}^{\alpha }\right]={∥z^{\prime }∥}_{\infty }\frac{{\tau }^{\alpha +1}}{\alpha }{\left(k+1\right)}^{\alpha }=\frac{1}{\alpha }{∥z^{\prime }∥}_{\infty }{t}^{\alpha +1}\tau .$$

□

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:

$$\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^{\prime \prime }∥}_{\infty }{t}_{k+1}^{\alpha }{\tau }^2,$$

(32)

Proof. 

The proof of estimate (32) is carried out by analogy with the proof of Lemma 1. □

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)$ and

$$\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 ,$$

(33)

$$\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,$$

(34)

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

$$\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.$$

(35)

Proof of Theorem 3. 

We will carry it out using the method of mathematical induction. Estimates (31) and (32) follow from Lemma 1 and Lemma 2, respectively. For simplicity, we will omit the index i, and also represent the right-hand sides of system (22) as a function $f\left(x,t\right)=\left\{f_1\left(x,t\right),f_2\left(x,t\right)\right\}$, which satisfies the Lipschitz condition:

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

where $L>0$ is the Lipschitz constant.

Let us show that the estimate holds:

$$\left|x\left(t_j\right)-x_j\right|\le C{\tau }^q,\forall j\in \left\{0,1,...,N\right\},$$

(36)

where C is an arbitrary constant. Then in the case when $j=0$ for (36) we get $\left|x\left(t_0\right)-a\right|=0$, because the initial conditions are specified exactly. Let us assume that estimate (36) is true for $j\in 0,1,...,k,k\le N$. Let us prove that estimate (36) is true for $j=k+1$. Let us estimate the error of the predictor $x_{k+1}^p$:

$$\left|x\left(t_{k+1}\right)-x_{k+1}^p\right|=\frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}f\left(x\left(t\right),t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}f\left(x_j,t_j\right)\right|=$$

$$=\frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}{\partial }_{0t}^{\alpha }x\left(t\right)dt+\left(\sum _{j=0}^k{\theta }_{j,k+1}f\left(x_j,t_j\right)-\sum _{j=0}^k{\theta }_{j,k+1}{\partial }_{0t}^{\alpha }x\left(t_j\right)\right)-\sum _{j=0}^k{\theta }_{j,k+1}f\left(x_j,t_j\right)\right|\le$$

$$\le \frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}{\partial }_{0t}^{\alpha }x\left(t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}{\partial }_{0t}^{\alpha }x\left(t_j\right)\right|+\frac{1}{\Gamma \left(\alpha \right)}\left|\sum _{j=0}^k{\theta }_{j,k+1}f\left(x\left(t_j\right),t_j\right)-\sum _{j=0}^k{\theta }_{j,k+1}f\left(x_j,t_j\right)\right|\le$$

$$=\frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}{\partial }_{0t}^{\alpha }x\left(t\right)dt-\sum _{j=0}^k{\theta }_{j,k+1}{\partial }_{0t}^{\alpha }x\left(t_j\right)\right|+\frac{1}{\Gamma \left(\alpha \right)}\sum _{j=0}^k{\theta }_{j,k+1}\left|f\left(x\left(t_j\right),t_j\right)-f\left(x_j,t_j\right)\right|.$$

Taking into account estimate (33) and taking into account that the function $f\left(x,t\right)$ is Lipschian and estimate (36) for $j=k$ is satisfied, then the last expression does not exceed the quantity:

$${\displaystyle \frac{C_1{t}_{k+1}^{\alpha }\tau }{\Gamma \left(\alpha \right)}}+{\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}\sum _{j=0}^k{\theta }_{j,k+1}LC{\tau }^q,$$

where $L>0$ is the Lipschitz constant. Note that the weighting coefficient has the property: $\sum _{j=0}^k{\theta }_{j,k+1}=\frac{1}{\alpha }{t}_{k+1}^{\alpha }\le \frac{T}{\alpha }$, since $t_{k+1}\le T$. Therefore, we finally have an estimate:

$$\left|x\left(t_{k+1}\right)-x_{k+1}^p\right|\le \frac{C_1{T}^{\alpha }\tau }{\Gamma \left(\alpha \right)}+\frac{LC{T}^{\alpha }{\tau }^q}{\Gamma \left(\alpha +1\right)}.$$

(37)

Let’s get an estimate for the corrector.

$$\left|x\left(t_{k+1}\right)-x_{k+1}\right|=\frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}f\left(x\left(t\right),t\right)dt-\sum _{j=0}^k{\rho }_{j,k+1}f\left(x_j,t_j\right)-{\rho }_{k+1,k+1}f\left(t_{k+1},x_{k+1}^p\right)\right|=$$

$$\begin{array}{cc}& =\frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}{\partial }_{0t}^{\alpha }x\left(t\right)dt+\left(\sum _{j=0}^k{\rho }_{j,k+1}f\left(x\left(t_j\right),t_j\right)+{\rho }_{k+1,k+1}f\left(x_{k+1},t_{k+1}\right)-\sum _{j=0}^{k+1}{\rho }_{j,k+1}f\left(x\left(t_j\right),t_j\right)\right)-\right.\hfill \\ & -\left.\sum _{j=0}^{k+1}{\rho }_{j,k+1}f\left(x_j,t_j\right)-{\rho }_{k+1,k+1}f\left(x_{k+1}^p,t_{k+1}\right)\right|\le \frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}{\partial }_{0t}^{\alpha }x\left(t\right)dt-\sum _{j=0}^{k+1}{\rho }_{j,k+1}f\left(x\left(t_j\right),t_j\right)\right|+\hfill \\ & +\frac{1}{\Gamma \left(\alpha \right)}\left|{\rho }_{k+1,k+1}f\left(x_{k+1},t_{k+1}\right)-{\rho }_{k+1,k+1}f\left(x_{k+1}^p,t_{k+1}\right)\right|+\frac{1}{\Gamma \left(\alpha \right)}\left|\sum _{j=0}^k{\rho }_{j,k+1}f\left(x\left(t_j\right),t_j\right)-\sum _{j=0}^k{\rho }_{j,k+1}f\left(x_j,t_j\right)\right|=\hfill \\ & =\frac{1}{\Gamma \left(\alpha \right)}\left|\underset{0}{\overset{t_{k+1}}{\int }}{\left(t_{k+1}-t\right)}^{\alpha -1}{\partial }_{0t}^{\alpha }x\left(t\right)dt-\sum _{j=0}^{k+1}{\rho }_{j,k+1}f\left(x\left(t_j\right),t_j\right)\right|+\frac{1}{\Gamma \left(\alpha \right)}{\rho }_{k+1,k+1}\left|f\left(x_{k+1},t_{k+1}\right)-f\left(x_{k+1}^p,t_{k+1}\right)\right|+\hfill \\ & +\frac{1}{\Gamma \left(\alpha \right)}\sum _{j=0}^k{\rho }_{j,k+1}\left|f\left(x\left(t_j\right),t_j\right)-f\left(x_j,t_j\right)\right|.\hfill \end{array}$$

Due to the fact that the weighting coefficient ${\rho }_{k+1,k+1}=1$ and for $j=k+1$ the function $f\left(x,t\right)$ satisfies Lipschitz condition, we obtain the following estimate:

$$\left|x\left(t_{k+1}\right)-x_{k+1}\right|\le \frac{C_2{t}_{k+1}^{\alpha }{\tau }^2}{\Gamma \left(\alpha \right)}+\frac{1}{\Gamma \left(\alpha \right)}\sum _{j-0}^k{\rho }_{j,k+1}L\left|x\left(t_j\right)-x_j\right|+\frac{{\tau }^{\alpha }}{\Gamma \left(\alpha +2\right)}L\left|x\left(t_{k+1}\right)-x_{k+1}^p\right|.$$

(38)

By estimate (37) and ${\displaystyle \sum _{j=0}^k}{\rho }_{j,k+1}=\frac{1}{\alpha }{t}_{k+1}^{\alpha }\le \frac{T^{\alpha }}{\alpha }$, we get:

$$\begin{array}{cc}& \left|x\left(t_{k+1}\right)-x_{k+1}\right|\le \frac{C_2{T}^{\alpha }{\tau }^2}{\Gamma \left(\alpha \right)}+\frac{LC{T}^{\alpha }{\tau }^q}{\Gamma \left(\alpha +1\right)}+\frac{L{\tau }^{\alpha }}{\Gamma \left(\alpha +2\right)}\left(\frac{C_1{T}^{\alpha }\tau }{\Gamma \left(\alpha \right)}+\frac{LC{T}^{\alpha }{\tau }^q}{\Gamma \left(\alpha +1\right)}\right)\le \hfill \\ & \le \left(\frac{C_2{T}^{\alpha }}{\Gamma \left(\alpha \right)}+\frac{LC{T}^{\alpha }}{\Gamma \left(\alpha +1\right)}+\frac{L{C}_1{T}^{\alpha }}{\Gamma \left(\alpha +2\right)\Gamma \left(\alpha \right)}+\frac{L^{2}C{T}^{\alpha }{\tau }^{\alpha }}{\Gamma \left(\alpha +1\right)\Gamma \left(\alpha +2\right)}\right){\tau }^q.\hfill \end{array}$$

The last inequality follows from the fact that $0<\alpha <1$ and $q>1+min\left({\alpha }_1,{\alpha }_2\right)$.

By choosing T small enough, we can be sure that the second term in parentheses will be limited to $C/2$. Having fixed this value for T, we can then make the sum of the remaining expressions in parentheses less than $C/2$ by choosing C large enough. Then we obtain the necessary estimate (36). □

Remark 19.

Note that in the case when ${\alpha }_1={\alpha }_2=1$ we have the second order of convergence according to estimate (35), 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), (28) using Runge’s rule (double recalculation method) [72]. For these purposes, we introduce the following calculation formulas:

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

(39)

where $q=1+min\left({\alpha }_1,{\alpha }_2\right)$ is the theoretical order of accuracy of the ABM method according to estimate (2.19), ${\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 formulas:

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

(40)

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’s 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^{*}=1.35,y^{*}=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’s consider the case ${\alpha }_1=0.9,{\alpha }_2=0.8$, 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), (28).

3.4. Simulation results

Let’s 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’s consider the case when in the numerical ABM algorithm (27), (28) the parameter values are: ${\alpha }_1={\alpha }_2=1$, ${\delta }_1={\delta }_2=0$ and n is a constant. This leads us to model (10), which was considered by S.V. Dubovsky in his work [14]. Let’s choose the values of the remaining parameters: $n=0.2,x^{*}=1.3,y^{*}=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^{*}=1.3,y^{*}=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 damp (dissipation) and on the phase plane the rest point $\left(x^{*},y^{*}\right)$ is a stable focus. Let’s give another example. Let’s fix the value ${\alpha }_1$ and change the values ${\alpha }_2=\left[0.8,0.9,1\right]$.

Here is the same as in Figure 12, in Figure 13. we see damped oscillations when changing the values of the parameter ${\alpha }_2$, rest point $\left(x^{*},y^{*}\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 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 [35].

Example 5.

Fractional mathematical model S.V. Dubovsky with parameter values [29]: $n=0.2,\lambda =2.25,$$x^{*}=1.3,y^{*}=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 here 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 get 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 system (22), i.e. It takes time to stabilize the oscillatory modes.

Let’s consider the following more general case (Figure 17). Let’s choose the 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 depends 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, 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 fo r 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.

In 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 [73], 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 the article [74], 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$:

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

(41)

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 (41) 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 [75,76].

4.2. Problem Statement and Solution Method

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

$$\left\{\begin{array}{cc}\hfill & {\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^{*}\right)+f_1\left(t\right),\hfill \\ \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^{*}\right)+f_2\left(t\right),\hfill \\ \hfill & x\left(0\right)=a,y\left(0\right)=b.\hfill \end{array}\right.$$

(42)

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 (42) becomes the Cauchy problem (22).

Remark 23.

We will call the Cauchy problem (42) 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 [77,78,79,80,81].

For the Cauchy problem (42), a generalization of the numerical ABM method (27), (28) is proposed [74]. To do 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:

$$\left\{\begin{array}{cc}\hfill & x_{k+1}^p=x_0+\frac{{\tau }^{{\alpha }_{1,k+1}}}{\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^{*}\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right),\hfill \\ \hfill & y_{k+1}^p=y_0+\frac{{\tau }^{{\alpha }_{2,k+1}}}{\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^{*}\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right),\hfill \\ \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.$$

(43)

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

$$\left\{\begin{array}{cc}\hfill & x_{k+1}=x_0+\frac{{\tau }^{{\alpha }_{1,k+1}}}{\Gamma \left({\alpha }_{1,k+1}+2\right)}\left(\begin{array}{cc}\hfill & \left(-\lambda n_{k+1}{x}_{k+1}^p\left(x_{k+1}^p-1\right)\left(y_{k+1}^p-y^{*}\right)+{\delta }_{1}cos\left({\omega }_1\tau \left(k+1\right)\right)\right)+\hfill \\ \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^{*}\right)+{\delta }_{1}cos\left({\omega }_{1}j\tau \right)\right)\hfill \end{array}\right),\hfill \\ \hfill & y_{k+1}=y_0+\frac{{\tau }^{{\alpha }_{2,k+1}}}{\Gamma \left({\alpha }_{2,k+1}+2\right)}\left(\begin{array}{cc}\hfill & \left(n_{k+1}\left(1-n_{k+1}\right){\left(y_{k+1}^p\right)}^2\left(x_{k+1}^p-x^{*}\right)+{\delta }_{2}cos\left({\omega }_2\tau \left(k+1\right)\right)\right)+\hfill \\ \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^{*}\right)+{\delta }_{2}cos\left({\omega }_{2}j\tau \right)\right)\hfill \end{array}\right),\hfill \end{array}\right.$$

(44)

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

$${\rho }_{j,k+1}^i=\left\{\begin{array}{cc}& 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.

Lemma 3.

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

$$\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 ,$$

(45)

where ${∥ ∥}_{\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:

$$\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^{\prime \prime }∥}_{\infty }{t}_{k+1}^{{\alpha }_{k+1}}{\tau }^2,$$

(46)

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

$$\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 ,$$

(47)

$$\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,$$

(48)

then the ABM method (43) and (44) converges with the order

$$\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).$$

(49)

Let us evaluate the computational accuracy of the ABM method (43), (44) using Runge’s rule (39) and (40, which is compared with the theoretical one. Good agreement between them 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:

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

(50)

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 (43), (44).

Example 8.

$l_1=0.6,l_2=1.1,{\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 the theoretical one within a finite number of iterations of the ABM method. The latter indicates the adequacy of estimate (49).

4.4. Simulation Results

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

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 (50). Parameter values: ${\delta }_1={\omega }_1=0$, ${\delta }_2=0.5$, ${\omega }_2=1.5$, $n=0.2$, $\lambda =2.25$, $x^{*}=1.3,y^{*}=0.5$, $a=1.35,b=0.5$. For mathematical modeling we choose: $T=500,\tau =0.4$. To find a numerical solution to system (42), a function was written in MATLAB (Appendix 6.1) [36]:

$$\mathrm{ABMDubovskiyFracLine}(\mathrm{T},\mathrm{tau},\mathrm{alpha},\mathrm{delta},\mathrm{phi},\mathrm{par},\mathrm{Init}),$$

(51)

which outputs oscillograms and phase trajectories based on the numerical solution of system (42) using the ABM method (43), (44). Application of function (51) 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 formulas (50).

In Figure 21 show oscillograms for Example 9. Here we can observe that over time the oscillations reach a constant level, and the rest point $\left(x^{*},y^{*}\right)$ for the corresponding closed phase trajectory is called 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’s consider the case when the parameter values are: ${\alpha }_{10}$, ${\alpha }_{20}$, $M_1$, $M_2$ are fixed in formula (50), 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’s consider the case when the parameters ${\delta }_1$ and ${\omega }_1$ change.

Here we also see that the phase trajectories become closed over time. Let’s consider the case for Example 9, when the initial conditions of system (42) 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’s 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 (50) 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 (42) taking into account (14) and (50) can be represented in MATLAB as the function ABMDubovskiyFracLineny.m (see Appendix 6.2) . Let’s consider the case when 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 formulas (50) 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.

In 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 when the initial conditions a and b change (Figure 28).

5. Conclusions

In conclusion, the following conclusions can be drawn from the research results. 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 long waves N.D. Kondratiev 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 find solutions to proposed mathematical models based on the Adams-Bashforth-Moulton method from the predictor-corrector family. For such methods, an error analysis was carried out, the results were formulated in the form of lemmas and theorems. Next, using test examples and Runge’s rule, estimates of the computational accuracy 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 has only damped oscillations, i.e. there are no self-oscillations. The presence of external influence makes it possible for the existence of a limit cycle, 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 function of the rate of accumulation on capital productivity can also influence the shape of the limit cycle. Complex oscillatory regimes deserve separate study and, accordingly, economic interpretation.

It is necessary to note that further development of the mathematical models proposed in this article by S.V. Dubovsky is to take into account the dependence of the orders of fractional derivatives not only on time, but also on the solution function by analogy with the article [82].