CC BY
42
МАТЕМАТИКА
УДК 531.19+519.21
А. В. Калинкин, А. М. Л а н г е,
А. В. Мастихин, А. А. Шапошников
ЧИСЛЕННЫЕ МЕТОДЫ МОНТЕ-КАРЛО ДЛЯ МОДЕЛИРОВАНИЯ СХЕМ ВЗАИМОДЕЙСТВИЙ ПРИ ДИСКРЕТНЫХ СОСТОЯНИЯХ
Рассмотрены стохастические модели систем взаимодействующих элементов, задаваемые кинетическими схемами. Представлен метод Монте-Карло численного исследования таких систем. В качестве примеров взяты модели бимолекулярной реакции, распостра-нения эпидемии и брюсселятора в виде марковских случайных процессов с дискретным фазовым пространством.
Сложные системы с взаимодействиями и превращениями составляющих их элементов часто задаются схемами взаимодействий. Общая схема взаимодействий, в которой участвуют элементы типов Т,..., Тп, имеет вид
e]Tl + 4T2 + ... + 4Tn ^ 7\Ti + y\T2 + ... + 7ПTn,
elT + 4T2 + ... + 4Tn ^ 7IT + y2T2 + ... + 7nTn, (1)
е!т + 4Т2 + ... + £1гТп ^ т!Т + Т2 + ... + тПТп,
где е], 7], г = 1, 2,..., /, ] = 1, 2,..., п, — целые неотрицательные числа. В физике и химии схемы вида (1) называются кинетическими схемами [4, 8]; в экологии и эпидемиологии при описании систем с взаимодействием используется специальная терминология, однако во многих случаях такие описания сводятся к схемам вида (1) [11, 12]. В технике, в математической теории надежности запись Т ^ 0 обозначает выход из строя элементов технических систем; в теории массового обслуживания схема 0 ^ Т, Т ^ 0 интерпретируется как дисциплина обслуживания М/М/то [17]. В публикациях по математическому мо-делировани систем рассматриваемого класса для отдельных элементов используется название частица, что обусловлено спецификой процессов их взаимодействия. ри различии физической природ , происхождения и масштабов таких частиц существенным является то, что
для исследования этих систем используются единые математические модели.
В первой части работы дается описание детерминированных и стохастических моделей, используемых при рассмотрении систем с взаимодействием, задаваемых схемами (1). При детерминированном подходе к схеме (1) определяют систему нелинейных дифференциальных уравнений (см. далее уравнения (6)). Вероятностный подход основан на понятиях и результатах теории марковских процессов со счетным множеством состояний [2]. Моделью системы со схемой взаимодействий (1) является обобщенный процесс рождения и гибели £(£), £ € [0, то), при дискретном множестве состояний Мп, N = {0,1, 2,...} [4, 5, 17]. Точка фазового пространства («1 ,...,ап) € Nn интерпретируется как такое состояние систем , в котором имеется совокупность частиц а1Т1 + ... + апТп, состоящая из а1 частиц типа Т1, ..., ап частиц типа Тп; переход случайного процесса в другое состояние — результат взаимодействия одного из комплексов частиц £г1Т1 + ... + егпТп согласно схеме (1). Процесс £(£) введен в работе [1] для кинетической схемы Т + Т^ ^ Тк + Т, i,j,k,l = 1,..., п. В работе [1] указано на связь между стохастической и детерминированной моделями при большом числе частиц (см. таюке [4, 5, 12]). В настоящей статье примеры "термодинамического предельного перехода" даны для трех схем взаимодействий.
налитический метод исследования марковских процессов основан на рассмотрении первой и второй систем дифференциальных уравнений Колмогорова для переходных вероятностей (см. далее уравнения (7)). Для процессов £(£), определяемых схемами взаимодействий (1), такое исследование сводится к рассмотрению уравнений в частных производных (9) и (10). Точные решения уравнений Колмогорова для некоторых схем взаимодействий, в случаях п < 2 и I < 2, получены в работе [17].
В случаях п > 2 или I > 2 ввиду сложности уравнений (9) и (10) целесообразно использовать метод Монте-Карло — способ решения математических задач путем моделирования случайных величин и построения статистических оценок [3]. Модели при дискретных состояниях, совпадающие с рассматриваемыми в настоящей статье, изучались через численный эксперимент в работах: [9], где детально исследована схема 2Т ^ кТ, к = 0,1; [10], где моделировалась химическая реакция 2Т ^ Т2, Т2 ^ 2ТЬ Т ^ 0, Т2 ^ Т3; [15], где методом статистических испытаний изучалась динамика конкурирующих популяций; [20] и др.
Во второй части статьи изло ен итерационный алгоритм моделирования на ЭВМ случайного процесса £(£), соответствующего схе-
ме (1). етодом статистических испытаний рассмотрен предельные поведения при £ ^ то стохастических систем с конкретными схемами взаимодействий: бимолекулярная схема 2Т ^ Т, 3Т; Т ^ 0, 2Т; 0 ^ Т; процесс эпидемии Т1 + Т2 ^ кТ1, к = 1, 2; Т1 ^ 0; брюс-селятор 2Т1 + Т2 ^ 3Т1; Т1 ^ 0,Т2; 0 ^ Т1. Реализация моделей схемы (1) в виде комплекса программ для ЭВМ позволяет проводить вычислительные эксперименты и получать соответствующие выводы об объекте моделирования.
Детерминированные модели для схем взаимодействий. Кинетические уравнения. В химической кинетике для мономолекулярной реакции используется запись Т1 ^ Т2. При детерминированном подходе реакция описывается количеством х1(£) реагента Т1 и количеством ж2(£) реагента Т2 в момент времени £ € [0, то). Полагают справедливым феноменологический закон [8]
ж1 = —Лж1, ж2 = Лж1, х1(0) = ж1, ж2(0) = х2, (2)
где Л > 0 — константа скорости реакции.
Для бимолекулярной реакции Т1 + Т2 ^ Т3 пусть х1(£), ж2(£), ж3(£) — количества реагентов типа Т1, Т2, Т3. Полагают справедливым закон действующих масс [8]:
ж1 = —Лж1ж2, ж2 = —Лж1ж2, ж3 = Лж1ж2,
х1(0) = ж0, ж2(0) = Х2, х3(0) = ж3. (3)
ядерной физике для цепной реакции размно ения нейтронов используется запись Т ^ кТ, к = 0,1, 2,... .Уравнение детерминированной модели имеет вид [4, 5]
х = Лж, ж(0) = ж0, (4)
где —то < Л < то.
Динамику экологической системы "хищник-жертва" описывают количеством ж1(£) "жертв" и количеством ж2(£) "хищников". Схеме взаимодействий Т1 + Т2 ^ кТ2, к = 0,2; Т ^ 2Т1; Т2 ^ 0 [4, 11, 20] соответствует система дифференциальных уравнений
ж1 = —Л1ж1ж2 + Л2ж1, ж2 = Л3ж1 ж2 — Л4ж2, ж1(0) = ж0, ж2(0) = ж2,
(5)
где Л1 > 0, Л2 > 0, Л3 > 0, Л4 > 0.
При детерминированном подходе, в общем случае схемы (1), вводят количество жД£) частиц типа Т^, i = 1,..., п. Функции ж1(£),..., жп(£)
удовлетворяют системе нелинейных дифференциальных уравнений
X 1 = /1 (х 1, . . . ,Хп), .................. (6)
хп /п(х^ . . . , хп)
с начальными условиями х 1 (0) = х1 , ..., хп(0) = хп (кинетические уравнения). Вид функций /1,..., /п определяется по схеме (1) исходя из законов формальной кинетики [8] на основе простейших случаев (2) -(5). В прикладных задачах, как правило, функции /1,..., /п являются многочленами степени не выше третьей. Схемой (1) может быть учтено поступление частиц извне (открытая система) и образование финального продукта.
Стохастические модели для схем взаимодействий. Уравнения Колмогорова. Пусть = {а = (а1,...,ап), а = 0,1,2,..., г = 1,..., п} — множество п-мерных векторов с неотрицательными целочисленными компонентами. Далее для а, в, 7 £ приняты обозначения: 7 = а — в, если 71 = а1 — въ ..., 7п = ап — вп, и т. п. Для однородного по времени марковского процесса£(£) = (^(£),..., £п(£)), £ € [0, то), на множестве состояний обозначим переходные вероятности Рав(Ь) = Р{£(£) = в | ¿(0) = а}, а, в € Такой случайный процесс задается плотностями переходных вероятностей = (¿)/^£)^=0+. Далее всегда выполнены условия, при ко-
торых переходне вероятности удовлетворя т первой (обратной) и второй (прямой) системам дифференциальных уравнений Колмогорова [2]
= £ а„7Rlß(t), а g Nn,
Y€N n
¿Рав (t)
(7)
dt
Y€N'
(t) а7в, в G Nn,
с начальными условиями Paa(0) = 1, Рав(0) = 0 при а = в.
Марковский процесс, соответствующий схеме взаимодействий (1), указывается конкретным видом плотностей переходных вероятностей |аав,а,в G Nn} [17]. Из схемы (1) имеем векторы ег = (el,..., еП), i = 1,..., /. Каждому вектору ег сопоставим распределение вероятностей на Nn, {р^ > 0, 7gN„ = 1, Ре* = 0}, и набор чисел l^a > 0, а G Nn; ^ = 0, если ак < ek при некотором k} по формуле
n
= A П aj(а,- - 1)... (а, - ej + 1), а G Nn, (8)
где Xi > 0, i = 1,..., n, — коэффициенты пропорциональности. Для марковского процесса £(t) по определению полагаем
= - ^ a*ß = Y1 ^aPß-a+eJ, а = в ^ в ^ N'
Процесс £(£) интерпретируется как стохастическая модель системы взаимодействующих частиц п типов Т,..., Тп. Событие (£(£) = в} есть такое состояние системы, в котором в момент времени £ имеется совокупность Бв частиц, состоящая из в1 частиц типа Т,..., @п частиц типа Тп: Бв = в1Т1 + ... + впТп. Зададим I комплексов взаимодействия частиц Бег, соответствующих векторам ег. Через случайное время т^,
Р(т^ < £} = 1 — происходит взаимодействие комплекса частиц
Бег. В этот момент из в1 частиц типа Т выбирается частиц, ..., из вп частиц типа Тп выбирается егп частиц и этот комплекс частиц Бег с распределением вероятностей (р^} заменяется совокупностью Б7 новых частиц. Система из состояния Бв, соответствующего вектору в, переходит в состояние Бв-ег+7, соответствующее вектору в — + 7 (см. рис. 1), и далее происходит аналогичная эволюция системы частиц. В состоянии Бв система находится случайное время те, пока не произойдет какое-либо одно из I взаимодействий, т.е. тв = тт(т!,..., т^). Предполагается, что случайные величины т|,..., т^ независимы. Тогда
Р(тв < £} = 1 — + ^ и вероятность события, что произошло взаимодействие комплекса частиц Бег при условии, что взаимодействие произошло, равна (X]!=1 ^в)-1. Возможные превращения частиц в такой системе представляются схемой взаимодействий (1), где случайный вектор 7г = (71,..., 7п) имеет распределение (р^}, г = 1,..., I.
Выбор значений (8) для ^в объясняется следующим образом. Пусть марковский процесс находится в состоянии в = (вь ... ,вп), что соответствует наличию совокупности частиц Бв = в1 Т1 + ... + впТп.
Предполагается, что за время ДЬ, ДЬ ^ 0, вероятность ^ДЬ + о(ДЬ)
£г
взаимодействия комплекса частиц £е» пропорциональна числу С^ сочетаний частиц типа Т1 из имеющихся в1 частиц типа Т1, ..., пропорциональна числу Сд" сочетаний ^п частиц типа Тп из имеющихся вп частиц типа Тп.
Для введенных марковских процессов уравнения Колмогорова (7) для переходных вероятностей записываются в компактном виде с использованием производящих функций. Многомерной производящей функцией ^ (й1, ..., зп), соответствующей целочисленному случайному вектору £ = (^1,..., £п) с распределением вероятностей {ра > 0, а € ^„ ра = 1}, называется функция
F (si,...,sn) = Msl
ix
Q^n _
n
E
aeN n
0ax
paQl . . . Qn
(см. [2]). Далее применяется запись 8 = (в1,..., зп), = з^1 ... . Для вектора используется обозначение 1, если все его компоненты равны единице; 181 обозначает вектор с компонентами | в» |. Математические ожидания компонент случайного вектора вычисляются по формуле
M6 =
df (q)
dQi
s=1
дисперсия — по формуле
D6 =
д (q)
ds2
+
df (q)
s=1
dQi
s=1
df (q)
dsi
s=1
ля свертки второй системы дифференциальных уравнений используются производящие функции (|з| < 1)
Fa(t; q) = M(Q«(i) | ^(0) = а) = E Рав(t)Qe, hi(s)
eeN n
YeN n
Теорема 1 [17]. Производящая функция переходных вероятностей в) марковского процесса £(£), соответствующего схеме взаимодействий (1),удовлетворяет при |з| < 1 линейному дифференциальному уравнению в частных производных
dfa(t; q) dt
Q) - Qe*) д£Fa(t; Q), Fa(0; q) = (9)
i=1
ds'
£
где Aj > 0, i = 1,..., n.
n
2
Для свертки первой системы дифференциальных уравнений вводятся экспоненциальные производящие функции и дифференциальные операторы
za / д \ д71+...+Yn
G (t; z) = E 0! 4 d*) = ^ pY dz/1 ...dznn ,
a€Nn ~f£Nn 1
где z = (zi,..., zn), za = z^1 ... z^n, a! = ai!... anl.
Теорема 2 [17]. Экспоненциальная производящая функция переходных вероятностей Ge (t; z) марковского процесса £(t), соответствую -щего схеме взаимодействий (1), удовлетворяет линейному дифференциальному уравнению в частных производных
^ = t ^ (Ч|) - £)Ge(t; z), Ge(0; z) = в. (10)
i=i
Функции Fa(t; s), h^(s) и G^(t; z) являются аналитическими в рассматриваемых областях.
Алгоритм моделирования на ЭВМ марковского процесса ^(t).
При рассмотрении детерминированной модели схемы взаимодействий (1) — нелинейной системы уравнений (6) для x1(t),..., xn(t) — применяется метод Рунге-Кутта построения приближенного решения системы обыкновенных дифференциальных уравнений.
ри исследовании вероятностной модели схемы взаимодействий (1) строятся реализации процесса £(t) = (&(t),..., £n(t)) методом Монте-К арло .Алгоритм моделирования следующий.
На пользовательский интерфейс вводятся исходные данные: [0, T] — время моделирования; а = (ai,...,an) — начальное состояние; е1, ..., £1 — комплексы взаимодействия; {p,},..., {py} — распределения вероятностей на Nn; Л1,..., Л; — набор чисел.
Используются переменные: m — номер итерации алгоритма; t(m) — m-й момент изменения состояния марковского процесса; p(m) = = (pi ',..., pn ) — состояние процесса на m-й итерации; т = = (т1,..., т1) — вектор случайных времен, ттщ — минимальная координата вектора т; y = (y1, ..., Yn) — вектор скачка, координаты которого вычисляются на m-й итерации; r = random[0,1] — датчик случайных чисел, равномерно распределенных на отрезке [0,1].
1. Начальные значения переменных: m = 0, t(0) = 0, Р(0) = а.
2. Реализация вектора т. Цикл по i = 1,..., / .
2.1. Если pkm) < £k при некотором k, k = 1,... ,n, то тг = то. Иначе r := random [0,1], вычислить ^(m) по формуле (8) и присвоить значение т = -(1/^(m))ln(1 - r).
3. Определить rmin = min(r1,..., т1). Пусть rmin = т¿. Если rmin = то, то конец алгоритма.
4. Задать r := random[0,1]. С помощью r реализовать случайный вектор y с распределением (p^,}.
5. Присвоить значения переменным: t(m+1) = t(m) + rmin, в(т+1) = = в(m) - ег + Y.
6. Если t(m+1) < T, то m := m +1 и переход к п. 2. Иначе конец алгоритма.
По полученному массиву значений переменных (t(m),e(m), m = = 0,1, 2,...} производится визуализация числовых данных:
— вывод графиков траекторий детерминированной модели ^¿(t) и реализаций стохастической модели <^(t) в зависимости от t, t Е [0,T], i = 1,... ,n;
— вывод графиков траекторий детерминированной модели на фазовой плоскости x^Oxj, i, j = 1,..., n, i = j;
— вывод графиков реализаций £»(t), (t) стохастической модели на фазовой плоскости N2, i, j = 1,..., n, i = j.
Замечание 1. В п. 2.1 алгоритма используется следующее: если случайная величина r равномерно распределена на отрезке [0,1], то случайная величина, определенная формулой т = —(1/Л) ln(1 — r), имеет показательное распределение P(r < t} = 1 — e-Ai.
Замечание 2. В п. 3 алгоритма случай rmin = то означает остановку случайного процесса в одном из поглощающих состояний.
Вычислительные эксперименты проводились в системах Matlab, Maple, Delphi и других с использованием стандартных датчиков случайных чисел. Пример копии экрана ЭВМ с интерфейсом для одной из программ представлен на рис. 7 (см. далее).
Бимолекулярная схема. Стационарное распределение. На множестве состояний N = (0,1, 2,...} рассматривается марковский процесс £(t), t Е [0, то), соответствующий схеме взаимодействий
2T ^ T, 3T; T ^ 0, 2T; 0 ^ T. (11)
Вторая система дифференциальных уравнений Колмогорова для переходных вероятностей при помощи производящей функции Fi(t; s) =
= Pj (t)sj, i Е N, |s| < 1, записывается в виде уравнения в частных
производных второго порядка:
s) w 2 з , 2 2n d2Fi(t; s)
""дТ" = A2 2s + pis -
s)
+ Ai(p1s2 + pi - s) ds 7 + Ao(s - 1)Fi(t; s), Fi(0, s) = s\ (12)
Р0 = 1 Яг Pi
О 1 2 i - 1 i i +1
Рис. 2. Скачки марковского процесса, соответствующего схеме (11)
где Л2 > 0, А1 > 0, А0 > 0, рЗ > 0, р2 > 0, рЗ + р2 = 1, р2 > 0, Ро > 0,
р2 + Ро = 1-
Рассматриваемый случайный процесс является процессом рождения и гибели квадратичного типа [2]. В состоянии г процесс находится случайное время т, Р{т < Ь} = 1 — е-(Ао+А1г+А2Ф-1))^ переходы процесса г ^ г — 1 или г ^ г + 1 происходят с вероятностями (см. рис. 2)
~ = рОЛог + рОЛ2г(г — 1) _ = Ао + р-^г + р2Л2г(г — 1) (13) г Ао + Лог + Л2г(г — 1)' г Ло + Лог + Л2г(г — 1)
Имеем детерминированную модель для схемы взаимодействий (11). Дифференцируя уравнение (12) по 8, получим
d2F(t; s) = dsdt
= A2(3p3s2 + - 2s)^^ + A2(p3s3 + p2s - s2)
d2F(t; s) . 2 3 2 2ч д3F(t; s) 3S + p - - + A2(P2s + Pls - s )—p--
+ Ло(2р28 — 1) ^^ + Л^2 + рО — 8) ^^^ +
+ ЛоР(Ь; 8) + Ло(8 — 1)^^^. (14)
Вводя обозначение для среднего числа частиц А(Ь) = (др(Ь; 8)/д8)|в=1 и учитывая, что 1) = 1, получаем из (14) равенство
^ = А2(2РЗ — 1) ^^^ + Л1(2Р2 — 1)А(*) + Ло.
Считая при г ^ то справедливым "предельный термодинамический переход" [1], [5]
д2Fi(t; s) / dF(t; s)
ds2
s=!
2
в=1 V
приходим к кинетическому уравнени
Х = Л2(2р3 — 1)х2 + Л1 (2р2 — 1)х + Ло, х(0) = хо, (15)
где х (Ь) — количество реагента в момент времени Ь для бимолекулярной реакции со схемой (11). В случае р2 < 1/2 уравнение (15) известно
Рис. 3. Детерминированная модель х(£) и стохастическая реализация £(£)
как уравнение популяционной динамики с внутривидовой конкуренцией [11].
На рис. 3 дана траектория детерминированной модели х(£) при начальном условии х(0) = 10 и приведен пример реализации стохастического процесса £(£) при начальном условии £(0) = 10. Значения параметров составляют: Л2 = 2, А1 = 4, 5, Л0 = 4, 5; р3 = 1/4, Р? = 3/4, = 2/3, р0 = 1/3; Т = 3,0688.
Схема (11) интерпретируется как система попарно взаимодействующих частиц с притоком частиц извне [4, 5]. При р3 < 1/2 вследствие формул (13) можно показать [2, гл. 7, § 4], что в системе при £ ^ то существует стационарное распределение, определяемое распределением вероятностей ^,] = 0,1, 2,...}, где ^ = Р^(£). Выражения
для стационарных вероятностей ^ известны [2], но малопригодны для асимптотического исследования. Уравнение для производящей функции стационарных вероятностей /(в) = ^^ qjв, |в| < 1, имеет вид
¿=0
л?(рЗв3 + Р1в - в2)+
+ Л1 (Р?в2 + р0 - в)/в" + Ао(в - 1)/(в) = 0. (16)
Уравнение (16) сводится к гипергеометрическому уравнению, однако исследовать свойства стационарного распределения, рассматривая функцию /(в), не удается [19]. Уравнение (16) при А1 = 0 и р = 1 является частным случаем рассмотренного в работе [19] более общего
уравнения для стационарных вероятностей; установлена асимптотическая нормальность стационарного распределения {qj } при большой интенсивности Ао поступления новых частиц.
Теорема 3 [19]. Рассмотрим схему взаимодействий 2T ^ Т, 0 ^ Т. Обозначим п = П (А0) целочисленную случайную величину с распределеним {qj, j = 0,1, 2,...}. Положим x G (-те, то). Тогда
lim P{ П(Ао) -C1 Ао < x} = J= Г e-У2/2 dy, Ло ^ 1 с2л/ А0 J V2W
где ci > 0, c2 > 0 — некоторые константы.
Для схемы взаимодействий (11) описанным выше методом статистического моделирования получены оценки вероятностей стационарного распределения qj, j G N, вычисленные как отношение суммарного времени нахождения процесса в состоянии j ко всему времени моде-
m
лирования Т: qj = ^^ rj/Т, rj — время нахождения процесса в состо-
1=1 ю
янии j при Z-м попадании в это состояние qj = 1. На рис. 4 при-
j=0
ведена гистограмма для {qj} при следующих значениях параметров: А2 = 2, А1 = 15, Ао = 50; p2 = 1/4, р2 = 3/4, p2 = 2/3, pi = 1/3; Т = 25348. На рис. 5 дана гистограмма для {qj} при значениях параметров А2 = 2, А1 = 15, А0 = 500000; p2 = 1/4, p2 = 3/4, p2 = 2/3, pi = 1/3; Т = 6, 5953.
Сравнивая гистограммы, приведенные на рис. 4 и рис. 5, видим, что при увеличении параметра А0 они близки к нормальной плотности распределения; выборочные средние и дисперсия — параметры кривой нормального закона на рис. 5 — вычислены по стандартным формулам. Аналогичные численные результаты получены при других значе-
Ч) 0,10 0,08 0,06 о,оь 0,02
Рис. 4. Гистограмма стационарного Рис. 5. Гистограмма стационарного распределения при Ло = 50 распределения при Ло = 500 000
ниях параметров. Таким образом, статистическое моделирование марковского процесса £(£) со схемой взаимодействий (11) позволяет сделать предположение об асимптотической нормальности стационарного распределения {д^} при Ао ^ то.
Проведенные методом Монте-Карло исследования моделей схем взаимодействий общего вида еТ ^ кеТ, ..., 2Т ^ к2Т, Т ^ к1Т, 0 ^ коТ таюке привели к предположению об асимптотической нормальности стационарного распределения при Ао ^ то и некоторых дополнительных условиях.
Процесс эпидемии. Финальное распределение. На множестве состояний N рассматривается марковский процесс £(£) = (^(¿), £2(£)), £ € [0, то), соответствующий схеме взаимодействий
Т1 + Т2 ^ кТ1, к = 1, 2; Т1 ^ 0. (17)
Первое уравнение для производящей функции переходных вероятностей С(в1,в2)(Ь; ¿2) в случае процесса при , > 0 имеет вид
dGß (t; z) dt
д2 д д2 ziz2( P2^ + Piö--о я ) +
+ 1 - д^) Gß(t; z)
с начальным условием Се(0; г) = /в!. Здесь р1 > 0, р2 > 0, р1 + +р2 = 1. Второе уравнение для производящей функции Р(а1,а2)(£; 81, 82) получает вид
дМ; 8) , 2 , д2Ра(¿; 8) . , дРа(Ь; 8) /1ЛЧ
с начальным условием Ра(0; в) =
Рассматриваемый случайный процесс является двухмерным процессом рождения и гибели квадратичного типа; пример реализации процесса изображен на рис.6, случай а. Событие {£(£) = (във2)} означает наличие в1 частиц типа Т1 и в2 частиц типа Т2 в момент времени ¿.В вероятностных моделях распространения эпидемий частицы типа Т1 интерпретиру тся как больные особи, частиц типа Т2 — как здоровые особи, восприимчивые к инфекционному заболеванию. Через случайное время т2^ ^), Р{т(в1 ^) < ¿} = 1 — е-в1в24, происходит контакт больной и здоровой особей. Здоровая особь становится больной и с вероятностью р2 остается в популяции, и процесс переходит в состояние, соответствующее вектору (в1 + 1,в2 — 1), или с вероятностью р1 удаляется из популяции, и процесс переходит в состояние,
соответствующее вектору (въвг — 1)-Кроме того, через случайное время
Чвъв2 )
, P{
Т
(в1 ,в2)
< t} = 1 - e-Meit,
больная особь умирает и процесс переходит в состояние, соответствующее вектору (в — 1, ) - Случай-
1 2
ные величины T(iвьв2), Т^ независимы; в состоянии (във2) процесс находится случайное время
т(в1,в2) = ттСгв.в,)))- Случай
(0,72)
£*1
Рис. 6. Скачки марковских процессов на N2
(в!,в2)' '(в1,в2) р = 1 называется моделью эпидемии
Вейса [13], а случай р2 = 1 называется
моделью эпидемии Бартлетта-Мак-Кендрика [14].
Чтобы получить детерминированное приближение для процесса
£(£), соответствующего схеме взаимодействий (17), дифференцируем
второе уравнение (18) по и 82; получаем
д2Fa(t; s)
ÖSlÖt
, ч д2 Fa (t; s)
(2P2S1 + P1 - S2) . a) +
+ (P2S1 + PlSi - SiS2)
д3Fa (t; s) ds1ds2
ds1ds2
dFa(t; s) , . ч д2Fa(t; s) ß—^--Ь 1 - sij
дs1
дs2
д2Fa(t; s) д2Fa (t; s)
= s1 ——---h
дs2дt
дs1дs2
+ , 2 + 4 д 3Fa(t; s) + , 4 д 2Fa(t; s)
+ (p2 s2 + p1s1 - +ß (1 - s1) ^sr
Подставив 8 = 1 и используя обозначения для среднего числа частиц типа Т, А(£) = (д^Ц^; )|5=1, г = 1, 2, получаем равенства
dA1(t) _ д2Fa(t; s) = p2
dt дs1дs2
dA2(t) д2Fa(t; s)
dt
дs1дs2
s=1
s=1
- ßA (t),
(19)
Предположим, что при большом начальном числе частиц а = = (па!, па2), п ^ то, выполнено условие "предельного перехода" [1]
д2Fa(t; s)
дs1 дs2
s=1
дFa(t; s) дFa (t; s)
s=1
дs1
дs2
s=1
тогда из (19) приходим к системе дифференциальных уравнений
xX 1 = p2x1 x2 - ßx1, xX2 = -x1x2, x1(0) = x0, x2(0) = x0,
£ Схема взаимодействий: T1+Т2 —> кТ1, к=1,2 ; Т11л> О
0111®
jf|_ ВЫХОД График будет сохранен в Файле |
[0; Т] П шаг параметр д параметр р2 Х1(0), Кси1 (0) Х2(0), Кси2(0) количество реализаций
Вывод графиков
Г ХЩХЗД
С фазовая плоскость Х1, Х2
Г стохастическая реализация Кси1 (1), Кси2(1)
С фазовая плоскость Кси1, Кси2
* совместно Х1, Х2; Кси1, Кси2 С совместная Фазовая плоскость Г финальные вероятности
Пуск
Очистить
Вернуть исходные данные
Сохранить
Рис.7. Детерминированная модель и стохастические реализации
ш, т
где х1(Ь) — количество больных особей, х2 (¿) — количество особей, восприимчивых к инфекционному заболевани в детерминированной модели эпидемии. заимосвязь вероятностного и детерминистического описаний для различных схем эпидемий рассматривалась в [12, 13] и др.
На рис. 7 дана траектория детерминированной модели х1(Ь), х2(Ь) при начальных условиях х1(0) = 10, х2(0) = 100, и приведен пример реализации стохастического процесса эпидемии (^(¿), £2(Ь)) при начальных условиях (^1(0),^2(0)) = (10,100). Значения параметров следующие: , = 3, р2 = 0,8, Т =1.
С вероятностью, равной 1, процесс остановится в одном из поглощающих состояний (0,72), 72 = 0,1, 2,... (остались только здоровые особи). Для финальных вероятностей д^.О)2'' = ро^г)^^) вводим производящие функции (181 < 1)
/(«1,«2)(8) = Е ^Г)2^72, ^(г1,г2; 8)= Е Л«1>«2)(8).
72=о а1,а2=о
Стационарное первое уравнение Колмогорова получает вид
( д2 д д2 \ ( д\п
4P2äü+P1 äzi- д^) + н1 - g(zi'z2; s) =0 (20)
с граничными условиями g(0,z2; s) = ez2s, g(z1, 0; s) = ez1. Рвение методом Римана задачи Дарбу-Пикара для гиперболического уравнения (20) дает в случае p1 = 1 [16] при а1 =0 выражение для производящей функции
1 Гх
f(ai,a2)(s) = (а _ 1)J ха1-1(1 - + se"x/^)a2 e-x dx. (21)
Пусть n(ai ,a2) — число финальных частиц типа T2, которые останутся после остановки процесса. Производящая функция (21) задает распределение случайной величины n(ai,a2) на состояниях {(0,y2), y2 = = 0,..., а2}. Из (21) получаем для среднего и дисперсии: Mn(ai,a2) = = c1a2, Dn(ai,a2) ~ c2a| при а2 ^ то, где c1 > 0, c2 > 0 — некоторые константа. Случай а2 ^ то представляет интерес для приложений: в начальный момент времени имелось большое число здоровых особей и некоторое число больных особей. Из (21) стандартным методом характеристических функций устанавливается
Теорема 4 [16]. Рассмотрим марковский процесс, соответствую -щий схеме взаимодействий T1 + T2 ^ T1, T1 ^ 0. Положим x G [0,1]. Тогда
( n(«i,«2) 4 lim Р\ n-< А = F(ai)(x), (22)
«2^^ I «2 J
ai —1
где F(ai)(x) = x^ ^ (_^ ln x)n/n!.
n=0
Предельные соотношения типа (22) называют "пороговыми теоремами", они используются для установления численного значения так называемого "эпидемического порога" при распостранении инфекционных заболеваний [12].
Методом Монте-Карло для схемы взаимодействий (17) построены эмпирические функции распределения Fiai,a2)(x) по выборкам, состоящим из n реализаций случайной величины n(ai,a2)/а2 (линейная по а2 нормировка). Статистические испытания показали сходимость функции F(ai,a2)(x) к некоторой функции распределения при увеличении n.
На рис. 8, случай 1, приведен график функции F(3)(x) = x^(1 — —^ ln x + ln2 x/2) при ^ = 40. При значениях параметров ^ = 40, p2 = 0 и начальных условиях a1 = 3, a2 = 200 было проведено n = 300 испытаний. Эмпирическая функция распределения F3(0o200)(x), представленная на рис. 8, случай 2, близка, в соответствии с теоремой 4, к функции распределения F(3)(x). Также на рис. 8 представлены графики функции F3(030200)(x) при значениях: 3 — p2 = 0,3; 4 — p2 = 0,5; 5 — p2 = 0,7 и ^ = 40.
ГП(Ы"*2>Ы 7
О
х 1
Рис. 8. Эмпирическая функция распределения
Разработанная программа моделирования (см. рис. 7) позволила провести численные эксперимента при различных значениях параметров р2 и начальных условиях аь а2. Полученные эмпирические функции распределения при больших значениях а2 аналогичны представленным на рис. 7 и позволяют сделать следующий вывод. Поскольку при стремлении параметра р2 к единице эмпирическая функция распределения сходится к ступенчатой с единственным разрывом в нуле, то в случае р2 = 1 имеет место нормировка для числа финальных частиц, отличная от линейной. Таким образом, в рассматриваемой стохастической эпидемической модели случай удаления р > 0 и случай неудаления р = 0 из популяции зараженных особей имеют качественное различие. В первом случае число выживших здоровых особей примерно пропорционально начальному их числу, а во втором случае число выживших здоровых особей — меньше.
Решение задачи о нахождении нормировки для п(а1'а2) в случае р2 = 1 и вывод предельной теоремы для числа финальных частиц возможны аналитическими методами через построение решения уравнения (20) [18].
Брюсселятор. Стационарное распределение. Рассматривается марковский процесс = («^(¿),£2(¿)), t Е [0, то), на множестве состояний N2, соответствующий схеме взаимодействий [6]
2T1 + T2 ^ 371; T1 ^ 0, T2; 0 ^ Ть
(23)
Второе уравнение Колмогорова для производящей функции переходных вероятностей Т(а11а2)(£; 81, 82) имеет вид
д^а(*; в) = . ( 3 , А д3ЗД 8) , = Л2 ( 81 — 8182)
д£ ^ 1 1 4 д281д82 + Л^Р282 + Ро - 81) д 8' + Ло (81 - Ч в), Та(0; в) =
1 (24)
где Л2 > 0, Л1 > 0, Л0 > 0, р0 > 0, р2 > 0, р0 + р2 = 1. Рассматриваемый марковский процесс является двумерным процессом рождения и гибели кубического типа; возможные скачки процесса изображены на рис. 6, случай б.
Детерминированная модель брюсселятора выводится аналогично, как в предшествующем пункте; из уравнения (24) приходим к дифференциальным уравнениям
Ж1 = Л2^2Ж2 - Л1Х1 + Ло, Ж2 = -Л2^2Ж2 + Л1Р2Х1,
Ж1(0) = Ж?, Ж2(0) = ж0, ()
где ж1 (£), ж2(£) — количество реагентов в момент времени £ для реакции с кинетической схемой (23). Исследованиям поведения решений системы уравнений (25) посвящена обширная литература (см., например, [6, 7]). Для приложений важны случаи поведения траекторий, представленные на рис. 9 (устойчивый фокус) и рис. 11 (предельный цикл).
На рис. 9 дана траектория детерминированной модели ж1 (£), ж2(£) на фазовой плоскости ж1 Ож2 при начальных условиях ж1(0) = 1000, ж2(0) = 3000 и приведен пример реализации стохастического процесса (6(£),6(£)) при начальных условиях (6(0),6(0)) = (1000,3000). Значения параметров составляют Л2 = 0,000001, Л1 = 7, Л0 = 4000, Р2 = 0, 7, Т = 10.
На рис. 11 дана траектория детерминированной модели ж1(£), ж2(£) на фазовой плоскости ж1 Ож2 при начальных условиях ж1(0) = 1000, ж2(0) = 3000 и приведен пример реализации стохастического процесса (6(£),&(£)) при начальных условиях (6(0),6(0)) = (1000,3000). Значения параметров составляют Л2 = 0,000001, Л1 = 8, Л0 = 4000, Р2 = 0,7, Т = 10.
Схема (23) интерпретируется как открытая система частиц с тройными взаимодействиями. При £ ^ то существует стационарное распределение, характеризуемое распределением вероятностей
{?(в1.02^ (в1,в2) Е N2}, где 9(в1 ,в2) = Ит^ Р((з111в22)(£).
1000 1200 то 1600 1800 2000 2100 2400 2600
Рис. 9. Детерминированная и стохастическая реализации в случае устойчивого фокуса
Рис. 10. Двумерная гистограмма стационарного распределения при Ах = 7
Уравнение для производящей функции стационарных вероятностей
то
I(-1 ,-2)= Е 9(в1 )£# , |-1|< 1, |-2 1< 1, 01,02 =0
имеет вид
Л2 (-1 - ^2 ) ^1 (5)
d2s1 ds2
+
f (s)
+ S2 + P0 _ S^ + A„(S1 — 1) f (s) = 0. (26)
4500
то
3500 3000 2500 2000
1500500 1000 1500 2000 ¿500 3000
Рис. 11. Детерминированная и стохастическая реализации в случае предельного цикла
Л
'/|,>У/>'.!>х10-6
Рис. 12. Двумерная гистограмма стационарного распределения при Ах =8
Исследование свойств стационарного распределения путем решения уравнения третьего порядка (26) аналитическими методами затруднительно. Также не представляет интереса решение линейного уравнения в частных производных (26) сеточными методами, так как приближенное решение не дает возможности определить коэффициенты разложения функции /($1, з2) в степенной ряд в нуле, т.е. искомые стационарные вероятности.
Методом статистического моделирования для схемы (23) получены
оценки вероятностей стационарного распределения д(в11(з2), (във2) £ € N2, вычисленные как отношение суммарного времени нахождения процесса в состоянии (във2) ко всему времени моделирования Т:
т
?(вьв2) = Е Т(в1,в2)/Т, т-(в1в2) — время нахождения процесса в состоянии (в1,в2) при /-м попадании в это состояние, ?(вьв2) = 1- На
^ в1,в2=0 рис. 10 приведена гистограмма для ?(в1,в2) при значениях параметров
Л2 = 0,000001, Л1 = 7, Л0 = 4000; р2 ' = 0,7; Т = 50000. На рис. 12 дана гистограмма для ?(в1,в2) при значениях параметров Л2 = 0,000001, Л1 = 8, Л0 = 4000; р2 = 0,7; Т = 50000. Возможность получения "кратерообразной поверхности" для стационарных вероятностей состояний обсуждается в работе [7] в связи с приложениями схемы (23) в теории диссипативных структур.
Статистическое моделирование марковского процесса £ (£) со схемой взаимодействий (23) при различных значениях параметров приводит к выводу, что стохастические реализации близки к детерминированным траекториям лишь при значениях параметров из очень небольшой области по отношению ко всему пространству параметров. В большинстве случаев стохастические реализации брюсселятора носят вырожденный характер, так как большую часть времени случайный процесс проводит в точках около границы множества N2. Варианты реализаций, близкие к представленным на рис. 9 и рис. 11, возможны при значениях параметров Л», значительно отличающихся по порядку. В некоторых случаях гистограмма стационарного распределения ?(в1,в2), (във2) € N2, оказывается близкой к плотности распределения двумерного нормального закона, что проверялось с помощью статистических критериев согласия. Детальному исследованию поведения стохастической модели брюсселятора и стационарного распределения в зависимости от параметров будет посвящена отдельная публикация.
СПИСОК ЛИТЕРАТУРЫ
1. Леонтович М. А. Основные уравнения кинетической теории газов с точки зрения теории случайных процессов // Журн. эксперим. и теорет. физики. -1935.- Т. 5, вып. 3-4.- С. 211-231.
2. Гихман И. И., Скороход А. В. Введение в теорию случайных процессов. -М.: Наука, 1977.-568 с.
3. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. Изд. 2-е.-М.: Наука, 1982.
4. Г а р д и н е р К. В. Стохастические методы в естественных науках. - М.: Наука, 1986.-528 с.
5. ВанКампен Н. Г. Стохастические процессы в физике и химии. - М.: Высшая школа, 1990. - 376 с.
6. Гленсдорф П., Пригожин И. Термодинамическая теория структуры, устойчивости и флуктуаций. Изд. 2-е. - М.: УРСС, 2003. - 280 с.
7. Малек-Мансур М., Николис Г., Пригожин И. Неравновесные фазовые переходы в химических системах // Термодинамика и кинетика биологических процессов. - М.: Наука, 1980. - C. 59-83.
8. Эмануэль Н. М., Кнорре Д. Г. Курс химической кинетики. - М.: Высшая школа, 1974. - 400 с.
9. McQuarrie D. A., Jachimowcki C. J., Russel M. E. Kinetic of small system. II // J. Chim. Phys. - 1964. - V. 40, № 10. - P. 2914-2921.
10. G i 11 e s p i e D. T. Approximate accelerated stochastic simulation of chemically reacting systems//J. Chim. Phys. - 2001. - V. 115, № 4.-P. 1716-1733.
11. Базыкин А. Д. Математическая биофизика взаимодействующих популяций. -М.: Наука, 1985. - 182 с.
12. Б е й л и Н. Математика в биологии и медицине. - М.: Мир, 1970. - 328 с.
13. W e i s s G. On the spread of epidemics by carries // Biometrics. - 1965. - V. 21, № 2.-P. 481-490.
14. Э п и д е м и и процесс // Математическая энциклопедия. Т. 5. - М.: Советская энциклопедия, 1985. - Кол. 1008.
15. Пичугин Б. Ю., Перцев Н. В. Статистическое моделирование популяций взаимодействующих частиц с произвольным распределением времени жизни // Матем. структуры и моделирование. - 2001. - Вып. 7. - С. 67-78.
16. Калинкин А. В. Финальные вероятности ветвящегося процесса с взаимодействием частиц и процесс эпидемии // Теория вероятн. и ее примен. - 1998. -Т. 43, вып. 4. - С. 773-780.
17. Калинкин А. В. Марковские ветвящиеся процессы с взаимодействием // Усп. матем. наук. - 2002. - Т. 57, вып. 2. - С. 23-84.
18. Мастихин А. В. Функция Римана для стационарного уравнения марковской эпидемии // Обозрение прикл. и промышл. математики. - 2003. - Т. 10, вып. 2. -C. 502.
19. Ланге А. М. Стационарное распределение в открытой стохастической системе с парными взаимодействиями частиц // Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". - 2005. - № 1. - C. 3-22.
20. Калинкин А. В., Шапошников А. А. Марковский процесс "хищник-жертва" и результаты экспериментов Г.Ф.Гаузе // Обозр. прикл. и промышл. матем. - 2005. - Т. 12 (В печати.).
Статья поступила в редакцию 10.03.2005
Александр Вячеславович Калинкин родился в 1956 г., окончил в 1978 г. МГУ им. М.В. Ломоносова. Д-р физ.-мат. наук, профессор кафедры "Высшая математика" МГТУ им. Н.Э. Баумана. Автор более 40 научных работ в области теории вероятностей и математического моделирования.
A.V. Kalinkin (b. 1956) graduated from the Moscow State University n.a. M.V. Lomonosov in 1978. D. Sc. (Phys.-Math.), professor of "Higher Mathematics" department of the Bauman Moscow State Technical University. Author of more than 40 publications in the field of theory of probabilities and mathematical simulation.
Андрей Михайлович Ланге родился в 1979 г., окончил в 2002 г. МГТУ им. Н.Э. Баумана. Ассистент кафедры "Высшая математика" МГТУ им. Н.Э. Баумана. Автор 3 научных работ в области математического моделирования.
A.M. Lange (b. 1979) graduated from the Bauman Moscow State Technical University in 2002. Assistant of "Higher Mathematics" department of the Bauman Moscow State Technical University. Author of 3 publications in the field of mathematical simulation.
Антон Вячеславович Мастихин родился в 1962 г., окончил в 1983 г. МГУ им. М.В. Ломоносова. Старший преподаватель кафедры "Высшая математика" МГТУ им. Н.Э. Баумана. Автор 4 научных работ.
A.V Mastikhin (b. 1962) graduated from the Moscow State University n. a. M.V Lomonosov in 1983. Senior teacher of "Higher Mathematics" department of the Bauman Moscow State Technical University. Author of 4 publications.
Андрей Анатольевич Шапошников родился в 1980 г., окончил в 2004 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Высшая математика" МГТУ им. Н.Э. Баумана. Специализируется в области статистического моделирования.
A.A. Shaposhnikov (b. 1980) graduated from the Bauman Moscow State Technical University in 2004. Post-draduate of "Higher Mathematics" department of the Bauman Moscow State Technical University. Specializes in the field of statistic simulation.
В издательстве МГТУ им. Н.Э. Баумана вышла книга
ТЕПЛОТЕХНИКА
Учебник для втузов / AM. Архаров, И.А. Архаров, В.Л. Бон-даренко, Б.П. Борисов и др.; Под общ. ред. А.М. Архарова, В.Н. Афанасьева. - 2-е изд., перераб. и доп. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. - 712 с.: ил.
Во втором, исправленном и дополненном издании (1-е изд. - под ред. В.И. Крутова 1986 г.) учебника, имеющего энциклопедический характер, рассмотрены основы термодинамики и теории теплообмена, топливо и его горение, схемы и элементы расчета котлов, промышленных печей, паро- и газотурбинных установок, холодильных установок и компрессоров, двигателей внутреннего и внешнего сгорания, ракетных и авиационных двигателей, атомных и плазменных энергоустановок. Приведены расчеты систем отопления, вентиляции и кондиционирования воздуха. Кроме того, включены важные разделы, касающиеся космических энергоустановок, теплообменных аппаратов, гидромашин, фотонных энергосистем, криогенных систем для ожиженных газов, разделения воздуха, получения неона, криптона и ксенона, термоэлектрических и термомагнитных низкотемпературных установок, а также систем регулирования. К работе над новыми разделами были привлечены известные специалисты. Большое внимание уделено вопросам экологии и защиты окружающей среды.
Содержание учебника соответствует курсам лекций, которые авторы читают в МГТУ им. Н.Э. Баумана и других крупных российских и зарубежных университетах.
Для студентов высших учебных заведений, обучающихся по направлению "Энергомашиностроение"
По вопросам приобретения обращаться по тел. 263-60-45;
e-mail: press@bmstu.ru
