CC BY
68
Том 19, № 02, 2016_Научный Вестник МГТУ ГА
Vol. 19, № 02, 2016 Civil Avition High TECHNOLOGIES
A STOCHASTIC MODEL OF EPIDEMIC
N.I. OVSYANNIKOVA
First let's construct the determined model of the dynamic of uncontrolled epidemiological process. The quotient p describes frequency of meetings of sick people with healthy and probability of infection. It is subject to action of random factors. Let's enter an item for it which will take into account the influence of a random destabilization. We'll have stochastic model of epidemic. Comparing the determined and stochastic models, we'll find admissible borders for destabilizing quotient c if the maximum deviation of dynamic variables can not be higher than 5 %.
Key words: The dynamic of epidemic, the determined uncontrolled model of epidemic, stochastic model of epidemic, destabilizing quotient, perturbed coefficient.
Determined model, describing uncontrolled process of the spread of epidemic, is described by a system of differential equations:
fx(t)=-№№) - MO+Л
U (()=M0y(0- (м+М+гШ, (1)
x(t) > 0, y(t) > 0, t e [0, T], x(0) = x0, y(0) = y0, (2)
where x(t) - the rate of change in the number of people exposed to the disease, yi ((t) - the rate (speed) of change in the number of infected people,
¡x(t)y(t) - function characterizing the number of meetings of people exposed to the disease and infected ones per unit of time,
Yy(t) - the number of people who regained their health per unit of time without the influence
of external means: quarantine, vaccination and others (y~l average time of natural healing),
5 - the growth coefficient, which characterizes the frequency of meetings of healthy people with infected people (in general case it can be considered as a function ¡( x(t), y (t))), / - the coefficient of natural mortality of people, ~ - the coefficient of mortality from this infection, A - average birthrate (reproduction).
The considered mathematical model is determined and allows to calculate in advance the change of a condition of the studied system, on an interesting time segment by solving the Cauchy problem (1)-(2). We can assume that the values of some of the coefficients of the system in the moment t e [0, T] are not uniquely defined, for example, because of their dependence on many unpredictable factors, and they can be regarded as random processes, the mathematical expectations of which are known.
Assume that the coefficient of growth has a random component5, i.e. it can be represented as:
¡(t) = m(t) + a-4(t ,a>), (3)
where m(t) - mathematical expectation of the coefficient ¡, set it permanent, i.e. m(t) = 5 = const; ^t,a) - random process; a - constant characterizing the degree of influence of the random perturbation on the value of the coefficient 5.
In this case, the mathematical model (1)-(2) takes the following form:
Civil Avition High TECHNOLOGIES
Vol. 19, № 02, 2016
dx
—=-ßxy-ßx+Л- oxy^(t, а), dt
rhi
dy = ßxy - (ß+ß+Y)y+oxy£(t,a), dt
(4)
x (0,а) = x0(a), y (0,а) = y0(a).
(5)
In this case, the state of the system (x(t), y(t)) is no longer a deterministic vector-function but is a vector random (stochastic) process (x(t, (), y(t, ()), t e [0, T]. In general, the system (4)-(5) can be written:
dX (t, а) = A( X, t )dt + B( X, t) df (t, а),
(6)
X (0,а) = X>(®),
(7)
where A : R2 x [0,T] — R2; B : R2 x [0, T] — R2x1; /(t,() - scalar Wiener process;
X (t)
x(t) ^
v y(t) j
e R2, A(X, t)
- ßxy -ßx + Л ^
ßxy - (ß + ß + y) у
e R 2, B(X, t)
^ -oxy ^ axy
e R2, ^(t, а) e R1.
The obtained stochastic differential equation will be solved numerically, for this we use a stochastic analogue of Taylor's formula. Apply the unified stochastic Taylor-Ito expansion in iterated stochastic integrals and also the approximation of iterated stochastic integrals by means of the polynomial system of functions [1].
Formulate a theorem on Ito process expansion rj(s) = R(x(s), s), where R: Rn x [0, T] — R1, in
the unified Taylor-Ito series in iterated stochastic integrals Qj (s,t).
Theorem 1. Let the process rj(s) = R( x(s), s) be Ito continuously differentiable r +1 times in the mean-square sense on [0, T] along trajectories of the equation (6). Then for all s, t e [0,T], s > t it decomposes into a unified Taylor-Ito series of the following type:
n(s) = £ (cDq П(Т)} % (s © t)Dq )+ Hr+Дs, t)
q = 0
(8)
( r+1)
and there exists such a constant Cr+1 < <» that tJm{Hr+1 (s,t)2} < Cr+1 (s-1) 2 , r = 0,1,..., where
def / \
Hr+i(s,t) = (C{n(T)}Ur % (s © t)Ur )+ Dr+1(s,t),
Dr+i (s,t) = J(r {n(r)}dr % (s © t)Ar)+ )
1 л HAr {n(r)}-df(T) % (s © t) Ar
t V4
(9) (10)
Vol. 19, № 02, 2016
Civil Avition High TECHNOLOGIES
£Aq{n(T)} ={l(*C(A~Jk)e Aq} , HAq{n(T)} = {л-П(г): (Jh-Jb)e A*} , CDq {n(T)} = {(*СЛ-П(г): J.../*) e Dq } (11)
(*CЛ-.'* {.} _
(*)
L {•} =
Gk..GkL{•} if * > 0 L {•} if * = 0!
L...L« if J > 0 if J = 0'
0),
Gp{•} = I(((4L{}- (1)LGp_1 {•}), p = 1,2,...
L{} —+5 (x,t ^+1 5 £ bJ(xt b (x,t) dx/a}
G0i {•} = ±bj, (x, t) I-1; i = 1,...,
j=1 dxj
m
( * Л
Dq = U*, J, /1,..., lk): * + 2
J+X/p = q; *, J, /1,. ., /* = 0,1,. |,
p=1 ;
tfr = J(*, J, /1,..., /*): * + j + X lp * r \,
p=1
^ = |(k, j, /1,..., ik): k+j+X 1p = q; k, j, ii,..., lk = 0,1,...j,
equality (8) is just (is fair, is true, takes place, holds) with probability 1, right parts of (8)-(10) exist in the mean-square sense.
We construct a unified Taylor-Ito expansion for the components of the solution X(t) of the
5/
system (4)-(5) for infinitesimal of order of O((s -1)/2), i.e. we will construct expansion of Ito process
n(t) = X (t):
x(s) = x(t)+(s -1)[(-¡xy -/x+A)+axy(p(y - x)+(s, t)+a2 xy(p(2xy - (x - y)2) (s -1 )2
+/(x - y)K0(s, t)] +^-L(j32xy(y - x)+2pxyx -Apy-A/+(/+ft + Y)pxy)
-oxy/0(s, t)+ 0 xy(y-x) />, t)+ 03 xy(2xy-(x-y)2) C(s, t)
+оу(Л-(м+ft +Y) x)/1(s, t)+04 xy (11xy (x - y) - (x3 - y3))C0(s, t)
-O y (Л(у- x) - xy(M+ft + Y)+Mx2)/110(s, t)+О xy((M+ft + Y)(y-1)+Л)/1 (s, t) + H5x (s, t),
Civil Avition High TECHNOLOGIES
Vol. 19, № 02, 2016
y (s) = y(t) + (s -1 )[(ßxy - (ß + fi + Y) y) + oxy(ß( x - y) - (ß + fi + Y)) I0 (s, t) +o2 xy(ß(-2 xy + (x - y)2) + (ß + ß + y)( y - x)) I°0 (s, t)] (s - t)2
+
-(ß2 xy( x - y) - 2ßxy(ß + ft + Y) - ßß xy + Лßy + (ß + ft + Y)2 y)
(13)
+ oxyIf (s, t) + o2xy(x - y)I™ (s, t) + o3xy(2xy - (x - y)2)I1C1»10(s, t) + oy(ßx - Л)I1 (s, t) -o4xy (11 xy (x-y)- (x3 -y3))IZ°(s,t) + o2y (Л(y-x) -xy(ß + ft + Y) + ßx2К>,t) +o2 xy (ß (x - y) - Л) I0 (s, t) + Hy (s, t),
Relations (12)-(13) on a uniform discrete grid {Tj =0 constructed for the segment [0, T ], such that Tj = jA, tn = NA = T are selected as a numerical method for modeling the system (6)-(7). Denote x(Tj) = Xj, y(Tj) = , and then by putting s = (k + 1)A, t = kA, k = 0,1,... in expansions
(12)-(13) and using the expansions of the iterated stochastic integrals /f(s, t), //(s, t),... in terms of a polynomial basis, the following expressions for the numerical method are obtained:
Xk+1 = Xk + A[((xy - JUX + A) + oxy(p(y - x) + /1 )/0(s, t) + a2xy(p(2xy - (x - y)2)
+ß (x - y)) I00 (s, t)] +—(ß2 xy( y - x) + 2ßxy + ß2 x - Лßy -Лß + (ß + ß + y)ßxy)
-oxyI0 (s, t) + o2xy(y - x)I0 (s, t) + o3xy(2xy - (x - y)2)I01 (s, t) + oy(A - (ß + ß + Y)x)I1 (s, t) + o4xy(11xy(x - y) - (x3 - y3))I100100(s, t)
-o2y(Л(y - x) - xy(ß + ß + Y) + ßx2)I10(s, t) + o2xy((ß + ß + Y)(y -1) + Л)If? (s, t),
(14)
yk+1 = yk + A[(ßxy - (ß + ß + Y) y) + oxy(ß( x - y) - (ß + ß + Y)) If (s, t) +o2 xy(ß(-2 xy + (x - y)2) + (ß + ß + y)( y - x)) I1010 (s, t)] + ^[ß2 xy( x - y)
-2ßxy( ß + ß + Y) - ßß xy + Aßy + (ß + ß + y)2 y ] + oxyIf (s, t)
+o2xy (x - y)I"0 (s, t) + o3xy (2xy - (x - y)2)С (s, t) + oy(ßx - Л)I1 (s, t)
-o4xy(11 xy(x - y) - (x3 - y3))I10101010 (s, t)
+o2y (Л(y - x) - xy (ß + ß + Y) + ßx2)I110 (s, t) + o2xy(ß(x - y) - Л)Iff (s, t),
(15)
where
If(Tk+1,Tk ) =
A
I1(Tk+l, Tk) =
(1) 0 '
1
С(Г) +J_Ar)
I„(Tk+1,Tk) = -
-J3
=-• [(fD2 -1],
1000 (т t ) =
"411 V Lk+n k /
A 32
(i?) )3 - 3f
I (T T ) = —
-'ш^Lk+n LkJ 24
(f01) )4 - 6 (f01) )2+3
(16)
2
Vol. 19, № 02, 2016
Civil Avition High TECHNOLOGIES
I01(т T ) = _A!t1/>0K2 c(D ,c(D +_L.C(D
-41VLk+\> lk) ^^.^KbO J 2-J5
+z-
1
Z=1
•7(2/ + 1)(2i + 5) (2/ + 3)
fei fei+2
1
(2/ _ 1)(2/ + 3)
(1)\2
(Л
/10(Т т ) = с(1)Сс(1) +—С(1))_/01Сг т ) = ^^2(с(1))2 +—с(1)-с(1)
41 V 'к+1> 'k/ 2 '0 vbo ' ^ ^^ \3
1
1 9
—С(1) -с(1) + У, ,_
3^5 1=11 \l(2i +1)(2/ + 5) (2/ + 3)
C(D ,c(i) + fei fei+2 ~
1
(2/ _ 1)(2/ + 3)
(if)2 К 2]
[g\J), i = 0,1,..., q + 2; j = 1} - a system of independent Gaussian random variables with zero mean
(expectation) and variance of one, which is generated on a step of integration with a number of k and is independent with the analogous systems of random variables that are generated on all the preceding steps of integration towards (with respect to) the step of integration with the number of k; A - step of integration of the numerical method; the number q is chosen from the condition ([1], p. 199):
M {((О (Tk+1, T) _ I ¿I), (rk+1, T) )2} = M {(I1011 (rk+1, Tk) _ I (01), (rk+1, T) ) }
4 f
<
16
3 П
16 (90 5 г+ ^2 ^fri 4/2 _ 1
(17)
< C A5
where constant C must be given. We choose it for the sake of simplicity to be unity (to be equal to one). The value of q increases with a decrease in the value of the step of integration. Consider the results of the choice of number q with the help of the relation (11). These results are placed in the following table.
A 0,004 0,001 0,0005
Q 1 2 4
I.e., it is enough to let q be equal to 1 that A would be 0,004, then the expansions for I^0 and I1011 take the form:
Ow*)=- Aif^)2^ • ^+f}5 ^ • +5^(1) • ^ 4^(1))2 -2], /Ж+1,т*)=-Aif^)2С •if & •if -^^f •^э(1)+-5(£(1))2+2].
Make the numerical modeling of the solution of the system (6)-(7) by means of relations (14)-(15) on the time interval T = 10 with the step A = 0,004 with the following input (initial) data:
в = 2• 10-6; м = 0,003; ~ = 0; Л = 20; y = 1; X0(o) = 380 000; Y0(a) = 2000; o = 0. The result of the numerical modeling is presented on fig. 1. Now introduce the stochastic perturbation о > 0. The evolution of processes x(t, a>), y(t, a>), which characterize the process of evolution of epidemic of the
system (6)-(7) for values o = 10-7;2•Ю-7; 5•Ю-7; 7• 10-7; 10-6 is presented on figs. 2-6, respectively.
Civil Avition High TECHNOLOGIES
Vol. 19, № 02, 2016
The values of the maximal (maximum) trajectories deviations of the perturbed system from the trajectories of the determined system are listed in table 1, from which direct dependence of maximal (maximum) deviations of the solution of a perturbed system on the value of the perturbed parameter a is well seen.
2 000 1 SOD 1600 1400 1200' 1000 800 600 400 200
It)
E а
О
i
{N.
[
1
0 1 2 3 4 . S 6 7 6 9
Fig. 1. Determined model of epidemic, a = 0
2 КО I«« "ЙЮ 14(0 1JW ISM 300
V.........................................
Л V
1 \ч
f '4N. j
N4 ; !
\>Ч 1
—Г" t weeks-—;*....... ;
Fig. 3. 7 = 1,5 10
-7
2НЮ-1SOO-1600140012001 ЛИ-S00' 600400-
1 ! ! ! ! !
ü }:
1 ] ! ! ! ! Г
--------:--------r --X weeks........^ТГГ
0 1 2 3 4 5 6 7
Fig. 2. Determined and stochastic models of epidemic, a = 10
-7
2000 1IM ISM HDD 1Й0 1000 ООО S» IM MO
Fig. 4. 7 = 2 10
%
\ ..........
1 л -
I V
К % j
l, tttdu "" ,
Fig. 5. 7 = 2,5 10
-7
3 4 S
Fig. 6. 7 = 3 10
6 7
-7
Table 1
7
Value of 0 Maximum deviation of the trajectories of the perturbed system from the trajectories of the determined system depending on the value of the parameter a
Deviation X, % Deviation Y, %
10-7 Practically none 2,5
1,5 10-7 Practically none 3,5
2 10-7 0,1 4
2,5 10-7 0,12 4,5
3 10-7 0,15 5
Vol. 19, № 02, 2016
Civil Avition High TECHNOLOGIES
For a< 10 8 stochastic model is almost equal to the determined one, so it is necessary to take the determined model for description of the system; for a> 3 10-7 the stochastic model significantly (by more than 5%) is different from the determined one; therefore the perturbed coefficient needs to be taken ranging from 10-8 to 3 • 10-7.
For different realizations of the system of independent Gaussian values ), i = 0,1,...,q + 2; j = 1,2} we obtain different realizations of the solution of the system of stochastic
differential equations (6). These trajectories for small perturbations lie inside of a tube constructed in a small neighbourhood of the solution of determined system (1). Find mean of the solution of system (6) in 5, 10 realizations for a = 5 • 10-7, in 5, 10, 15 and 50 realizations for a = 10-6.
On figs. 7-12 it is shown a comparison of means with the solution of the determined system of differential equations (1).
One can conclude that M {X ((;«)} ^ X (t), which is confirmed by numerical experiments, where X(t) = (x(t),y(t)) - solution of determined system (1), (2), X(t;m) = (x(t,m),y(t,m)) - solution of stochastic system (4), (5).
23456789 Fig. 7. Determined model and mean of the solution
_7
of system (6, 7) found in 5 realizations (a =
4M'
400
; ; ;
и
} и n..........
s V ! ;
>
1 1
------- ц wcols , ,
0(21456789 Fig. 8. Determined model and mean of the solution
_7
of system (6, 7) found in 10 realizations (a =
200
и yV 7!. \ M
1 \y (L .......
8 i
>
\ V
\
\ t: weeks "W
0 1 2 3 4 5 6 7 3 s Fig. 9. Determined model and mean of the solution
of system (6, 7) found in 5 realizations (a =
\ 1
¡j \ }
а rv\
8
1 1 | Г1^4-^-.
1 , Г Ц WCCUb , ,
0 1 2 3 4 5 6 7 3 9 Fig. 10. Determined model and mean of the solution
of system (6, 7) found in 10 realizations (a =
Civil Avition High TECHNOLOGIES
Vol. 19, № 02, 2016
Fig. 11. Determined model and mean of the solution Fig. 12. Determined model and mean of the solution
of system (6, 7) found in 15 realizations (О = 10-6) of system (6, 7) found in 50 realizations (о = 10-6)
REFERENCES
1. Kuznetsov D.F. Numerical modeling of stochastic differential equations and stochastic integrals St. Petersburg: Science. 1999. 459 p. (in Russian).
2. Dmitrieva O.N. A stochastic model of the dynamics of the forests. Collection of proceedings. Tver. 2006. 187 p. (in Russian).
СТОХАСТИЧЕСКАЯ МОДЕЛЬ ДИНАМИКИ ЭПИДЕМИИ
Овсянникова Н.И.
Строится детерминированная модель динамики неуправляемого эпидемиологического процесса. Затем, считая, что коэффициент в, характеризующий частоту встреч и вероятность заражения при встрече, подвержен воздействию случайных факторов, введем для него слагаемое, учитывающее влияние случайного возмущения. Получим стохастическую модель эпидемии. Сравнивая детерминированную и стохастическую модели, найдем допустимые границы для возмущенного коэффициента О при условии, что максимальное отклонение динамических переменных не должно превышать 5 %.
Ключевые слова: Динамика эпидемии, детерминированная неконтролируемая модель эпидемии, стохастическая модель эпидемии, дестабилизирующий фактор, возмущенный коэффициент.
СВЕДЕНИЯ ОБ АВТОРЕ
Овсянникова Наталья Игоревна, окончила Поморский государственный университет имени М.В. Ломоносова (1990), кандидат физико-математических наук, доцент кафедры прикладной математики МГТУ ГА, автор 30 научных статей, монографии, сфера научных интересов - математическое моделирование сложных динамических систем, оптимальное управление, численные методы, вариационное исчисление, обыкновенные и стохастические дифференциальные уравнения, математическая статистика, электронный адрес: natmat68@mail.ru.
