CC BY
53
UDC 004.94, 519.21
One-Step Stochastic Processes Simulation Software Package
E. G. Eferina, A. V. Korolkova, M. N. Gevorkyan, D. S. Kulyabov, L. A. Sevastyanov
Department of Applied Informatics and Probability Theory Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, Russian Federation, 117198
It is assumed that the introduction of stochastic in mathematical model makes it more adequate. But there is virtually no methods of coordinated (depended on structure of the system) stochastic introduction into deterministic models. Authors have improved the method of stochastic models construction for the class of one-step processes and illustrated by models of population dynamics. Population dynamics was chosen for study because its deterministic models were sufficiently well explored that allows to compare the results with already known ones.
To optimize the models creation as much as possible some routine operations should be automated. In this case, the process of drawing up the model equations can be algorithmized and implemented in the computer algebra system. Furthermore, on the basis of these results a set of programs for numerical experiment can be obtained.
The computer algebra system Axiom is used for analytical calculations implementation. To perform the numerical experiment FORTRAN and Julia languages are used. The Runge-Kutta method for stochastic differential equations is used as numerical method.
The program complex for creating stochastic one-step processes models is constructed. Its application is illustrated by the predator-prey population dynamic system.
Computer algebra systems are very convenient for the purposes of rapid prototyping in mathematical models design and analysis.
Key words and phrases: stochastic differential equations; "predator-prey" model; master equation; Fokker-Planck equation; computer algebra software; Axiom system.
1. Introduction
This work corresponds our research on mathematical models stochastization. This item is interesting due to the following problems: the construction of population models from first principles and the introduction of the stochastic into such models (the population dynamics is studied because of similar models introduction in other areas).
The problem of stochastic term introduction arises during mathematical models stochastization. There are several ways to solve this problem. The easiest option is an in the deterministic equation. But when additive stochastic term is introduced some free parameters that require further definition appears. Furthermore, these stochastic terms usually interpreted as an external (rather than structural) random impact. In this regard, we used and improved the stochastic one-step processes models construction method, based on master equation [1,2]. Stochastic differential equation is considered as its approximate form. It allows to get the model equations from general principles. Furthermore, deterministic and stochastic parts are derived from the one equation so we can regard it as stochastic and deterministic parts consistency.
The aim of this work is the software complex development for rapid prototyping construction of stochastic one-step processes models. This complex consists of two blocks. The first block generates the equations of dynamic stochastic process model on the principles similar to chemical kinetic relations describing the investigated process. This block is implemented by means of the computer algebra system — system FriCAS, which is offshoot of Axiom.
The second block is used for the numerical analysis of the resulting model. For numerical solution of deterministic and stochastic models equations some Runge-Kutta different orders methods [3,4] are used.
Received 25th June, 2014.
To illustrate the developed system the well-known population model predator-prey is used [5-7].
The structure of the paper is as follows. The basic notation and conventions are introduced in Section 2. Section 3 is devoted to brief introduction to the one-step processes stochastization method. Further, in the Section 4 the model under investigation is described. In the subsection 4.1 there is a brief reference to standard (deterministic) approach, and in the subsection 4.2 the stochastic extension of our model with the help of the one-step processes stochastization method is obtained.
In the Section 5.1 we justify selection of the system, which implements the model equations generating unit. The actual interface of this part of the program complex is described in the Section 5.2.
The possibility of applying Runge-Kutta methods for the analysis of stochastic differential equations is considered in the Section 6. The software interface of the model equations numerical analysis unit is also described in this section. Calculations example is based on the predator-prey model.
2. Notations and Conventions
1. We use abstract indices notation [8]. In this notation tensor as a whole object is denoted just as an index (e.g., xz), components are denoted by underlined index (e.g^ x-).
2. We will adhere to the following agreements. Latin indices of the middle of the alphabet (i, j, k) will apply to the space of the system state vectors. Latin indices from the beginning of the alphabet (a) will relate to the Wiener process space. Latin indices from the end of the alphabet (p, q) will refer to the indices of the Runge-Kutta method. Greek indices (a) will set a number of different interactions in kinetic equations.
3. A Dot over a symbol denotes differentiation with respect to time.
4. The comma in the index denotes partial derivative with respect to corresponding coordinate.
3. One-Step Processes Modeling
Let's briefly review the method of one-step processes stochastization on the basis of [9].
We understand one-step processes as Markov processes with continuous time with values in the domain of integers, which transition matrix allows only transitions between neighbouring portions. Also, these processes are known as birth-and-death processes.
One-step processes are subject to the following conditions:
1. If at the moment t the system is in state i G Z^0,then the probability of transition to state i + 1 in time interval [t, t + At] is equal to k+At + o(At).
2. If at time moment t the system is in state i G Z+,then the probability of transition to state i — 1 in the time interval [t, t + At] is equal to k-At + o(At).
3. The probability of transition to a state other than the neighbouring is equal to o(At).
4. The probability to remain in the same state is equal to 1 — (k+ + k-)At + o(At).
5. State i = 0 is an absorbing boundary.
The idea of the one-step processes stochastization method is as follows. Based on the patterns of interaction we construct a master kinetic equation, expand it into a series, leaving only the terms up to and including the second derivative. The resulting equation is the Fokker-Planck equation. In order to get more convenient model we record corresponding Langevin equation. In fact, as we shall see , from the patterns of interaction we will immediately obtain the coefficients of the Fokker-Planck equation
(and accordingly, the Langevin equation), so for practical use of the method there is no need to construct the master kinetic equation.
3.1. Interaction Schemes
We will describe the state of the system by a state vector x1 £ Rn, where n is the system dimension (the state vector is considered as the set of mathematical values, fully describing system). The operator nj £ Z>0 x Z>0 defines the state of the system before the interaction and the operator mj £ Z>0 x Z>0 — after the interaction. The result of interaction is a system transition to another state [1,10].
There are s kinds of different interactions that may happen in the system, s £ Z+. So, instead of nj and mj let's consider the operators nja £ Z>0 x Z>0 x Z>0 and mf £ Z>0 x Z>0 x Z>0. "
System elements interaction will be described with the interaction schemes similar to chemical kinetic schemes [11]:
k+
nfxj ^ mfxj. (1)
k-
ka
Here Greek indices specify the number of interactions and Latin ones specify dimensionality of the system. The state change is given by the operator
rja = mja — nf. (2)
Thus, one-step interaction a in forward and opposite directions can be written as
xl —> x1 + rl~Xj,
(3)
i i i ^ j X ^ x — r~ XJ .
We can write (1) not in the form of vector equations, but in the more traditional form sums:
.. fc+ . .
nfx3Si ^ mjaxjSi, (4)
k-
where Si = (1,..., 1).
Also, we will use the following notation:
nia := nf5j, mia := mf5j, ria := rf5j. (5)
3.2. Master Equation
Transition probabilities per unit of time from the state x' to the state xi + r^xj
(to the state xi — r~xj) are proportional to the number of xi combination from a set of nia elements (of x1 — combinations from a set of mia) and are given by:
n
Sa = ^a H , - iaM , X , (x- - n- )■
i=1 v '
- (6)
n i, '
- = U- IT X ■
sa = / - iax . .
, (x- - m-)■
Thus, the general form of the master kinetic equation for the states vector x1 (it changes by r%—xJ per step), takes the form:
dp(xi ,t) dt
m
L
a=1
s-(xi + r% t)p(xi + r% t) - s^(xi)p(xi,t)
+
+
s+(xi - ri»,t)p(xi - ri«,t) - s--(xi)p(xi,t) \. (7)
3.3. Fokker-Planck Equation
With the help of the Kramers-Moyal expansion, the Fokker-Planck equation [11] is obtained. For this purpose we will make several assumptions:
1) there are only small jumps, ie sa(xl) is a slowly varying function with the change
of xl]
2) p(xl,t) also slowly changes with the change of xz.
Then in Fokker-Planck equation (7) one can shift from the point (x1 ± r%~xJ) to the
point x1, and by expanding the right-hand side in a Taylor series and dropping terms of order higher than the second, we obtain Fokker-Planck equation:
^ = —dz [A*p] + 1 dzd, [B*p] , (8)
where _ _
(^k
Ai := Ai(xk,t) = r^
Sa Sa
Bi j := Bij(xk,t)= r^
Sa Sa
(9)
a = 1, m.
As seen from (9), the coefficients of the Fokker-Planck equation can be obtained directly from (2) and (6), i.e. in practical calculations, there is no need to write the master equation.
3.4. Langevin Equation
The Langevin equation corresponds to the Fokker-Planck equation:
d^ = aidt + bldWa, (10)
where a1 := al(xk,t), bla := bla(xk,t), x1 G 1" — is the system state vector, Wa G rto — m-dimensional Wiener process. Wiener process is implemented as d^ = sVdt, where s ~ N(0,1) is normal distribution with average 0 and variance 1. Latin indices from the middle of the alphabet denote the values related to the state vectors (dimension of the space is n), and Latin indices from the beginning of the alphabet denote the values related to the Wiener process vector (dimension of the space is m < n).
The connection between the equation (8) and (10) expressed by the following relationships:
Ai = a\ Bij = VaVa. (11)
We will use Ito interpretation. Under the Ito interpretation, differential of complex functions does not obey the standard formulas of analysis. To calculate it rule or Ito lemma are used.
Let f := f (xk, t) is a function of a random process xk(t), f G C2. Then the formula of the differential is [12]:
df
dtf + a*/,. + ttb^aftij
dt + bïf,i d^a, (12)
where f := f(xk,t), a* := ai(xk,t), Va := bàa(xk,t), and d^a := d^a(i).
4. Predator—Prey Model 4.1. Deterministic Predator-Prey Model
Systems with the interaction of two predator-prey populations types are extensively studied and there are a lot of various models for these systems. The very first predator-prey model is considered to be a model which was obtained independently by A. Lotka and V. Volterra. Lotka in [13] described some hypothetical chemical reaction:
A X k Y k B, (13)
where X,Y are intermediates substances, coefficients k\,k2,k3 are rates of chemical reactions, A is a initial reagent, and B is a resultant. As a result was a system of differential equations:
{% = hx - k2xy, (14)
\y = k2xy - k3y.
This system is identical to the system of differential equations, obtained by Volterra, who considered the growth mechanism of two populations with predator-prey interaction type. In order to get equations [5] Volterra made a series of idealized assumptions about nature of intraspecific and interspecific relationships in the predator-prey system.
4.2. Stochastic Predator-Prey Model
Consider a model of predator-prey system, consisting of two individuals species, one of which hunts, second is provided with inexhaustible food resources. Let's introduce the notation, where X is a prey and Y is a predator, then we can write the possible processes (4) for the state vector x- = (X, Y)t [14-17]:
X -b 2X, r-1 = (1, 0)T,
X + Y % 2Y, r-2 = (-1,1)T, (15)
Y k 0, r-3 = (0,-1)T,
which have the following interpretation. The first relation means that the prey which eats the food unit immediately reproduced. The second relation describes the case when predator absorbs the prey and then it is instantaneously reproduced. Only such possibility of prey death is considered. Last ratio is a natural predator death.
All processes are irreversible, so sa =0, and
, , x\ y\ Sl(x,V) = k17-7T7 = k1^
(x - 1)\y\
x\ y\
4(x, y) = k2 x - -1)\ =k2xy, (16)
l x\ \
s l(x, y) = k3 x.J^Ty. = k3 y-
With the help of the formula (8) we have the Fokker-Planck equation:
= -di (Ai(x, y)p(x, y)) + \d,d0 (B"(x, y)p(x, y)) , (17)
where
Ai(x, y) = s+(x, y) ria, Bij(x, y) = s+(x, y)riarja.
As a result:
(18)
(19)
A-(x, y) = Q k\x + ^ k2xy +
/ 0 \ (kxx - k2xy\
+ {-l) k3V =Wy-k,y),
B^(x, y) = Q (1, 0)kix + ((-1,1)k2xy + . f0\in fk 1x + k2xy ~k2xy \
+ -1)k3y = ^ _ k2xy k2xy + k3y).
In order to write a stochastic differential equation in Langevin form (10) for predator-prey model, it is enough to take the square root of the resulting matrix Bi:> in Fokker-Planck equation
ix\ = fkix - k2xy\ a idW^ d\y) = \k2xy - k3V) d + H dW 2) , (20)
b-bia = Ba = fkix +k2xy k2xy x
a ^ -k2xy k2xy + k3y^
It should be noted that the specific form of the matrix bla is not written out because of the extreme awkwardness of the expression. However, with further studies we will need not actually matrix bla, but its square, i.e. the matrix BiJ.
5. Implementation of the One-Step Stochastic Processes Model in the Computer Algebra System
5.1. Justification of the Computer Algebra System Choice
Let's consider systems of analytical calculations, Maxima and Axiom. Maxima is the first system of analytical calculations and it is written in Lisp. Maxima successfully runs on all modern operating systems: Windows, Linux and UNIX, Mac OS and even on PDA running Windows CE/Mobile. Documentation is integrated into the program as a handbook with search. There is no distinction between objects and data
in Maxima, and there is no clear distinction between the operator and function. There is no integrated graphics rendering in the system.
Unlike Maxima Axiom language is strongly typified for better mathematical objects and relationships display. The mathematical basis is written in Spad language. Axiom portability is slightly worse: the system runs under Linux, UNIX, and graphs does not work under Windows. Axiom has its own graphics subsystem.
In 2007 two Axiom open source forks appeared: OpenAxiom and FriCAS. Open Axiom is developed by adhering to the ideology of Axiom, problems that occurred in the Axiom are eliminated. FriCAS developers reorganized the assembly process, expanded functionality. Furthermore, FriCAS supports not only GCL, which operates on limited number of platforms, but ECL, Clisp, sbcl or openmcl, that allows to run FriCAS under wider range of platforms.
5.2. Implementation Description in the Axiom Computer Algebra
System
Method of one-step processes randomization is organized as a module for the FriCAS computer algebra system. To display all the calculations on the screen the variable SHOWCALC:=true is used. To call the method you need to use the main function, which has the following view:
osp(Matrix(Integer), Matrix(Integer), Vector, Vector, Vector)
where the first argument is before interaction states matrix nj, the second argument is after interaction states matrix mj, the third argument is the vector k+, the fourth argument is the vector k-, the fifth argument is the state vector x1. Let's consider the features of the language FriCAS on auxiliary functions. For example, the function calcProd is used to simplify the calculations s+ and s-. In the implementation of the function operator of the condition and built-in function reduce are used:
calcProd : (Matrix(Integer), Vector, Integer, Integer) -> Void calcProd (n, x, a, i) ==
nai:Integer := n(a,i)
if nai = 0 then 1 else reduce(*,[x(i) - j for j in 0..(nai-1)])
In the function Bi intermediate calculations for elements of the matrix B-j are made:
Bi (rv, sp, sm, i) == rv(i) * (transpose rv(i)) * (sp(i) + sm(i))
In order to use the module for predator-prey system model, we call the function with the following arguments:
osp ([[1,0],[1,1],[0,1]],[[2,0],[0,2],[0,0]],
vector([k1,k2,k3]), vector([0,0,0]),vector([x,y]))
Fig. 1 represents the result obtained in TgXmacs shell. In fact, we repeated the results obtained in (19).
Buffer File Edit Insert Text Paragraph Document Project Options Help
j*
T y/* y» *. ** » E (| ) a ® -! p B C 5 B R!
osp ([[1,0] , [1,1]. [0,1]], [[2,0] , [0,2], [0,0]] .vector([kl,k2,k3]) .vector( [0,0, 0]).vector( [x,y]))
FriCAS will attempt to step through and interpret the code.
Cannot compile map: calcProd
We will attempt to interpret the code.
Compiling function Ai with type (ListCHatrix(Integer)),List(
Polynomial(Integer)).List(Polynomial(Integer)).NonNegativelnteger ) -> Matrix(Polynomial(Integer)) Compiling function A with type (List(Matrix(Integer)), List(
Polynomial(Integer)).List(Polynomial(Integer)) , NonNegativelnteger ) -> Matrix(Polynomial(Float))
Compiling function Bi with type (List(Hatrix(Integer)),List(
Polynomial(Integer)),List(Polynomial(Integer)).NonNegativelnteger ) -> Matrix(Polynomial(Integer)) Compiling function B with type (List(Matrix (Integer)), List (
Polynomial (Integer)).List(Polynomial(Integer)) , NonNegativelnteger ) -> Matrix(Polynomial(Float))
( —k2xy + klx \ ( k2xy + klx -A:2xy \ V -Jt3 + it2xy /'V -k2xy k3y + k2xy J
Figure 1. The output of the module for predator—prey model in graphic
TXmacs shell
6. Numerical Experiment for the Program Complex 6.1. Stochastic Runge-Kutta Methods
Euler-Maruyama method is one of well-known numerical methods for solving SDE, it is a special case of a more general Stochastic Runge-Kutta method. Classical Runge-Kutta method can be generalized to the case of the SDE system (10) in the following manner [3,4]:
X- =4 + hB--ai(X[,... + R-Jabi(X[,.. .,X?), x\ =x0 + hrl-ai (X-,..., XI1) + rl-Jabi (X-,..., X.
Indexes k = 1,..., s and I = 1,... ,n refer to stochastic Runge-Kutta method. J ~ N(0, h) or J ~ \fhe, £ ~ N(0,1) are normal distributed random variables. Such a choice of these numerical values for approximation is made because the Wiener process is implemented as d^ = e\/dt. You should also pay attention to double summation in the third term of both numerical scheme formulas as well as the fact that each number J1,..., Jn should generated separately.
The method coefficients, as well as for the classical analogue, can be grouped into a table called the Butcher table:
Rij Rij
ro
For calculations we used a method with the table
0 0 0 0 0 0
2/3 0 0 2/3 0 0
— 1 1 0 — 1 1 0
0 3/4 1/4 0 3/4 1/4
6.2. Software Implementation Description
The purposes of the programs complex were to automate the SDE coefficients A1 and B-j computation with the help of general principles described above, and to find a numerical solution of the equation obtained by means of stochastic Runge-Kutta methods. From a programming standpoint we can derive three subtasks:
1. coefficients A1 and B-j generation using the computer algebra system;
2. generation of source code in languages Fortran and Julia, implementing the SDE on the basis of the coefficients, saved as a text file;
3. writing subroutines/functions implementing stochastic Runge-Kutta methods in Fortran and Julia, and their subsequent compilation together with automatically generated source codes.
As a result of its work Axiom module creates a text file which contains the coefficients A- and B-j in the following form:
# A
A[1]
A[N]
# B
B[1,1] B[1,2] .. B[1,N]
B[N,1] B[N,2] .. B[N,N]
Matrix b\ is calculated numerically with the help of the singular
value matrix decomposition (a subroutine DGESVD from library LAPACK is used).
For the second subtask scripting language Python was chosen (version 3). This language has a wide set of tools to work with strings and text files. Except matrices A- and B-j additional information about the mathematical model was specified as dictionary (standard data type in Python), with model name, list of variables, list of parameters, initial values of variables, parameters values and parameters of the numerical method (integration section and step size).
On the basis of these data, the script automatically generates two files functions.f90 and main.f90, where the first is a module with functions defining the SDE, and the second one is a main program file. While compiling these files the third additional module with auxiliary procedures with Stochastic Runge-Kutta method is added.
6.3. The Numerical Experiment Description
For the programs complex work verification a well-known predator-prey model was chosen with vector a- components
and matrix BiJ:
a1 = ax — ßxy, a2 = —7y + öxy
ax + ßxy —ßxy —ßxy ßxy + 7x
(22) (23)
x is the number of preys, y is the number of predators. Coefficients also have the following physical (biological) meaning: a is the growth rate of the prey population, fi
is a frequency of predators and prey meetings, 7 is an intensity of predators death or migration in a lack of preys, 5 is a predator population growth rate on the assumption of the excess of the prey.
During numerical simulations it was taken into account that the value of variables x,y could not be less than zero (program stop working when one of the variables becomes equal to zero).
Numerical simulation shows that the addition of stochastic to the classical predator-prey model leads to the fact that after a certain time death of one of the competing species comes. So, for the following parameters: a = 10, ft = 1.5, 7 = 8.5, 5 = 1.8 and the initial values: x = 9.7, x = 6.77, victims are first to die, and after that predators die due to lack. This case is illustrated in Fig. 2. For comparison in Fig. 3 is a graph for the deterministic case.
Figure 2. Stochastic predator—prey model, prey die
Figure 3. Deterministic predator—prey model
Under other conditions (a = 10, ft = 1.5, 7 = 8.5, 5 = 0.5, x = 22, y = 6.76) predators die, and the number of victims is increasing rapidly, as for their model assumes an infinite source of food. Graphics for this case are shown in Fig. 4, and Fig. 5 shows for comparison with deterministic case.
Figure 4. Stochastic predator—prey model, predators die
Figure 5. Deterministic predator—prey model
7. Conclusions
This work demonstrates the application of the developed initial physical system formalization method. The system is presented in the form of one or more one-step processes. Formalization of the system is done by introducing the evolution operator. Wherein the analytical description of the model requires a lot of routine operations. To simplify the work we propose to use the computer algebra system (Axiom fork FriCAS).
We have developed an analytical software package block that receives inlet evolution operator and produces the SDE, which describes the original model. For numerical studies of obtained SDE system a second software unit that converts the resulting system of equations into the program code in Fortran and gives its numerical solution was developed. Thus, the software system is applicable for both analytical and numerical study of the original model.
Currently the software package does not cover all possibilities, incorporated in the proposed method of formalizing the original physical system. Since the original system description uses ODE, we should introduce the boundary conditions by ties
or indicator functions. Partial differential equations can help to solve this problem. Further objective is the development of a complete software complex for a method of one-step original physical system model construction.
References
1. A. V. Demidova, D. S. Kulyabov, The Introduction of an Agreed Term in the Equation of Stochastic Population Model, Bulletin of Peoples Friendship University of Russia. Series "Mathematics. Information Sciences. Physics" (3) (2012) 69-78, in Russian.
2. A. V. Demidova, L. A. Sevastyanov, D. S. Kulyabov, Application of Stochastic Differencial Equations to Model Population Systems, in: Third International Conference on Mathematical Modelling of Social and Economical Dynamics MMSED-2010, Russian State Social University, 2010, pp. 92-94, in Russian.
3. K. Debrabant, A. Robler, Classification of Stochastic Runge-Kutta Methods for the Weak Approximation of Stochastic Differential Equations, Mathematics and Computers in Simulation 77 (4) (2008) 408-420. doi:http://dx.doi.org/10.1016/j.matcom.2007.04.016.
4. A. Tocino, R. Ardanuy, Runge-Kutta Methods for Numerical Solution of Stochastic Differential Equations, Journal of Computational and Applied Mathematics 138 (2) (2002) 219-241. doi:http://dx.doi.org/10.1016/S0377-0427(01)00380-6. URL http://www.sciencedirect.com/science/article/pii/ S0377042701003806
5. V. Volterra, Mathematical Theory of the Struggle for Existence, Nauka, 1976.
6. A. D. Bazykin, Nonlinear Dynamics of Interacting Populations, World Scientific, 1998.
7. A. S. Bratus, A. S. Novozhilov, A. P. Platonov, Dynamical Systems and Biology Models, Fizmatlit, M., 2010, in Russian.
8. Roger Penrose and Wolfgang Rindler, Spinors and Space-Time: Two-Spinor Calculus and Relativistic Fields, Vol. 1, Cambridge University Press, 1984.
9. A. V. Demidova, M. N. Gevorkyan, A. D. Egorov, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastyanov, Influence of Stochastization to One-Step Model, Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics" (1) (2014) 71-85, in Russian.
10. A. V. Demidova, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastyanov, The Method of Stochastization of One-Step Processes, in: Mathematical Modeling and Computational Physics, JINR, 2013, p. 67.
11. C. W. Gardiner, Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences, Springer Series in Synergetics, 1985.
12. B. K. 0ksendal, Stochastic Differential Equations: An Introduction with Applications, Springer, Berlin, 2003.
13. A. J. Lotka, Elements of Physical Biology, BiblioBazaar, 2011. URL http://books.google.ru/books?id=tFN9pwAACAAJ
14. A. V. Demidova, The Equations of Population Dynamics in the Form of Stochastic Differential Equations, Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics" (1) (2013) 67-76, in Russian.
15. A. V. Demidova, The Method of Stochastization of Mathematical Models for the Example of the "Predator-Prey", in: A Scientific Session NRNU MEPHI-2013, 2013, p. 127, in Russian.
16. A. V. Korolkova, D. S. Kulyabov, Methods of Stochastization of Mathematical Models on the Example of Peer to Peer Networks, in: A Scientific Session NRNU MEPHI-2013, MEPHI, Moscow, 2013, p. 131, in Russian.
17. A. V. Demidova, D. S. Kulyabov, L. A. Sevastyanov, The Agreed Stochastic Term in Population Models, in: XI Belarusian Mathematical Conference, Institute of Mathematics of the National Academy of Sciences of Belarus, Minsk, 2012, p. 39, in Russian.
УДК 004.94, 519.21
Программный комплекс стохастического моделирования одношаговых процессов
Е. Г. Еферина, А. В. Королькова, М. Н. Геворкян, Д. С. Кулябов, Л. А. Севастьянов
Кафедра прикладной информатики и теории вероятностей Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198
Нашим коллективом разработана методика согласованного (зависящего от структуры системы) введения стохастики в детерминистические модели. На данном этапе методика ограничена классом одношаговых процессов.
Для оптимизации работы по созданию моделей следует автоматизировать как можно больше рутинных операций. В данном случае процесс составления уравнений модели можно алгоритмизировать и реализовать в системе компьютерной алгебры. Кроме того, на базе этих результатов можно получить и набор программ для проведения численного эксперимента.
Для реализации аналитических расчётов используется система компьютерной алгебры Axiom. Для проведения численного эксперимента используются языка FORTRAN и Julia. В качестве численного метода используется метод Рунге—Кутты для стохастических дифференциальных уравнений.
Разработан программный комплекс для создания стохастических моделей одношаго-вых процессов. Проиллюстрировано его применение на примере системы популяционной динамики типа «хищник—жертва». Детерминистические модели для таких процессов достаточно хорошо исследованы, что позволяет сравнить полученные результаты с уже известными.
Системы компьютерной алгебры очень удобны для целей быстрого прототипирования при создании и исследовании математических моделей.
Ключевые слова: стохастические дифференциальные уравнения; модель «хищник-жертва»; основное кинетическое уравнения; уравнение Фоккера-Планка; системы компьютерной алгебры; система Axiom.
Литература
1. Кулябов Д. С., Демидова А. В. Введение согласованного стохастического члена в уравнение модели роста популяций // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2012. — № 3. — С. 69-78.
2. Demidova A. V., Sevastianov L. A, Kulyabov D. S. Application of Stochastic Differencial Equations to Model Population Systems // Труды Третьей Международной конференции «Математическое моделирование социальной и экономической динамики (MMSED-2010)» / Российский государственный социальный университет. — 2010. — С. 92-94.
3. Debrabant K., Robler A. Classification of Stochastic Runge-Kutta Methods for the Weak Approximation of Stochastic Differential Equations // Mathematics and Computers in Simulation. — 2008. — Vol. 77, No 4. — Pp. 408-420. — ISSN 03784754.
4. Tocino A., Ardanuy R. Runge-Kutta Methods for Numerical Solution of Stochastic Differential Equations // Journal of Computational and Applied Mathematics. — 2002. — Vol. 138, No 2. — Pp. 219-241. — ISSN 0377-0427. — http://www.sciencedirect.com/science/article/pii/S0377042701003806.
5. Volterra V. Mathematical Theory of the Struggle for Existence. — Nauka, 1976.
6. Базыкин А. Д. Нелиненая динамика взаимодействующих популяций. — Москва-Ижевск: Институт компьютерных исследований, 2003.
7. Братусь А. С., Новожилов А. С., Платонов А. П. Динамические системы и модели биологии. — М.: Физматлит, 2010.
8. Пенроуз Р., Риндлер В. Спиноры и пространство-время. Два-спинорное исчисление и релятивистские поля. — М.: Мир, 1987. — Т. 1, 528 с.
9. Влияние стохастизации на одношаговые модели / Анастасия Вячеславовна Демидова, Мигран Нельсонович Геворкян, Александр Дмитриевич Егоров и др. // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2014. — № 1. — С. 71-85.
10. The Method of Stochastization of One-Step Processes / A. V. Demidova, A. V. Ko-rolkova, D. S. Kulyabov, L. A. Sevastianov // Mathematical Modeling and Computational Physics. — JINR, 2013. — P. 67.
11. Гардинер К. В. Стохастические методы в естественных науках. — Мир, 1986.
12. 0ksendal B. K. Stochastic Differential Equations: An Introduction with Applications. — Berlin: Springer, 2003.
13. Lotka A. J. Elements of Physical Biology. — BiblioBazaar, 2011. — ISBN 9781178508116, 492 p. — http://books.google.ru/books?id=tFN9pwAACAAJ.
14. Демидова А. В. Уравнения динамики популяций в форме стохастических дифференциальных уравнений // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2013. — № 1. — С. 67-76.
15. Демидова А. В. Метод стохастизации математических моделей на примере системы «хищник-жертва» // Научная сессия НИЯУ МИФИ-2013. — 2013. — С. 127.
16. Королькова А. В., Кулябов Д. С. Методы стохастизации математических моделей на примере пиринговых сетей // Научная сессия НИЯУ МИФИ-2013. Аннотации докладов. В 3 томах. — Москва: МИФИ, 2013. — С. 131.
17. Демидова А. В., Кулябов Д. С., Севастьянов Л. А. Согласованный стохастический член в популяционных моделях // XI Белорусская математическая конференция. — Минск: Институт математики НАН Беларуси, 2012. — С. 39.
