Введение согласованного стохастического члена в уравнение модели роста популяций Текст научной статьи по специальности «Математика»

УДК 517.958

Введение согласованного стохастического члена в уравнение модели роста популяций

А. В. Демидова, Д. С. Кулябов

Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198

В работе рассматривается механизм получения стохастического дифференциального уравнения с согласованными стохастической и детерминистической частями.

Ключевые слова: стохастические популяционные модели, стохастическое дифференциальное уравнение, уравнение Фоккера—Планка.

1. Введение

Учёт случайных воздействий более адекватно позволяет описывать многие системы, такие как биологические, химические, экологические и т.д. Встаёт принципиальный вопрос о механизме ввода стохастических членов в детерминистическое уравнение. Обычно ввод осуществляется произвольным образом. Представляется более адекватным введение стохастических частей, согласованных с детерминистическими. Одним из возможных способов согласования является вывод стохастической и детерминистической частей из одного и того же уравнения. Наиболее удобным для исследования является стохастическое дифференциальное уравнение в форме Ланжевена, где стохастическая и детерминистическая части разделены. В свою очередь СДУ Ланжевена можно вывести из уравнения Фоккера-Планка, его, в свою очередь, можно получить по средствам разложения управляющего уравнения. Предметом изучения являются уравнения типа «хищник-жертва». Целью является записать их в виде СДУ и исследовать влияние стохастической части на поведение решения уравнения, описывающего попу-ляционную систему. Исследуются уравнения экспоненциального роста (уравнение Мальтуса) и модель «хищник-жертва».

2. Кинетические уравнения

Зачастую уравнения для популяционных моделей получают ad hoc, из разных пространных рассуждений. Однако представляется желательным получение этих уравнений из общих соображений. Одним из таких методов является метод химической кинетики [1].

В рамках этого метода рассматриваются процессы, состоящие из большого числа взаимодействий (чтобы можно было говорить о вероятности взаимодействия). При записи уравнений реакции используются концентрации1 объектов.

Пусть взаимодействуют два вещества: А и В. В результате получаем вещество Z. Реакция обратимая.

Запишем формальную схему этой реакции.

X + Y ^ Z. (1)

Статья поступила в редакцию 6 марта 2012 г.

1 Понятие концентрации употребляется здесь в расширенном смысле и зависит от моделируемой системы. Например, в популяционных моделях под концентрацией понимают число особей.

В рамках химической кинетики считается, что скорость взаимодействия веществ пропорциональна произведению концентраций этих веществ, т.е.

= к+Сх Су. (2)

Сх и Су — концентрации, коэффициент к+ называется абсолютной скоростью реакции и характеризует среднюю эффективность взаимодействия.

Количество взаимодействующих веществ задаёт порядок реакции, т.е. данная реакция есть реакция второго порядка.

Исходя из этих соображений записывают уравнения для концентраций.

-к+ху + к-г,

-к+ху + к-г, (3)

к+ху - к-г.

Положительные члены описывают прибыль соответствующих концентраций, а отрицательные — убыль.

3. Комбинаторная кинетика

Комбинаторная кинетика [2,3] расширяет методику химической кинетики. Она вносит следующие основные изменения:

— рассматриваются не скорости реакций, а вероятности перехода;

— детерминистические уравнения заменяются стохастическими уравнениями.

3.1. Уравнение Колмогорова-Чепмена

Достаточно часто поведение системы, которая содержит частицы или индивидуальные объекты, весьма правдоподобно описывается с помощью скачкообразного процесса.

Уравнение Колмогорова-Чепмена гласит, что вероятность нахождения частицы в точке х в момент I + т равна сумме вероятностей всех возможных смещений А из положениях х + А, умноженной на вероятность нахождения в точке х + А в момент ¿. Это предположение основано на независимости скачка А от какой-либо предыстории движения: необходимо знать начальное положение частицы только в момент ¿, а не в какие-либо предшествующие моменты. Это и есть постулат Маркова, и уравнение Колмогорова-Чепмена является основным динамическим уравнением всех марковских процессов.

Рассмотрим последовательные моменты времени <12 < ^ и марковский процесс У(¿), тогда пользуясь определением условной вероятности и марковского процесса можно получить

Р(У1,Ь; у2,г2] уз^з) = Р(У\М; У2М)р(узМу1М; У2М) =

= р(У1 ,^Му2МуииМузМу2(4)

Интегрируя по у2 получаем

^ дУ_

Ж ¿г

р(УъЬ; У2,Ь; уз^з)ёУ2 = Р(УъЬ; уз^з) =

= РЫ,^ Р(У2^21У1,Ь)р(Уз^з1У2^2)ёУ2- (5)

После деления на р(у1,Ь\) получим уравнение Колмогорова-Чепмена:

Р(У3^з\УъЬ) = I р(У3^з\У2^2)р(У2^2\У1,к)Ау2. (6)

Этому тождеству удовлетворяет вероятность перехода любого марковского процесса. Кроме того, уравнение также верно и в случае, если у является вектором с г компонентами. В том случае, когда у принимает только дискретные значения, интеграл заменяется суммой:

р (шЛК^з) = ^ р (пиЬ\П2М)Р (П2 М-пзМ). (7)

"2

3.2. Управляющее уравнение

Управляющее уравнение [4, с. 661], или основное кинетическое уравнение, описывает эволюцию распределения вероятностей в цепи Маркова с непрерывным временем и представляет собой разновидность уравнения Колмогорова-Чепмена для марковских процессов.

При соответствующих предположениях

limp(x,t + At\z,t) = ö(x - z) (8)

и

p(x,t + At\z,t) .

lim At = W (x\z>f) (9)

уравнение Колмогорова-Чепмена может быть сведено к управляющему уравнению:

dP (y,t)

dt

= (W(y\y')P(y',t) - W(y'\y)P(y,t))dy', (10)

где Ш(у\у') — вероятность перехода системы из состояния у' в состояние у за единицу времени. Если множество возможных значений является дискретным множеством состояний, то управляющее уравнение сводится к выражению:

^¡^ = р* ю - ™п>,прпт (11)

п'

Таким образом, основное кинетическое уравнение представляет собой уравнение баланса для вероятности каждого состояния п. Первый член соответствует возрастанию вероятности из-за переходов из других состояний п', а второй — уменьшению вероятности из-за переходов в другие состояния.

4. Процессы рождения—гибели

Широкий круг явлений можно моделировать специальным классом процессов, называемых процессами рождения-гибели. Название происходит из рассмотрения популяций людей или животных, в которых отдельные индивидуумы рождаются и умирают.

Одномерный прототип всех систем типа рождения-гибели представляет собой систему особей вида X, количество которых х может принимать неотрицательные целочисленные значения. Обычно принимается, что в ходе одного события может возникать (рождаться) или исчезать (погибать) лишь конечное число особей. В простейшем случае рождения и гибели особей X являются единичными событиями с зависящей от времени вероятностью, так что вероятность перехода

может быть записана в виде

№ (х№,*) = г+(х)'6х,х'+1 + Г(х)'6х,х'-1. (12)

Здесь имеют место два процесса:

х ^ х + 1 : Ь+(х) = вероятность перехода в единицу времени,

х ^ х - 1 : ¿- (х) = вероятность перехода в единицу времени. Тогда общее управляющее уравнение примет вид:

дгР(х,Цх',г') = г+(х - 1)Р(х - ) + г(х + 1)Р(х + ) -

- [г+(х) + г(х)] (13)

4.1. Многомерные процессы рождения-гибели

Эволюция во времени многомерных систем рождения-гибели может быть рассмотрена как результат индивидуальных взаимодействий между элементами некоторого множества. К таким системам можно отнести химические реакции (реакции взаимодействия молекул), экологические системы, где особи гибнут, дают потомство и уничтожают друг друга, системы, описывающие эпидемии, в которых заболевание передаётся от одного индивидуума другому при контакте, и т.д. Такие системы можно моделировать с помощью метода комбинаторной кинетики. В основе данного метода лежит предположение, что вероятность перехода из одного состояния в другое, являющегося следствием взаимодействия, пропорциональна числу возможных взаимодействий данного типа.

Рассмотрим систему из п компонентов, в которой происходят 8 различных взаимодействий:

£ ^ХакЛ £ МАХа (А = 1, 2,..., 8). (14)

а кА а

Коэффициент N А при Ха есть число компонентов типа Ха в левой части уравнения, а МА — соответственно в правой. Запишем в векторном представлении:

X = (Х1,Х2, ...,Хп),

^ = (ХА,ХА,...,ХА), (15)

ША =

где ха — число компонентов вида Ха. Также определим

И = ША - (16)

Таким образом, один шаг взаимодействия А в прямом и обратном направлениях соответственно можно записать как

х —> х + гА,

А (17)

х —> х - г .

Вероятности перехода в единицу времени пропорциональны соответственно числу способов выбора комбинации МА или МА из х компонентов и определяются

выражениями:

*А (х) = к+Л П (Ха -хАу

Л ха! (18)

Г-(х) = кЛ П (Ха -МЛ)!

Таким образом, общий вид однородного управляющего уравнения для цело-

Л

численной переменной х, изменяющейся шагами длины г Л, принимает вид: дР(*,€) ^ гг+_ , л

дЬ

= Т,{ + гЛ, *)Р(х + г Л,1) - *Л(х)Р(х, €)] +

+ ^Л(х - гЛ, г)Р(х- гЛ,г) -Г-(х)Р(х,г)]}. (19)

5. Стохастические дифференциальные уравнения

Стохастическое дифференциальное уравнение (СДУ) — дифференциальное уравнение, в котором один член или более имеют стохастическую природу, т.е. представляют собой стохастический процесс (другое название — случайный процесс) [5]. Таким образом, решения уравнения также оказываются стохастическими процессами. В физике СДУ традиционно записывают в форме уравнения Лан-жевена.

СДУ в форме уравнения Ланжевена состоит из обычного нестохастического дифференциального уравнения и дополнительной части, описывающей белый шум. Запишем нестрогое уравнение Ланжевена в виде:

ёх

- = фл+ыт. (20)

где х — интересующая нас переменная, а(х, Ь) и Ь(х, Ь) — некоторые известные функции, а — быстро флуктуирующий случайный член (т.е. при £ = значения и ) статистически независимы). Далее приведём интерпретацию этого уравнения в виде стохастического интегрального уравнения:

х(г) - х(0) = !д£а [х(1!),1!] + !ш(¿)Ь [х(г'),г']. (21)

0 0

Определение. Стохастическая величина х(1) подчиняется стохастическому дифференциальному уравнению Ито, записанному в виде

ёх(Ь) = а [х(Ь),Ь] + Ь [x(t),t]dW(Ь), (22)

если для всех £ и 10

х(г) = х(Ьо) + I [х(£] + I (^)Ь [х(£]. (23)

4о 4о

5.1. Уравнение Фоккера—Планка

Уравнение Фоккера-Планка, которое представляет собой уравнение в частных производных и описывает эволюцию плотности вероятности во времени. Уравнение Фоккера-Планка — одно из стохастических дифференциальных уравнений,

которое описывает временную эволюцию функции плотности вероятности координат и импульса частиц в процессах, где важна стохастическая природа явления. Названо в честь нидерландского и немецкого физиков Адриана Фоккера и Макса Планка, также известно как прямое уравнение Колмогорова.

Впервые уравнение было использовано для статистического описания броуновского движения частиц в воде. Для этого вводится функция плотности вероятности Ш(у, ¿), описывающая вероятность того, что частица имеет скорость в интервале (V, V + dv), если в момент времени 0 она имела начальную скорость Vo. Общая форма уравнения Фоккера-Планка для N переменных:

дШ

~дГ

N

д

N N

д2

1=1

дх;дхл г°

г=1 J = 1

Б2(х1, ..., хм)

W, (24)

где В1 — вектор сноса и В2 — тензор диффузии, причём диффузия вызвана действием сил стохастической природы.

Рассмотрим следующее стохастическое дифференциальное уравнение

¿х = а(х, г) ¿г + в(х, ^ато", (25)

где X £ — функция состояния системы, а W € — стандартное N-мерное броуновское движение.

Можно показать, что уравнение Фоккера-Планка для плотности условной вероятности р(х, £|хо, ¿о) = р имеет вид:

дР (х г)

Ш

N

д

N N

1=1

д2

• 1-1

г=1 ] = 1

В(х, ^В(х, t)J

(26)

5.2. Разложение Крамерса-Мойала

Разложение Крамерса-Мойала используется для перехода от управляющего уравнения к эквивалентному ему уравнению Фоккера-Планка [2].

Преобразуем управляющее уравнение (19) к уравнению Фоккера-Планка. Для этого делается несколько предположений. Во-первых, предполагается, что имеют место только малые скачки, т.е. Ьа(х) является функцией медленно изменяющейся с изменением х. Второе предположение говорит о том, что Р(х^) также медленно изменяется с изменением х. Тогда можно выполнить сдвиг из точки (х ± гл) в точку х, разложив правую часть в ряд Тейлора:

дР (х г)

ш

£

П=1

(гАу)п _

П\

*~Л(х)Р (х,()

+ £

п=1

(-гАУ)

п\

*а(х)р (х,г)

(27)

и отбросив члены порядка выше второго получим уравнение Фоккера-Планка:

(28)

дРI) = - £ ^ [Мх)р(х, г)]+1_ £ [ВаЪ(х)р(х, т ,

т

а.Ь

где

Аа(х) = ^ ^ №(х) - ГА(х)]

А

ваЬ(х) = £ г^ [г+(х) - г+(х)

п

6. Популяционные модели и комбинаторная кинетика

Популяционные системы можно рассматривать как системы взаимодействующих элементов. Таким образом, эволюция этих систем является результатом взаимодействия элементов. Процессы, происходящие в системе, можно описать с помощью схем взаимодействия следующего вида

£хЛха% £мЛха (Л = 1,2,...,з). (30)

а кА а

На основе этих схем можно записать управляющее уравнение, описывающее процессы, происходящие в системе, и с помощью разложения Крамерса-Мойала получить стохастическое дифференциальное уравнение в виде уравнения Фок-кера-Планка. Рассмотри простую систему популяционной динамики, описывающую экспоненциальный рост.

6.1. Модель экспоненциального роста

Экспоненциальный рост имеет место в тех популяциях, в которых прирост численности (число рождений минус число смертей) пропорционален числу особей популяции, а также делается предположение, что рост популяции происходит при неограниченных ресурсах и отсутствии гибели от болезней, хищников и т.п. Таким образом, в модели неограниченного роста возможны два процесса: рождение и гибель особи. Запишем схему для данной модели:

X 2Х

(31

X 0.

Первая строка описывает рождение новой особи с коэффициентом рождаемости к\, а вторая — гибель с коэффициентом смертности Вероятность перехода из состояния х в х' имеет вид:

№ (х\х',г) = (х)' 5Х>Х>+1 + Г(х)' 5Х>Х> Ль (32)

где

х! х!

Ь+(х) = к1~.-—т = к1 х, Ь+(х) = к-2--—гт = к2х. (33)

1 (х - 1)! 2 (х - 1)!

Все процессы необратимы, поэтому 1Л(х) = 0. Далее, пользуясь формулой разложения Крамерса-Мойала, получим УФП для модели неограниченного роста популяций. Коэффициенты щ:

а1 = - х')№'(х'\х) = к1х - к2х,

Х

а2 = - х')2'№(х'\х) = к1х + к2х.

(34)

Тогда УФП имеет вид:

яр (г +) й 1 В2

= - — ((к1Х - к2х)Р(х, 1)) + 2 ((к1Х + к2х)Р(х, I)), (35)

что эквивалентно СДУ в форме Ланжевена:

ёх(?) = (к1Х - к2х)<И + \/(к1Х + к2х)(Ж(г). (36)

6.2. Стохастическая модель «хищник—жертва»

Рассмотрим модель системы «хищник-жертва» [6-8], состоящей из животных двух видов, причём один из них охотится за другими, которые обеспечены неисчерпаемыми пищевыми ресурсами. Введя обозначения X — жертва, У — хищник, А — пища жертв, можно записать возможные процессы:

X + А -Ь 2Х, X + У 2У, У % В, (37)

которые имеют следующую интерпретацию. Первое соотношение означает, что жертва, которая съедает единицу пищи, немедленно репродуцируется. Второе соотношение описывает поглощение жертвы хищником и мгновенное репродуцирование хищника. Это единственная рассматриваемая возможность гибели жертвы. Последнее соотношение — это естественная смерть хищника.

Зададим распределение вероятности Р(х,у,Ь) для числа особей в заданный момент и найдём вероятностный закон для изменений (37). При этом предполагаем, что для бесконечно малого временного интервала А! имеют место следующие формулы для вероятностей переходов:

Р (х ^ х + 1; у ^ у) = к1ахА1,

Р(х ^ х - 1; у ^ у + 1) = k2xyАt, (38)

Р(х ^ х; у ^ у - 1) = к3уАг, Р(х ^ х; у ^ у) = 1 - (к1ах + к2ху + к3у)А1.

Далее получим уравнение Колмогорова-Чепмена. Для этого запишем вероятность в момент времени £ + А£ как сумму членов, каждый из которых вероятность предыдущего состояния, умноженная на вероятность перехода в состояние (х,у).

Р(х,у,г + А£) = к1й(х - 1)Р(х - 1,у,г)Аг + к2(х + 1)(у - 1) х х р(х + 1,у - 1,г)Аг + к3(у + 1)р(х,у + 1,г)Аг +

+ Р(х, у, Ь) - (к1ах + к2ху + к3у)Р(х, у, ¿)АЬ (39)

к1а(х - 1)Р(х - 1,у,£) + к2(х + 1)(у - 1) х

или

р(х,у,г + Аг) - р(х,у,г) АЪ

х Р (х + 1, у - 1,г) + к3(у + 1)Р (х, у + 1,г) -

- (к1ах + к2ху + кзу)Р(х,у,Ь), (40)

и, переходя к пределу А! ^ 0, получаем

9Р^^ = к1 а(х - 1)Р(х - 1, у, г) + к2(х + 1)(у - 1)Р(х + 1,у - 1,г) +

+ к3(у + 1)Р(х, у + 1, г) - (к1ах + к2ху + кзу)Р(х, у, г) (41)

Стоит отметить, что вероятностные законы (38) верны при условии марковости процесса. Для данной системы эти предположения выполняются до тех пор, пока различные особи одного вида можно считать идентичными. Конкретный вид переходных вероятностей принимается на различных основаниях и зависит от степени знания деталей рассматриваемых процессов рождения-гибели.

Запишем разложение Крамерса-Мойала для управляющего уравнения системы «хищник-жертва». Как уже отмечалось, в системе существуют следующие

процессы:

X + А ^ 2Хг1 = (1, 0), X + У ^ 2Уг2 = (-1,1), У ^ Вг3 = (0,-1). (42) Все процессы не обратимы, поэтому Ь-А(х) = 0, при А = 1, 2, 3, а

, , , х\ у\

£Г(х) = к1а--— — = к1ах,

(х - 1)\у\

г+(х) = к2 = к2Ху, (43) , , , х\ у\.

Ч(х) = к3 ^ = Ьу.

и управляющее уравнение имеет вид:

9Р^ У) = (к1а(х - 1)Р(х - 1,у) + к2(х + 1)(у - 1)Р(х + 1,у - 1) +

+ кз(у + 1)Р(х, у + 1) - (Ьах + к2ху + кзу)Р(х, у) (44) Воспользуемся для него разложением Крамерса-Мойала и получим:

дР<x,t) '< )п

Е (Цг^ <*Я

A ^ /

dt \ п

А,п

(45)

Отбросив члены порядка выше второго имеем уравнение Фоккера-Планка, соответствующее исходному управляющему уравнению:

Щ^И = - £ да <Аа<х, у)Р<х, у))+ 1 £ дадь <ВаЪ<х, у)Р<х, у)). (46)

а а,Ъ

где

Аа<Х,у)= ^ tA -x)^, ВаЪ <Х,у)= ^ ^ ' (47)

А=1,2,3 А=1,2,3

Таким образом имеем

(1 \ (—1\ / 0 \ fkiax — к2ху\

А<х) = ^ к,ах +{ i ) к^ху + ^—i) кзу = к2ху — кзу),

B<x)=fy<1 0) к\ах +( -1^-—1 1) к2ху + ( °1j<0 —1) кзу =

/к\ах + к2ху —к2ху \ = \ —к2ху к2ху + кзу)'

(48)

7. Заключение

Для описания эволюции популяционных систем в большинстве случаев используются детерминистские модели. Но они имеют свои недостатки, например не учитывают вероятностный характер процессов рождения и гибели индивидуумов популяций. Некоторые из этих недостатков позволяют решить стохастические модели, поэтому именно им в работе уделено основное внимание. Наиболее удобным для исследования является стохастическое дифференциальное уравнение в форме Ланжевена, так как в нем стохастическая и детерминистическая

части разделены. В работе был рассмотрен механизм получения стохастического дифференциального уравнения с согласованными стохастической и детерминистической частями.

В качестве уравнения, описывающего систему, предлагается использовать уравнение Колмогорова-Чепмена, из которого с помощью разложения Крамерса-Мой-ала можно получить уравнение Фоккера-Планка. Уравнение Фоккера-Планка является стохастическим дифференциальным уравнением, но в нем нет чёткого разделения на стохастическую и детерминистическую части, однако его можно привести к эквивалентному стохастическому дифференциальному уравнению в форме Ланжевена. Кроме того, в качестве примера рассмотрены модели экспоненциального роста популяции и модель «хищник-жертва» Лотки-Вольтерра.

Литература

1. Рубин А. Б., Пытьева Н. Ф., Ризниченко Г. Ю. Кинетика биологических процессов. — Изд-во Моск. ун-та, 1977. [Rubin A. B., Pihtjeva N. F., Riznichenko G. Yu. Kinetika biologicheskikh processov. — Izd-vo Mosk. un-ta, 1977. ]

2. Гардинер К. В. Стохастические методы в естественных науках. — Мир, 1986. — 526 с. [Gardiner K. V. Stokhasticheskie metodih v estestvennihkh naukakh. — Mir, 1986. — 526 s. ]

3. Ван-Кампен Н. Г. Стохастические процессы в физике и химии. — М.: Высшая школа, 1990. [Van-Kampen N. G. Stokhasticheskie processih v fizike i khimii. — M.: Vihsshaya shkola, 1990. ]

4. Паули В. Труды по квантовой теории. — М.: Наука, 1975. — Т. 1. [Pauli V. Trudih po kvantovoyj teorii. — M.: Nauka, 1975. — T. 1. ]

5. Оксендаль Б. Стохастические дифференциальные уравнения. Введение в теорию и приложения. — М.: Мир, 2003. [Oksendalj B. Stokhasticheskie differencialjnihe uravneniya. Vvedenie v teoriyu i prilozheniya. — M.: Mir, 2003. ]

6. Вольтерра В. Математическая теория борьбы за существование. — Наука, 1976. [Voljterra V. Matematicheskaya teoriya borjbih za suthestvovanie. — Nauka, 1976. ]

7. Lotka A. J. Elements of Physical Biology. — Baltimore: Williams and Wilkins, 1925.

8. Ризниченко Г. Ю., Рубин А. Б. Математические модели биологических продукционных процессов. — Изд-во МГУ, 1988. [Riznichenko G. Yu., Rubin A. B. Matematicheskie modeli biologicheskikh produkcionnihkh processov. — Izd-vo MGU, 1988. ]

UDC 517.958

Introduction of Self-Consistent Term in Stochastic Population

Model Equation

A.V. Demidova, D. S. Kulyabov

Telecommunication Systems Department Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, Russia, 117198

In this paper we consider a mechanism to obtain stochastic differential equation with a consistent stochastic and deterministic parts.

Key words and phrases: stochastic population model, stochastic differential equation, Fokker-Planck equation.