Modeling of fracture concentration by Sel’kov fractional dynamic system

Microseismic phenomena are studied by a Sel’kov generalized nonlinear dynamic system. This system is mainly applied in biology to describe substrate and product glycolytic oscillations. Thus, Sel’kov dynamic system can also describe interaction of two types of fractures in an elastic-friable medium. The first type includes seed fractures with lower energy and the second type are large fractures which generate microseisms. The first type of fractures are triggers for the second type of fractures. However opposite transition is possible. For example, when large fractures lose their energy and partially become seed ones. After their concentration increase, the process repeats providing auto oscillation character of microseism sources. Generalization of Sel’kov dynamic system is its analogue which is based on hereditarity. Hereditarity is studied within hereditary mechanics and it shows that a dynamic system can “remember” for some time the impact which was made upon it. It is typical for viscoelastic and yielding mediums. The Sel’kov generalized dynamic system will be called Sel’kov fractional dynamic system as long as from the point of view of mathematical description, it can be represented in the form of a system of differential equations with fractional derivatives. Fractional derivative orders are associated with system hereditarity and are responsible for energy dissipation intensity emitted by firstand second-type fractures. In the paper, the Sel’kov fractional dynamic model was numerically solved by Adams-Bashforth-Moulton method. Oscillograms and phase trajectories were plotted. It was shown that fractional dynamic model may have relaxation and damped oscillations.


Introduction
In the paper [1] the authors suggested an interesting approach to describe fraction interactions in an elastic-friable medium. The approach is based on application of Sel'kov nonlinear dynamic system which is investigated in Biology [2].
To the author's opinion, the Sel'kov dynamic system describes well the interaction of two types of fractures. The first type is seed fractures with lower energy which develop into the second type of fractures of greater energy when they reach a critical level of their concertation. Fractures of the second type are the sources of microseismic phenomena (oscillations) and they partially change to seed fracture after energy output. Then this auto oscillation process repeats.
In the present paper we consider the Sel'kov generalized dynamic system which takes into account the hereditary effect. Hereditarity means that a system can "remember" for some time the impact which was made upon it. Thus, the process of interaction between two types of fractures may be quite slow. Hereditary effects were firstly considered in Hereditary Mechanics to describe viscoelastic and yielding mediums [3].
Mathematical description of hereditarity is based on Volter-type integro-differential equations with difference kernels in the expression under integral sign which are sometimes called hereditary functions [4]. Hereditary functions can be selected on the basis of experimental data, medium properties or a process under investigation. However, in a first approximation, power hereditary functions are usually selected to describe the hereditary effect the influence from which decreases with time. Power hereditary functions give the possibility to apply mathematical apparatus of fractional calculus and to turn from integro-differential equations to the equations with fractional derivatives [5][6][7]. Dynamic systems which are described by fractional derivative are called fractional dynamic systems or dynamic systems of fractional order [8].
Owing to the necessity to study the auto oscillation character of interaction of two types of fractures, the main aim of the paper is to state the possibility of existence of auto oscillation modes within the framework of Sel'kov fractional dynamic system (SFDS).

Some information from the theory of fractional calculus
In this section we consider the main definitions from the theory of fractional calculus. Its aspects can be studied in detail in the books [5][6][7]. Definition 1. Riemann-Liouville fractional integral of order α: where Γ (·) is the gamma function. Riemann-Liouville fractional integral (1) has the following properties: . Definitiion 2. Gerasimov-Caputo fractional derivative of order α has the form: Operator (2) has the following properties:

Problem statement
We consider the following nonlinear dynamic system: where x (t) is the function determining the concentration of seed first-type fractures; y (t) is the function determining the concentration of second-type fractures which generate microseisms, t ∈ [0, T ] is the coordinate responsible for current time of the process, T > 0 is the constant, modeling time;x 0 , y 0 , v, a, b are given positive constants; fractional differentiation operators are understood in the sense of Gerasimov-Caputo of orders 0 < α 1 , α 2 < 1 and are determined according to (2). Remark 1. Dynamic system (3) is a generalization of the known Sel'kov dynamic system which is used in Biology to describe substrate and product glycolytic oscillations. When the parameters α 1 = α 2 = 1, the systems coincide.
The aim of the investigation of the paper is to find the SFDS (3) solution for given values of the parameters x 0 , y 0 , v, a, b and to determine the conditions for existence of auto oscillations.
Before we turn to the method of SFDS (3) solution, we introduce some definitions form the theory of qualitative analysis of dynamic systems to study asymptotic stability of its equilibrium points.

Asymptotic stability of equilibrium points
Definition 4. Sel'kov dynamic system is commensurable if fractional derivative orders are equal between each other α 1 = α 2 = α, otherwise it is incommensurable.
Definition 5. SFDS (3) has one equilibrium point E = (x * , y * ) which is the solution of the following system: and has the form: Definition 6. Jacobi matrix of SFDS (3) is the matrix of the form: We should note that taking into account equilibrium point (4), Jacobian (5) has the form: Examine equilibrium point (4) for asymptotic stability. In order to do that, we apply the following theorems suggested in the paper [9]. Assume that m is the least common multiple of fractional order denominators and are calculated by the standard equation For the incommensurable SFDS the following theorem is true. Theorem 2. Equilibrium point (4) is asymptotically stable for the incommensurable system (3), where λ i are calculated from the standard equation (8).
Assume that f 1 (x, y) , f 2 (x, y) are right parts of system (3) under consideration. Proposition 1. (Criterion for the existence of a closed trajectory (Bendikson criterion)). If in some simply connected region the function f (x, y) = ∂ f 1 ∂x + ∂ f 2 ∂y does not change its sign and is not identically zero, than there are no closed contours composed of trajectories in this region.
Proposition 2.(Criterion for existence of chaotic modes [9]). If the condition M = π 2m − min i arg |λ i | < 0 is fulfilled, than there are no chaotic modes in the dynamic system. Proposition 2 is a sequence of Theorem 2.

Solution method
As a method for solution of SFDS (3) we choose Adams-Bashforth-Moulton (ABM) which refers to numerical predictor-corrector methods. ABM method was studied in detail and discussed in the papers [10][11][12]. We adapt this method to solve SFDS (3). To do that we use definitions (1), (2) and their properties. On a uniform mesh with the step τ = T /N we introduce the functions x p n+1 , y p n+1 , n = 0, . . . , N − 1 which are determined by Adams-Bashforth formula (predictor) as and the functions x n+1 , y n+1 which are determined by Adams-Moulton formula for the corrector where the weighting coefficients in (10) are determined by the formula It is known that for ABM method (x 1 = x (t) , x 2 = y (t) , i = 1, 2) the error estimation max 1≤ j≤k Remark 2. We should note that in a classical case α i = 1 we obtain a classical ABM method of the second-order accuracy.
In this section, applying qualitative analysis and numerical modeling, we show that auto oscillation modes are possible for SFDS (3). Example 1. (Sel'kov classical dynamic system) When α 1 = α 2 = 1 in SFDS (3), than we come to the classical model of Sel'kov dynamic system [1]. We choose the rest parameters of system (3) from the paper [1]: a = 0.116, b = 0.602, v = 0.64421. Parameter v was estimated from the trace of the Jacobian (6) equal to zero to obtain sustained oscillations. In this case, standard equation (6) for m = 1 has pure imaginary roots λ 1,2 = ±0.6048437164I, thus, the stationary point E = (0.6442125706, 1.760933066) may be a center of a focus. We can show that according to Bendikson criterion, the function f (x, y) = −0.602x 2 + 1.204xy − 1.116 changes its sign in a simply connected region. Thus, the stationary point E = (0.6442125706, 1.760933066) is the center.
It is known that a stationary point of center type is stable according to Lyapunov but it is not asymptotically stable. Indeed, according to Theorem 1, condition (7) is not fulfilled for the commensurable system (3) when α 1 = α 2 = 1. However, phase trajectories reach stable boundary cycle just like in the paper [1].   (7). Thus, the equilibrium point is asymptotically stable. We can show that according to Bendikson criterion, the function f (x, y) = −1.3x 2 + 2.6xy − 1.03 gives sign change in a simply connected domain, and the measure of existence of chaotic modes is negative, M = −0.0050548591 < 0. Thus, we can make a conclusion that the system has a closed phase trajectory and does not have chaotic modes. Owing to the fact that the point is asymptotically stable, closed trajectory is a stable boundary cycle (Figure 2). Example 3.(SFDS incommensurable case). We consider an incommensurable SFDS (3) α 1 = 1, α 2 = 0.9. We choose other parameter values to be a = 0.12, b = 1.1, v = 0.6. Then the equilibrium point (4) has the coordinates E = (0.6, 1.162790698). Investigate this point for asymptotic stability. According to Theorem 1, standard equation (7), taking into account that α 1 = 10 10 , α 2 = 9 10 i.e. m = 10, takes the form λ 19 − 0.0188837209λ 9 + 0.516 = 0. All the roots of this equation satisfy condition (9) of Theorem 2. Thus, the equilibrium point is asymptotically stable.
Owing to Bendikson criterion, the function f (x, y) = −1.1x 2 + 2.2xy − 1.12 changes its sign in a simply connected domain and, consequently, the system may have a closed trajectory. According to the criterion for existence of chaotic modes for Example 3, M = −0.0068640576 < 0, that is why they are absent. Thus, we can make a conclusion that there is a stable closed trajectory, stable boundary cycle.   (10), (11) for N = 2000. We see the boundary cycle. We can show that this cycle is stable. A similar boundary cycle was obtained by the authors of the paper [1] but within the framework of Sel'kov classical dynamic system.

Acknowledgment
The work was carried out according to the Subject AAAA-A17-117080110043-4" «Dynamics of physical processes in the active zones of near space and geospheres» IKIR FEB RAS.

Conclusion
In the present paper, we have showed that with the help of elements of qualitative analysis of fractional dynamic systems and numerical modeling, Sel'kov fractional dynamic system may have auto oscillation modes. It is also interesting to investigate chaotic modes, for example, by constructing Lyapunov maximum exponent spectra similar to the papers [13][14][15].