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

\(Ш))

I \QQ/ *

УДК 519.6 551.7

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ТРЕХ-ПОПУЛЯЦИОННОЙ СИСТЕМЫ

Суханов А. И., Зайцева Е. Ю.

Донской государственный технический университет, Ростов-на-Дону, Российская федерация sukhinov@gmail .com katerinazaytseva77@gmail.com Рассматривается модель динамики трех-популяционной системы, включающей два вида жертв, отличающихся друг от друга скоростью размножения и показателями максимально возможной плотности, и один вид хищника. Рассматривается нуль-мерная модель (пространственно-однородная),

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

UDC 519.6 551.7 MATHEMATICAL MODELING OF THE DYNAMICS OF THREE-POPULATION SYSTEM

Sukhinov A. I., Zaytseva E. Yu.

Don State Technical University, Rostov-on-Don, Russian Federation

sukhinov@gmail .com katerinazaytseva77@gmail.com In this paper, we consider a model of the dynamics of a three-population system, including two species of prey, differing from each other in their reproduction rate and indicators of the maximum possible density, and one predator species. A zero-dimensional model (spatially homogeneous), formulated as a Cauchy problem for a system of three ordinary differential equations with quadratic nonlinearities in the right-hand sides, is considered. For this model, its discrete analogue is built - an explicit difference scheme, having the first order of accuracy, relative to the time grid step. The limitations on the time step that ensure the stability of the constructed difference scheme are determined. Numerical experiments were performed for various modes of model operation: a monotonous exit to the stationary mode of individual populations, a monotonous increase in the populations of prey, a highly oscillatory mode of the number of hare populations, an oscillatory change in the number of prey and predators with different characteristic times of oscillation and the establishment of a stationary mode.

Ключевые слова: хищник, жертва, система, Keywords: predator, prey, system, model, модель, схема, плотность, популяция, scheme, density, population, number, mode. численность, режим.

Введение. Математические модели популяционной динамики давно являются предметом исследования специалистов в области математической биологии, системного анализа, математического моделирования [1]. Первая из известных моделей типа «хищник-жертва», носящая имя Т. Р. Мальтуса, была построена более трехсот лет назад [2]. За прошедшее время

г(шУ

Iл

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

В данной работе рассматривается модель динамики трех-популяционной системы, включающей два вида жертв, отличающихся друг от друга скоростью размножения и показателями максимально возможной плотности, и один вид хищника. Рассматривается нульмерная модель, формулируемая как задача Коши для системы трех обыкновенных дифференциальных уравнений с квадратичными нелинейностями в правых частях. Для данной модели построен ее дискретный аналог — явная разностная схема, имеющая первый порядок точности, относительно шага временной сетки [3]. Определены ограничения на шаг по времени, обеспечивающие устойчивость построенной разностной схемы.

Выполнены численные эксперименты для различных режимов функционирования модели:

- монотонного выхода на стационарный режим численностей отдельных популяций;

- монотонного увеличения популяций жертв;

- сильно-колебательного режима численности популяции зайцев.

- колебательного изменения численности жертв и хищников с различными характерными временами колебаний и установления стационарного режима.

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

Описание модели. Математическая постановка задачи. Рассматривается биологическая система, состоящая из трех видов популяций, которые условно назовем зайцы — их численность обозначается символом Ъ = Z(t), косули, численность которых есть К = К(1:) и популяция волков, численностью Ш = где t — время, t > 0. Первые две из перечисленных популяций будем называть жертвами, а волков — хищниками. Для жертв будем принимать во внимание процесс естественного размножения и их убыли вследствие уменьшения численности ввиду наличия хищников, а также максимально возможную численность, обусловленную наличием корма и биологической емкостью среды для каждого из видов. Скорость размножения жертв будем считать пропорциональной их численности с коэффициентом, зависящим от отклонения их численности от максимально возможной, с учетом биологической емкости среды. Естественную смертность жертв учитывать не будем. Скорость убывания жертв пропорциональна численности хищников и численности жертв. Скорость размножения хищников предполагается пропорциональной численности каждого вида жертв. Убыль хищников происходит за счет их естественной смертности, которая в данном случае пропорциональна численности волков.

С учетом сделанных предположений, изменения численности трех видов популяций описывается следующей системой обыкновенных дифференциальных уравнений (ОДУ), (1)-(3):

= а'з(Ш) • (Уэ • К + у2 • Ъ) - уз • Ш

(1)

(2) (3)

где:

где К, Ъ — число жертв (косули, зайцы); Ш — число хищников (волки); ах, а2, а3 — коэффициенты размножения косуль, зайцев, волков; рх, р2 — коэффициенты уменьшения численности косуль, зайцев; ух, у2 — коэффициент усвоения волками пищи (косуль и зайцев); у3 — коэффициент уменьшения(смертности) численности волков; К тах, Ътах, Штах — максимальные значения численности популяций косуль, зайцев и волков; К уравнениям (1)-(3), для выделения единственного решения, необходимо добавить соответствующие начальные условия (4)-(6).

Я^о) = (4)

гао) = Вх, (5)

Ж^) = С!, (6)

где — начальный момент времени.

Аналитическое исследование математической модели. Поскольку система ОДУ является нелинейной, содержит произведения численностей популяций в правых частях, то получение решения в замкнутом виде, как некоторых конечных формул, является весьма затруднительным [4]. С точки зрения исследования поведения этой системы отметим, что правые части системы ОДУ являются достаточно гладкими функциями относительно переменных, которыми являются численности отдельных видов популяций: производные первого и второго порядков существуют и являются непрерывными функциями, а производные более высоких порядков обращаются в 0 [5, 6]. Для краткости изложения отметим, что условия существования и единственности решения задачи (1)-(6) выполняются. Установим условия, при которых решение является стационарным, т.е. не зависит от времени. Для этого приравняем правые части уравнений (1 )—(3) к 0. Получим систему (7)-(9):

К • Ш = 0

1тах 2

-= 0

(7)

(8) (9)

+ 72^)-Кз^ = 0.

Получить в общем виде решение нелинейной системы (7)-(9) представляется затруднительным, т.к. она, в общем случае, сводится к алгебраическому уравнению шестого порядка, для решения которого общая формула не существует. Будем искать равновесные значения одним из несложных итерационных методов — методом половинного деления, учитывая тот факт, что, исходя из биологического смысла задачи, каждое из равновесных значений Ъ, К, W — является положительным числом и не превосходит, соответственно, значения:

^тах = 1000;

■'тах

Аналитическое исследование

= 10000; Жтах = 100; математической

модели. После 10 шагов методом

половинного деления получаем значения равновесных концентраций отдельных видов популяций с точностью не хуже 0,1%, что достаточно для прогностического моделирования. Различные режимы функционирования модели будут определяться на основе задания начальных условий (4)-(6), которые могут отличаться от равновесных значений концентраций, как в меньшую, так и в большую сторону, и даже превосходить в начальный момент их максимально возможные значения:

К(*;0) « 505;

7(Г0) « 8748;

Ж(^) « 88;

Построение разностной схемы и алгоритма вычислений. Будем предполагать, что решение поставленной задачи Коши (1)-(6) требуется найти для Т > I >0. На данном временном промежутке построим равномерную временную сетку [7]:

т

to

0 tl t2... tn t„-l...T

► t

Рис. 1. Временная сетка: т — шаг по времени, Т — максимальное значение времени, 10. — начало

отсчета (начальный момент времени), = 0

Явная схема Эйлера для поставленной задачи имеет вид: Для уравнения (1):

йК

аг

Для уравнения (2):

dZ dt

Для уравнения (3)

K(tn+l)-K(tn) — ai * K(tn) *(l -¡¡.Щ-^* K(tn) * W (tn);

~ V Лmax'

Z(tn+1)-Z(tn) — a2 * Z(tn) *(l -Ш)-р2* Z(tn) * W(tn);

~ V ¿max'

dW

dt

W(tn+i) - W(tn)

= a3(W) • (Yi • K(tn) + Y2 • Z(tn)) - Ys • W(tn);

Выразим:

K(tn+i) — Kn+i;

K(tn+i) — K(tn) + ai*T* K(tn) * (l - Pi * т * K(tn) * W(tn);

^ Лmax'

Kn+i — Kn + ai*T*(l-^)-Pi*T*Kn*Wn;

^ ^max'

Z(^n+i) — Zn+i;

Z(tn+i) — Z(tn) + a2*T* Z(tn) * (l - -P2*T* Z(tn) * W (tn);

V ¿max'

Zn+i — Zn + a2*T*(l - ) - 02

V ¿max'

*x*Zn*Wn;

W(tn+i) — Wn+i;

W(tn+i) — W(tn) + a3 (W) * (yi * K(tn) + Y2 * Z(tn)) -Ys*** W(tn); Wn+i — Wn + a3(W) * (Yl *Kn + Y2*Zn) - Y3*T * Wn; Получены расчетные формулы для каждого из уравнений (10)-(12):

Kn+i — Кп + a1*T*(l-Jjb-)-^

^ ^max'

Zn+i — Zn + a2*T*(l - ) - 02

V ¿шах'

*T*Kn*Wn ;

*x*Zn*Wn ;

Шп+1 = Шп + •т• (у^Кп + У2^п) -Уз^^п ;

Введем значения коэффициентов: Максимальные значения для каждой популяции:

Кщах = I00,

гтах = 10000;

Щпах = 100; Начальные значения популяций (4)-(6):

К(10)= 30;

(10) (11) (12)

2(*о)= 50; Ж (^)= 20;

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

Приближенное определение а1 (при р1 =0) из уравнения (1):

^ ~ м ' Т = ДГ « 8766 часов;

1 к

2Л ( К \ 4- = -—);

1 Л ^ 1 1 = *1*(1 - 4) =47? = 4;

«1 « - « 0,5 * 10-3; Приближенное определение в 1 (при а1 =0) из уравнения (1):

А*^* |ДК| ДГ = ^Тр« 103 Жтах = 100 = 10

^1тт =

00005

_ 1 /1 1 \ 1 /1 1 \ 1 _ = 4ТД7 * (Щпах + Й^"/ = 2 * 103 * (10 + 100) = 2 * 104 = 0,

Приближенное определение а2 (при в 2=0) из уравнения (2): 1

47 7 — = «2*^* (1 -)

1 Л ^ 1 1 4*7 = а2*(1 - 4) = 47? = а2* 4 1

«2 « - « 0,5 * 10-2

Приближенное определение в2 (при а 2=0) из уравнения (2):

|Ду = |Д2| ДГ = 102 = 100 = 10

^2тт =

_ 1 /1 1 \ 1 /1 1 \ 1 _ = 4ТД7 * щЗ + ж„/„) = 4 * 102 * (10 +100) = 4 * 103 = 0,0005;

Приближенное определение а3: а3=1;

Приближенные определения у1, у2, у3: у1=0,1; у2=0,002; у3=0,1.

Определение допустимого шага по времени. Введем безразмерную временную переменную по формуле.

t

= г

Т = 365 • ^ суток • 44 часа = 8766 часов

ж ¿к

—— = • ^ = 0<^<1 х«1

С учетом новой переменной система (1)-(3) примет вид:

^ _ ^ _ ^ 1 ^ = йГ dt = йГ Т

£==л;

£=М-(1 =/2;

^ = (а'з • (У1-К + У2 • г) - Гз • И0 'Г = /з;

Используя полученные соотношения, определим допустимую величину временного шага.

^ 1 ^ 1 1 Т1 < ¡э/ТГ, т2 < ¡э/ТТ, тз < '

(14)

(15)

I ЭК I 0<С<1

I эг I 0<с'<1

|Э/3| 0<С<1

Т1 <

т; (16)

0<с'<1

^2 <

г,, тах|7(с')|+2 ,а шах|г({')| + „ тах|^')к; (17) ( 2 0<С'<1 2 2тах Р2 0<С'<1 )

Тз <

Т • (а'зГз-таХ|^Г (С')|)' ( 3 0<С'<1 )

(18)

Допустимое значение шага должно быть:

Т < ш/п(т1, т2, тз).

Проведение вычислительных экспериментов для различных режимов функционирования модели. Монотонный выход на стационарный режим численностей отдельных популяций (с насыщением для популяции волков) изображен на рис. 2.

Значения параметров:

Рис. 2. Монотонный выход на стационарный режим

а1 = 1,7; а2 = 3,5;

= 1;

№

«з = «з '(1 -)

& = 0,01; & = 0,005; 71 = 0,1; 72 =0,002; Уз = 0,1; г = 0,045.

Монотонное увеличение популяций жертв изображено на рис. 3.

Рис. 3. Монотонное увеличение популяции жертв

Значения параметров:

а1 = 1,7; а2 = 3,5; из = 1;

а' 3 = а3;

ßi = 0,01; ß2 = 0,005; y1 = 0,0001; у2 =0,005; Уз = 0,1; т = 0,025.

Колебательный режим численности популяций с выходом на стационарный режим представлен на рис. 4.

е

:

1 /

еМ . 1 ±1

_ Волки

Зайцы

— Косули

Значения параметров:

Рис. 4. Колебательный режим с выходом на стационар

а1 = 1,7; а2 = 3,5; аз = 1; а'з = аз;

ßi = 0,01; ß2 = 0,005; Yi = 0,3; Y2 =0,0001; Уз = 0,1; т = 0,025. k(t0) = 2000; Z(t0) = 20000;

ВД = 50.

Сильно-колебательный режим численности популяции зайцев продемонстрирован на рис.

0

Ы.1 .1

Волки Зайцы Косули

Рис. 5. Сильно-колебательный режим

Значения параметров:

а1 = 1,7; а2 = 3,5;

аз = 1; а'з = аз;

в1 = 0,01; Р2 = 0,005; У1 = 0,3; У2 =0,0001; Уз = 0,009; т = 0,045. k(t 0) = 4000; Ъ(1:0) = 40000; ВД = 50.

Колебательный режим численностей всех популяций изображен на рис. 6.

I-

— Волки

Значения параметров:

Рис. 6. Колебательный режим численностей всех популяций а1 = 1,7; а2 = 3,5;

аз = 1; а'з = аз;

в1 = 0,01; Р2 = 0,005; У1 = 0,3; У2 =0,001; Уз = 0,009; т = 0,045.

г(шУ

Iл

k(t0) = 2000; 2(г0) = 20000; кОю) = 50.

Заключение. В ходе данной работы были решены следующие задачи:

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

- определено положение равновесия соответствующей динамической системы;

- построена явная разностная схема Эйлера для численного решения соответствующей задачи Коши;

- выполнена оценка допустимо возможного шага по времени разностной схемы.

- на основе вычислительных экспериментов выявлены пять, содержательно отличных друг от друга, режимов функционирования динамической модели для системы с тремя видами популяций.

Библиографический список

1. Самарский, А. А. Математическое моделирование: идеи, методы, примеры / А. А. Самарский, А. П. Михайлов. — 2-е изд., испр. — Москва : Физматлит, 2001. — 320 с.

2. Моделирование сложных систем. Часть 1 / А. И. Сухинов [и др.]. — Ростов-на-Дону: ДГТУ-Принт, 2019. — 241 с.

3. Ильин, П. П. Математическое моделирование, часть 1: Непрерывные и дискретные модели / П. П. Ильин. — Новосибирск : Изд-во СО РАН, 2017. — 423 с.

4. Базыкин, А. Д. Математическая биофизика взаимодействующих популяций / А. Д. Базыкин. — Москва : Наука, 1985. — 181 с.

5. Вольтерра, В. Математическая теория борьбы за существование / В. Вольтерра; [пер. с франц. О. Н. Бондаренко]; под ред. и послесл. Ю. М. Свирежева. — Москва : Наука, 1976. —

6. Базыкин, А. Д. Нелинейная динамика взаимодействующих популяций / А. Д. Базыкин. — Москва-Ижевск : Институт компьютерных технологий, 2003. — 368 с.

7. Сухинов, А. И. Решение задач вычислительной математики на ЭВМ : учебное пособие / А. И. Сухинов, А. Е. Чистяков, Е. А. Проценко. — Ростов-на-Дону : ДГТУ-Принт, 2019. — 120 с.

286 с.