CC BY
67
УДК 519.21;51-76
Влияние стохастизации на одношаговые модели
А. В. Демидова*, М. Н. Геворкян*, А. Д. Егоров^, Д. С. Кулябов*, А. В. Королькова*, Л. А. Севастьянов*
* Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198
^ Институт математики НАНБ ул. Сурганова, д. 11, г. Минск, Белоруссия, 220072
Предполагается, что введение стохастики в математическую модель делает её более адекватной. При этом практически отсутствуют методы согласованного (зависящего от структуры системы) введения стохастики в детерминистические модели. Авторами была усовершенствована методика построения стохастических моделей для класса одноша-говых процессов и проиллюстрирована на примере моделей популяционной динамики. Популяционная динамика была выбрана для исследования потому, что её детерминистические модели достаточно хорошо исследованы, что позволяет сравнить полученные результаты с уже известными.
В работе изучено влияние введения стохастики в детерминистические модели на примере системы популяционной динамики типа «хищник—жертва». Полученные ранее стохастические дифференциальные уравнения исследуются методами качественной теории дифференциальных уравнений. Получено стационарное состояние и первый интеграл системы. Для демонстрации результатов производится численное моделирование на основе метода Рунге—Кутты для стохастических дифференциальных уравнений. Первый интеграл детерминистической системы (фазовый объём) в стохастическом случае не сохраняется, а возрастает, что в конечном итоге приводит к гибели одной или обеих популяций.
Одним из недостатков классической системы типа «хищник—жертва» считается сохранение амплитуды колебаний популяций. В стохастической же модели процесс завершается гибелью одной или обеих популяций, что, с точки зрения авторов, делает модель более адекватной.
Ключевые слова: стохастические дифференциальные уравнения, модель «хищник-жертва», основное кинетическое уравнения, уравнение Фоккера-Планка.
1. Введение
Данная работа находится в русле проводимых нами исследований по стохастизации математических моделей. Наш интерес к данной тематике проистекает из следующих проблем: построение популяционных моделей из первых принципов и введение стохастики в модели данного вида (следует заметить, что обратились мы к популяционной динамике потому, что нами исследовались схожие модели в других областях).
При стохастизации математических моделей возникает проблема введения стохастического члена. Существует несколько способов сделать это. Самый простой вариант — аддитивное добавление стохастического члена к детерминистическому уравнению. Однако при таком введении возникают свободные параметры, требующие дальнейшего определения из дополнительных соображений. Кроме того, данные стохастические члены обычно интерпретируются как внешнее (а не структурное) случайное воздействие. В связи с этим мы использовали и усовершенствовали метод построения стохастических моделей одношаговых процессов на основе основного кинетического уравнения [1,2]. Стохастическое дифференциальное уравнение рассматривается как его приближённая форма. Это позволяет получить модельные уравнения из общих принципов. Кроме того, детерминистическая и стохастическая части получаются из одного и того же уравнения, что
Статья поступила в редакцию 24 декабря 2013 г.
Работа выполнена при частичной финансовой поддержке гранта БРФФИ-ОИЯИ №198 (от 06.002.2012 г.)
рассматривается нами как согласованность стохастической и детерминистической частей. Полученный стохастический член интерпретируется нами как структурная стохастика.
Целью данной работы является выработка методики исследования влияния введения стохастики на исходную детерминистическую модель. Для верификации результатов предлагается использовать различные численные методы построения решений стохастических дифференциальных уравнений, а также сравнение с известными результатами, которые получены и представлены в литературе [3-5]. Для численного решения уравнений детерминированной и стохастической моделей используются методы Рунге-Кутты разных порядков [6,7].
Структура статьи следующая. В разделе 2 введены основные обозначения и соглашения. В разделе 3 даётся краткое введение в метод стохастизации одноша-говых процессов. Далее, в разделе 4 описывается исследуемая модель. При этом в подразделе 4.1 даётся краткая справка по стандартному (детерминистическому) подходу, а в подразделе 4.2 мы получаем стохастическое расширение данной модели по методу стохастизации одношаговых процессов. В разделе 5 рассматривается методика исследования влияния стохастического члена на поведение решений модели на примере системы типа «хищник-жертва». В разделе 6 рассматривается возможность применения методов Рунге-Кутты для анализа стохастических дифференциальных уравнений и проводится численное исследование полученных ранее уравнений, служащее иллюстрацией к аналитическим выкладкам.
2. Обозначения и соглашения
1. В работе используется нотация абстрактных индексов. В данной нотации тензор как целостный объект обозначается просто индексом (например, хг), компоненты обозначаются подчёркнутым индексом (например, х-).
2. Будем придерживаться следующих соглашений. Латинские индексы из середины алфавита (г, ], к) будут относиться к пространству векторов состояний системы. Латинские индексы из начала алфавита (а) будут относиться к пространству винеровского процесса. Латинские индексы из конца алфавита (р, О) будут относиться к индексам метода Рунге-Кутты. Греческие индексы (а) будут задавать количество разных взаимодействий в кинетических уравнениях.
3. Точкой над символом обозначается дифференцирование по времени.
4. Запятой в индексе обозначается частная производная по соответствующей координате.
3. Моделирование одношаговых процессов
Под одношаговыми процессами мы будем понимать марковские процессы с непрерывным временем, принимающие значения в области целых чисел, матрица перехода которых допускает только переходы между соседними состояниями. Также эти процессы известны как процессы рождения-гибели. Одношаговые процессы подчиняются следующим условиям:
1) если в момент времени £ система находится в состоянии % £ то вероятность перехода в состояние % + 1 в интервале времени [1, £ + ДЦ равна + о(Д£);
2) если в момент времени £ система находится в состоянии % £ Ъ+, то вероятность перехода в состояние г — 1 в интервале времени [1,1 + Д£] равна к- АЪ + о(АЪ);
3) вероятность перехода в состояние, отличное от соседних, равна о(АЪ);
4) вероятность сохранения прежнего состояния равна 1 — (к+ + к-)ДЬ + о(ДЬ);
5) состояние % = 0 есть поглощающая граница.
Идея метода введения стохастики в модель состоит в следующем. На основании схем взаимодействия мы строим основное кинетическое уравнение, раскладываем его в ряд, оставляя только члены до второй производной включительно.
Получившееся уравнение будет уравнением Фоккера-Планка. Для получения более привычного вида модели записываем соответствующее ему уравнение Ланже-вена. На самом деле, как мы увидим, из схем взаимодействия мы сразу получаем коэффициенты уравнения Фоккера-Планка (и, соответственно, уравнения Лан-жевена), поэтому при практическом применении метода строить основное кинетическое уравнение нет необходимости.
3.1. Схемы взаимодействия
Состояние системы будем описывать вектором состояния хг £ К п, где п — размерность системы (под вектором состояния будем понимать множество математических величин, полностью описывающих систему). Оператор п^ £ задаёт состояние системы до взаимодействия, оператор т^ £ — после. В
результате взаимодействия происходит переход системы в другое состояние [1,8].
В системе может происходить 8 видов различных взаимодействий, где 8 £ Поэтому вместо операторов п* и т^ будем рассматривать операторы пг0а £ Щ>0 х х и т,га £ х х ^0. Взаимодействие элементов системы будем описывать с помощью схем взаимодействия, подобным схемам химической кинетики [9]:
к+
пгххг ^ тгххг, (1)
г к-
ка
здесь греческие индексы задают количество взаимодействий, а латинские — размерность системы.
Изменение состояния будет задаваться оператором
гТ=тТ — *гт- (2)
Таким образом, один шаг взаимодействия а в прямом и обратном направлениях соответственно можно записать в виде:
хг ^ хг + г^хг, хг ^ хг — грхг. (3)
Мы также можем записать (1) не в виде векторных уравнений, а в виде более традиционных сумм:
к+
П^х^5г ^ т^х^ 5г, (4)
ка
где 6г = (1,..., 1).
Будем использовать также следующие обозначения:
пга := п^Р, тга := т^Р, гга := г)аЬг. (5)
3.2. Основное кинетическое уравнение
Вероятности перехода в единицу времени из состояния хг в состояние хг+г~хг
(в состояние хг — г~х^) пропорциональны соответственно числу способов выбора комбинации из хг по пга (из хг по тга) и определяются выражениями:
"" г! п Тг-!
*+П (Хг — п»)!, = к- П (Ж1 — тг-)!. (6)
г=1 у ' г=1 у '
Таким образом, общий вид основного кинетического уравнения для вектора
- г 1а 7
состояний ж', изменяющегося шагами длины г~ж^, принимает вид:
др(хг ,t) dt
rn
L
a=1
s-(xi + r™, t)p(xi + ri-, t) - s++(xi)p(xi,t)
+
+
8+(х* - ri^1t)p(xi - r%t) - s'-(xi) p(xi,t) . (7)
3.3. Уравнение Фоккера-Планка
Далее, используя разложение Крамерса-Мойала, получаем уравнение Фоккера-Планка [9]. Для этого делается несколько предположений:
1) имеют место только малые скачки, т.е. 8а(хг) является функцией медленно изменяющейся с изменением хг;
2) р(хг,1) также медленно изменяется с изменением хг.
Можно выполнить сдвиг в (7) из точки (хг ± г~х^) в точку хг и, разложив правую часть в ряд Тейлора и отбросив члены порядка выше второго, получим уравнение Фоккера-Планка:
^ = -а, [А'р] + 2didj [В*р]
где
Ai :=A\xk,t) = г^ Bij := Bij (хк, t) = r^r^
S a S о
S a S о
(8)
(9)
a = l,m.
Как видно из (9), коэффициенты уравнения Фоккера-Планка можно получить сразу из (2) и (6), т. е. в практических расчётах записывать основное кинетическое уравнение нет необходимости.
3.4. Уравнение Ланжевена
Уравнению Фоккера-Планка соответствует уравнение Ланжевена:
d^ = aidt + b{dWа, (10)
где аг := аг(хк, t), Ьга := Ьга(хк, t), хг £ Rn — вектор состояния системы, Wа £ rto — m-мерный винеровский процесс. Винеровский процесс реализуется как d W = eVdt, где £ ~ N(0,1) — нормальное распределение со средним 0 и дисперсией 1. Латинскими индексами из середины алфавита обозначаются величины, относящиеся к векторам состояний (размерность пространства п), а латинскими индексами из начала алфавита обозначаются величины, относящиеся к вектору винеровского процесса (размерность пространства m < п).
При этом связь между уравнениями (8) и (10) выражается следующими соотношениями:
Ai = ai, Bij = blbja. (11)
При интегрировании уравнения Ланжевена с непостоянными коэффициентами у винеровского процесса (Ьга = const) возникает неопределённость учёта скачка. Данная проблема снимается с введением определённой интерпретации стохастического уравнения (интерпретация Ито или интерпретация Стратановича). Мы будем использовать интерпретацию Ито.
В рамках интерпретации Ито дифференциал от сложной функции не подчиняется стандартным формулам анализа. Для его вычисления используется правило или лемма Ито:
Лемма (см. [10]). Пусть / := /(хк,1) — функция от случайного процесса хк(£), / £ С2. Тогда формула дифференциала будет выглядеть следующим образом:
= к/ + аг!г + 1 Ь\Ьг 1 М + ЪЦг (12)
где } := }(хк,Ь), аг := аг(хк,Ь), Ьа := Ьаа(хк,Ь) и dWa := <Ша(г).
4. Модель «хищник—жертва» 4.1. Детерминистическая модель «хищник—жертва»
Системы с взаимодействием двух видов популяций типа «хищник-жертва» широко исследованы, и для таких систем существует большое количество разнообразных моделей. Самой первой моделью «хищник-жертва» принято считать модель, полученную независимо друг от друга А. Лоткой и В. Вольтеррой. Лотка в своей работе [11] описывал некоторую гипотетическую химическую реакцию:
А к X —^ У —Ч В, (13)
где Х,У — промежуточные вещества, коэффициенты ^1,^2,^3 — скорости химических реакций, А — исходный реагент, а В — продукт реакции. В результате была получена система дифференциальных уравнений вида:
(х = кцх — к2ху, (14)
= к2ху — кзу.
Эта система полностью совпадает с системой дифференциальных уравнений, полученной Вольтеррой, который рассматривал механизм роста численности двух популяций с взаимодействием типа «хищник-жертва». Для получения уравнений Вольтерра в своей работе [3] делает следующие идеализированные предположения о характере внутривидовых и межвидовых отношений в системе «хищник-жертва»:
— для того чтобы охарактеризовать одним числом некоторую популяцию в ограниченной области, делается допущение об однородности (не учитывается возраст и размер), также предполагается, что тип индивидуума не меняется со временем;
— численность каждого вида увеличивается и уменьшается непрерывным образом, а функция, описывающая численность вида, является непрерывно дифференцируемой функцией;
— в отсутствие хищника популяция жертвы размножается в соответствии с принципом Мальтуса — экспоненциально, а популяция хищника в отсутствие жертвы экспоненциально вымирает;
— суммарное количество жертвы, потребляемое популяцией хищника в единицу времени, линейно зависит и от численности популяции жертвы, и от численности популяции хищника, а потреблённая хищником биомасса жертвы с постоянным коэффициентом перерабатывается в биомассу хищника;
— дополнительные факторы, оказывающие влияние на динамику популяций (такие как внутривидовая конкуренция, ограниченность ресурсов и т.д.), не учитываются.
Рассматривается случай, когда в ограниченной среде сосуществуют два вида, один из которых (хищник) питается за счёт второго вида (жертва). Очевидно,
что численность жертв будет увеличиваться тем медленнее, чем больше существует хищников, а хищников тем быстрее, чем многочисленнее жертвы. Воль-терра получил следующую систему дифференциальных уравнений для описания численности видов:
1^1 = (£ 1
= (-е 2 + 72^1)^2. ( )
где N1 — численность жертв, N2 — численность хищников, £1 и £2 — положительные постоянные коэффициенты, отражающие естественную рождаемость и смертность жертв и хищников соответственно, а 71 и 72 — положительные постоянные коэффициенты, описывающие межвидовое взаимодействие. Данная система эквивалентна системе (14) с точностью до обозначений.
4.2. Стохастическая модель «хищник—жертва»
Рассмотрим модель системы «хищник-жертва», состоящую из особей двух видов, причём один из них охотится, второй — обеспечен неисчерпаемыми пищевыми ресурсами. Введя обозначения X — жертва, У — хищник, можно записать возможные процессы (4) для вектора состояния х- = (Х,У)Т [12-15]:
X 2Х, г1-1 = (1, 0)Т, X + У 2У, г*2 = (-1,1)Т, (16)
у 0, г13 = (0,-1)Т,
которые имеют следующую интерпретацию. Первое соотношение означает, что жертва, которая съедает единицу пищи, немедленно репродуцируется. Второе соотношение описывает поглощение жертвы хищником и мгновенное репродуцирование хищника. Это единственная рассматриваемая возможность гибели жертвы. Последнее соотношение — естественная смерть хищника.
Все процессы необратимы, поэтому в - = 0, а
+1 х! У!
81 (х, У) = к1 (х _ 1)! = К]_х, X! у!
4(х, У) = к2 ^ - ^ ^ - ^ =к2ху, (17)
+ , х! у!
8 +(Х, У) = кз=кзУ
Воспользовавшись формулой (8), получим уравнение Фоккера-Планка:
= -дг {А\х, у)р(х, у)) + 1дгдэ (ВV (Х, у)р(х, у)) , (18)
где
Аг(х, у) = з+(х, у) гга, В" (X, у) = з+(х, Ууаг*а. (19)
Таким образом мы получили:
»)=(о) м+(-1) к2*»+кзу=(к:;--^), (20)
в^(X, у) = (1,0)к1Х + ( (—1,1)к2ХУ + (—^ (0, —1)кзу
(к1Х + к2ху —к2ху \ ^ — к2ху к2ху + кзу). ( )
Для того чтобы записать стохастическое дифференциальное уравнение в форме Ланжевена (10) для модели «хищник-жертва», достаточно извлечь квадратный корень из полученной матрицы Вгз в уравнении Фоккера-Планка:
(х\ = (к1Х — к2ху\ л ^ й\у) = \k2xy — кзу) ё + Н2 ,
ФЫ-а = Вг£ = (к1Х + к2ХУ —к2Ху 4 а ^ —к2ху к2ху + кзУ;
Следует заметить, что конкретный вид матрицы га не выписан из-за крайней громоздкости выражения. Впрочем, при дальнейших исследованиях нам понадобится не собственно матрица Ьга, а её квадрат, то есть матрица Вг3.
5. Исследование влияния стохастического члена 5.1. Детерминистическая модель
Качественное исследование детерминистической модели «хищник-жертва» (14) описано в многочисленной литературе [16-18] и др. Приведём основные результаты.
5.1.1. Репараметризация уравнений
Для уменьшения количества свободных параметров сделаем репараметриза-цию изучаемых уравнений.
Пусть х(1) = Ри(т), у(£) = Qv (т), £ = Тт, Р, Q иТ — калибровочные константы. Зададим вид репараметризованной системы:
{и = и(1 — V), V = ау(и — 1),
где точкой обозначена производная по т (как и везде в этом пункте). Получим:
{и = и (к1Т — к^Ть), у = у (к2ТРи — кзТ).
Сравнивая (23) и (24), получаем явный вид для замены переменных:
(23)
(24)
и(т) = ^х(г), ь(т) = ^ у(г), Т = к1 г, а = . (25)
Кз К1 К1
5.1.2. Стационарные точки
Найдём стационарные состояния системы (23), которые являются решением системы уравнений:
!и = 0 = и(1 — V), , ч
• 0 \ 1 (26) V = 0 = ау (и — 1),
Система (23) имеет два стационарных состояния: (0, 0) и (1,1). Точка (0,0) определяет положение равновесия, которое характеризуется полным истреблением жертв и вымиранием хищников. Точка (1, 1) отражает стационарный режим сосуществования хищников и жертв с некоторыми ненулевыми численностями.
5.1.3. Исследование линеаризованной устойчивости
Линеаризуем систему (23). Пусть и = и + £, V = и + г], где и и V — координаты точки равновесия, а и — малые возмущения:
= - V) - 'ПЩ (27)
\ 'г] = щ (и - 1) +
Запишем линеаризованную систему в окрестности каждой из точек равновесия. Рассмотрим точку (0, 0), тогда система (27) примет вид:
И = ^ (28)
\г] = - ощ.
Найдём собственные значения характеристического уравнения
(1 -Х)(-а -Х) = 0, (29)
которое соответствует системе (28). Таким образом, имеем А1 = 1, А2 = -о, собственные значения действительные и имеют разные знаки, значит, точка (0, 0) — седло.
Для точки (1, 1) линеаризованная система будет иметь вид:
{« = 7 (30)
[V =
Характеристическое уравнение
А2 + а = 0 (31)
имеет корни А1 = \у/а и А2 = - 1л/а. Так как собственные значения — мнимые, то точка (1, 1) является центром, т.е. по крайней мере в окрестности этой точки существуют замкнутые траектории.
Системы (14) и (23) — консервативные. Чтобы показать это, составим первый интеграл системы и построим её фазовый портрет.
Фазовые кривые системы (23) являются интегральными кривыми уравнения
ёи = и(1 - у) ё V ау(и - 1)
После решения получим первый интеграл системы:
1(и, у) = ои + у - 1п уиа. (33)
Действительно, можно проверить, что производная Ли от него равна нулю:
£ / I = 1,г! г = 0, (34)
где векторное поле /г есть правая часть системы (23): иг = /г(ик), иХ = (и, у)Т.
Фазовый портрет системы в окрестности стационарной точки (1,1) представляет собой замкнутые эллиптические орбиты (рис. 1).
Детерминированная модель Лотки-Вольтерры (фазовый портрет)
хШ — численность жертв
Рис. 1. Фазовый портрет для системы «хищник—жертва» в окрестности
точки (1, 1)
Таким образом изменение численности обоих видов происходит по периодическому закону с амплитудой колебаний, определяемой начальными значениями и и . Решения имеют осциллирующие зависимости, показанные на рис. 2.
Детерминированная модель Лотки-Вольтерры (графики решений)
I — время
Рис. 2. Изменение численности хищников N1 и жертв N2 во времени
Циклы на рис. 2 повторяются неограниченно долго и качественно отражают свойства многих реальных систем «хищник-жертва» [16,19].
Основная особенность модели Вольтерра, благодаря которой она стала классической для многих последующих моделей математической экологии, состоит в том, что на основе очень упрощённых представлений о характере закономерностей, описывающих поведение системы сугубо математическими средствами, было выведено заключение о качественном характере поведения такой системы и о наличии в системе колебаний численности популяций.
В то же время этой системе присущи два принципиальных и взаимосвязанных недостатка — один математического, другой биологического характера. С математической точки зрения система (14) — негрубая и консервативная. Это означает, что включение в модель каких бы то ни было дополнительных факторов качественным образом меняет её поведение. С биологической точки зрения недостатком модели является то, что в рамки её не включены принципиальные свойства, характерные для любой пары взаимодействующих по принципу
«хищник-жертва» популяций: эффект насыщения хищника, ограниченность ресурсов жертвы и хищника даже при избытке жертвы и т. п.
5.2. Стохастическая модель
Рассмотрим изменение качественного поведения системы (14) при введении стохастического члена (система (22)).
Используя (33), запишем первый интеграл для детерминистической части системы (22):
1(х, у) = к2(х + у) - к3 1пх - к11п у. (35)
Также можно проверить, что производная Ли от него равна нулю:
С ¡1 = г = 0, (36)
где векторное поле /г есть правая часть системы (14): хг = /г(хк), х- = (х, у)Т. Далее воспользуемся формулой Ито (12) для функции ё 1(х, у):
ё 1(х, у) = 1,гаг<И + 2Ь1Ьга1гё + ЬЦ,гё^а
(37)
Учитывая, что первый член в (37) равен нулю согласно (36), а также то, что среднее от винеровского процесса равно нулю ((ё ) = 0), запишем формулу для среднего изменения фазового объёма:
1
»» = 2 (В' 2(1)>
11 к3
+ В
22
к
2( )
) = 1 (
1 /кхкз + к2кзу + к^х + (38)
Поскольку х,у £ К^о, то видно, что в стохастической модели фазовый объём в среднем монотонно возрастает. В конце концов задевается одна из осей, что приводит к гибели одной или обеих популяций.
Данное поведение проиллюстрировано на рис. 3. При этом временная зависимость имеет вид, приведённый на рис. 4.
4.8 4.6 4.4 4.2 4 3.8 3.6 3.4
3.2
10
Стохастическая модель Лотки-Вольтерры
Хищник-Жертва
20 25
хШ — численность жертв
15
30
35
Рис. 3. Фазовый портрет для стохастической модели «хищник—жертва» в
окрестности точки (1,1)
Стохастическая модель Лотки-Вольтерры
I — время
Рис. 4. Стохастическая модель «хищник—жертва»
6. Численное моделирование
В целях визуализации полученных результатов было проведено численное моделирование. Даётся краткое описание подходов к численной реализации стохастических дифференциальных уравнений и даётся обоснование конкретного метода для модели типа «хищник-жертва».
6.1. Стохастические методы Рунге-Кутты с сильным порядком точности, равным единице
Для численного решения стохастических дифференциальных уравнений (СДУ) существует множество подходов. Рассмотрим один из них, заключающийся в распространении методов Рунге-Кутты на случай СДУ. При этом следует различать два вида решений: сильное и слабое [20]. Сильное решение означает, что решение СДУ полностью определяется в момент времени заданной траекторией вине-ровского процесса Ша на отрезке [0, Ц и начальным условием х0. Слабое решение означает качественное построение вероятностной модели, где определены некоторый винеровский процесс Ша и процесс хг, удовлетворяющий стохастическому уравнению.
Иначе говоря, в случае сильного решения интересуются собственно траекториями пары процессов, а в случае слабого решения — распределением этой пары.
Ещё одной важной особенностью, о которой следует помнить при работе с численными методами для СДУ, является то, что, используя лишь приращение А Шк винеровского процесса, нельзя получить схемы сильного порядка точности выше первого [7,20].
Перейдём теперь к изложению основных сведений о стохастических методах Рунге-Кутты для СДУ вида:
Ах1 = а\хк + Ьа(хк )АШа.
Классический метод Рунге-Кутты можно распространить на случай СДУ следующим образом [21]:
Г 4 = х0 + (сф + ЗахЩ (др, I х\ = х0 + ^^(ф +
= 1, е.
где Jl ~ N(0, К) или ~ у/Ке, е ~ N(0,1), р,д — индексы метода Рунге-Кутты. Коэффициенты метода можно сгруппировать в таблицу, называемую таблицей Батчера:
Rl R4p
гя f-q
Конкретный метод полностью определяется своей таблицей Батчера [22-24]. Теория помеченных деревьев, используемая для определения уравнений порядка для коэффициентов классических методов Рунге-Кутты, была распространена на случай стохастических методов. Приведём здесь два конкретных стохастических метода Рунге-Кутты [21]:
0 0 0 2/3 0 0 -1 1 0 0 0 0 2/3 0 0 , IM: -1 1 0 0 0 0 1/4 1/4 0 0 1 0 0 0 0 1/4 1/4 0 0 1 0
0 3/4 1/4 0 3/4 1/4 1/6 2/3 1/6 1/6 2/3 1/6
Согласно численному эксперименту, проведённому в статье [21], наименьшую глобальную погрешность вычисления обеспечивает метод Ш, однако это неявный метод, что затрудняет реализацию его на компьютере. Из явных методов наилучшую точность даёт метод БИ2.
6.2. Численное решение уравнений стохастической модели
хищник—жертва
Следует отметить, что для детерминированной модели «хищник-жертва» метод Эйлера даёт слишком большую погрешность даже вблизи стационарных точек. Из этого следует неприменимость метода Эйлера-Маруямы к стохастической модели «хищник-жертва», так как он является распространением метода Эйлера на случай СДУ.
Для численного решения СДУ модели «хищник-жертва» были использованы вышеизложенные методы Рунге-Кутты, которые были реализованы на языке Фортран.
В дальнейшем предполагается реализовать и использовать стохастические методы Рунге-Кутты сильного порядка точности больше 1. В таких методах используются винеровские приращения большей кратности. Однако количество коэффициентов резко возрастает, а следовательно, резко возрастает сложность их вычисления [6,7,25].
6.3. Генерация псевдослучайных, нормально распределённых чисел
Для генерации нормально распределённых псевдослучайных чисел использовалась функция random_normal(), написанная на языке Фортран. Функция входит в состав модуля random, доступного по адресу http://www.netlib.org/ random/index.html (автор Alan Miller). Алгоритм, используемый в этой функции, описан в статье [26].
7. Заключение
Для описания эволюции популяционных систем в большинстве случаев используются детерминистские модели. Но они имеют свои недостатки, например,
не учитывают вероятностный характер процессов рождения и гибели индивидуумов популяций. От некоторых из этих недостатков позволяют избавиться стохастические модели, поэтому именно им в работе уделено основное внимание. Наиболее удобным для исследования является стохастическое дифференциальное уравнение в форме Ланжевена, так как в нём стохастическая и детерминистическая части разделены.
В работе описан метод комбинаторной кинетики и получено уравнение Фоккера-Планка для модели «хищник-жертва» и описан механизм получения для этих уравнений СДУ в форме Ланжевена.
Для модели «хищник-жертва» в работе показано применение описанного подхода построения стохастической модели, который позволяет получить универсальные правила записи стохастических дифференциальных уравнений для систем, процессы в которых представимы как одношаговые. Кроме того, анализ полученных уравнений позволяет сделать вывод, что стохастические модели дают более реалистичное описание поведения системы. В частности, для системы «хищник-жертва» в детерминистическом случае решения уравнений имеют периодический вид и фазовый объём сохраняется, в то время как введение стохастики в модель даёт монотонное возрастание фазового объёма, что говорит о неизбежной гибели одной либо обеих популяций.
Литература
1. Кулябов Д. С., Демидова А. В. Введение согласованного стохастического члена в уравнение модели роста популяций // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2012. — № 3. — С. 69-78. [Demidova A. V., Kulyabov D. S. 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". — 2012. — No. 3. — P. 69-78. — (in russian). ]
2. Demidova A. V., Sevastianov L. A., Kulyabov D. S. Application of Stochastic Differencial Equations to Model Population Systems // Third International Conference on Mathematical Modelling of Social and Economical Dynamics MMSED-2010 / Russian State Social University. — 2010. — Pp. 92-94.
3. Вольтерра В. Математическая теория борьбы за существование. — М.: Наука, 1976. [Volterra V. Mathematical Theory of the Struggle for Existence. — Moscow: Nauka, 1976. — (in russian). ]
4. Базыкин А. Д. Нелиненая динамика взаимодействующих популяций. — Москва-Ижевск: Институт компьютерных исследований, 2003. [Bazykin A. D. Nonlinear Dynamics of Interacting Populations. — World Scientific, 1998. ]
5. Братусь А. С., Новожилов А. С., Платонов А. П. Динамические системы и модели биологии. — М.: Физматлит, 2010. [Bratus A. S., Novozhilov A. S. and Platonov A. P. Dynamical Systems and Biology Models. — Moscow: Fizmatlit, 2010. — (in russian). ]
6. 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.
7. 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.
8. 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.
9. Гардинер К. В. Стохастические методы в естественных науках. — Мир, 1986. [Gardiner C. W. Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences. — Springer Series in Synergetics, 1985. ]
10. 0ksendal B. K. Stochastic Differential Equations: An Introduction with Applications. — Berlin: Springer, 2003.
11. Lotka A. J. Elements of Physical Biology. — BiblioBazaar, 2011. — ISBN 9781178508116, 492 p. — http://books.google.ru/books?id=tFN9pwAACAAJ.
12. Демидова А. В. Уравнения динамики популяций в форме стохастических дифференциальных уравнений // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2013. — № 1. — С. 67-76. [Demidova A. V. The Equations of Population Dynamics in the Form of Stochastic Differential Equations // Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics". — 2013. — No. 1. — P. 67-76. — (in russian). ]
13. Демидова А. В. Метод стохастизации математических моделей на примере системы «хищник-жертва» // Научная сессия НИЯУ МИФИ-2013. — 2013. — С. 127. [Demidova A. V. The Method of Stochastization of Mathematical Models for the Example of the "Predator-Prey" // A Scientific Session NRNU MEPHI-2013. — 2013. — P. 127. — (in russian). ]
14. Королькова А. В., Кулябов Д. С. Методы стохастизации математических моделей на примере пиринговых сетей // Научная сессия НИЯУ МИФИ-2013. Аннотации докладов. В 3 томах. — Москва: МИФИ, 2013. — С. 131. [Korolkova A.V., Kulyabov D.S. Methods of Stochastization of Mathematical Models on the Example of Peer to Peer Networks // A Scientific Session NRNU MEPHI-2013. — 2013. — P. 131. — (in russian). ]
15. Демидова А. В., Кулябов Д. С., Севастьянов Л. А. Согласованный стохастический член в популяционных моделях // XI Белорусская математическая конференция. — Минск: Институт математики НАН Беларуси, 2012. — С. 39. [Demidova A. V., Kulyabov D. S., Sevastianov L. A. The Agreed Stochastic Term in Population Models // XI Belarusian Mathematical Conference. — Minsk,
2012. — P. 39. — (in russian). ]
16. Марри Д. Нелинейные дифференциальные уравнения в биологии. Лекции о моделях. — М.: Мир, 1983. [Murray J. D. Lectures on Nonlinear-Differential-Equation Models in Biology. — Oxford: Clarendon Press, 1977. ]
17. Свирежев Ю. М., Логофет Д. О. Устойчивость биологических сообществ. — М.: Мир, 1983. [Svirezhev Y. M., Logofet D. O. Stability of Biological Communities. — Moscow: Mir, 1983. — (in russian). ]
18. Эрроусмит Д., Плейс К. Обыкновенные дифференциальные уравнения. Качественная теория с приложениями. — М.: Мир, 1986. [Arrowsmith D., Place C. Ordinary Differential Equations: A Qualitative Approach with Applications. — Chapman and Hall Mathematics, 1982. ]
19. Одум Ю. Основы экологии. — Москва: Мир, 1975. [Odum E. P. Fundamentals of Ecology. — Philadelphia: W. B. Saunders, 1953. ]
20. Лукшин А. В., Смирнов С. Н. Численные методы решения стохастических дифференциальных уравнений // Математическое моделирование. — 1990. — Т. 2, № 11. — С. 108-121. [Lukshin A.V., Smirnov S.N. Numerical Methods of Solving Stochastic Differential Equations // Mathematic Modeling. — 1990. — Vol. 2, No. 11. — P. 108-121. — (in russian). ]
21. Soheili A. R., Namjoo M. Strong Approximation of Stochastic Differential Equations with Runge-Kutta Methods // World Journal of Modelling and Simulation. — 2008. — Vol. 4, No 2. — Pp. 83-93.
22. Геворкян М. Н. Конкретные реализации симплектических численных методов // Вестник РУДН. Серия «Математика. Информатика. Физика». —
2013. — № 1. — С. 77-89. [Gevorkyan M. N. Specific Implementations of Symplectic Numerical Methods // Bulletin of Peoples' Friendship University of
Russia. Series "Mathematics. Information Sciences. Physics". — 2013. — No. 1. — P. 77-89. — (in russian). ]
23. Геворкян М. Н. Условие явности и диагональной неявности при композиции метода Рунге-Кутты со своим присоединенным // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2012. — № 3. — С. 87-96. [Gevorkyan M. N. The Condition of Explicitness and the Diagonal Implicitness for Composition of Runge-Kutta Method with its Adjoint One // Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics". — 2012. — No. 3. — P. 87-96. — (in russian). ]
24. Gevorkyan M. N., Gladysheva J. V. Symplectic Integrators and the Problem of Wave Propagation in Layered Media // Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics". — 2012. — No 1. — Pp. 50-60.
25. Logmani G. B. High Strong Order Implicit Runge-Kutta Methods for Stochastic Ordinary Differential Equations // System Dynamics Society. Proceedings of the 22nd International Conference. — Oxford, England, UK: 2004. — http://www. systemdynamics.org/conferences/2004/SDS_2004/PAPERS/109LOGHM.pdf.
26. Kinderman A. J., Mona,ha,n J. F. Computer Generation of Random Variables Using the Ratio of Uniform Deviates // ACM Transactions on Mathematical Software. — 1977. — Vol. 3, No 3. — Pp. 257260. — http://stevereads.com/papers_to_read/computer_generation_of_ random_variables_using_the_ratio_of_uniform_deviates.pdf.
UDC 519.21;51-76
Influence of Stochastization on One-Step Models
A. V. Demidova*, M. N. Gevorkyan*, A. D. Egorov^, D. S. Kulyabov*, A. V. Korolkova*, L. A. Sevastyanov*
* Telecommunication Systems Department Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, Russia, 117198 ^ Institute of Mathematics, NASB 11, Suruganov str., Minsk, Belarus, 220072
It is assumed that the introduction of probability in mathematical model makes it more adequate. There are practically no methods of the agreed (depending on structure of the system) introduction of probability in deterministic models. Authors have improved the method of constructing stochastic models for the class of one-step processes and illustrated it by models of population dynamics. Population dynamics was chosen for study because its deterministic models are sufficiently well explored that allows to compare the obtained results with the results already known.
We have examined the impact of the introduction of stochastics in the deterministic model, on the example of population dynamics system of type "predator-prey". Previously obtained stochastic differential equations are studied by the methods of the qualitative theory of differential equations. Stationary state and first integral of the system are obtained. To demonstrate the results the numerical simulations on the basis of Runge-Kutta method for stochastic differential equations are performed. The first integral of deterministic system (phase volume) in the stochastic case does not remain unchanged, but increases, which ultimately leads to the death of one or both populations.
One of the disadvantages of the classical system of type "predator-prey" is preservation of the amplitude of populations oscillations. In the stochastic model the process terminates with the death of one or both populations, which from the authors' point of view makes the model more adequate.
Key words and phrases: stochastic differential equations, "predator-prey" model, master equation, Fokker-Planck equation.
