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

Вестник СамГУ — Естественнонаучная серия. 2009. № 6(72) 139

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 517.96: 574.34

ДВ-МОДЕЛЬ СИСТЕМЫ ’’ХИЩНИК-ЖЕРТВА”

© 2009 В.В. Зайцев, А.В. Карлов-младший, С.С. Телегин1

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

Ключевые слова: модель Вольтерра, запаздывающие связи, автоколебания, дискретное время.

Введение

Математическое моделирование процессов в био- и экосистемах является одним из основных методов их теоретического исследования. Знаменитая модель Вольтерра [1], описывающая систему ”хищник-жертва”, известна уже более семидесяти лет. Тем не менее, исследования в рамках различных модификаций этой классической модели активно проводились и в последнее время [2]. Известны недостатки модели, главный из которых — ее консервативность, не позволяющие проводить количественных оценок характеристик процессов в системе ”хищник-жертва”. Согласно сложившимся в последние десятилетия представлениям о динамике системы ее основным движением являются устойчивые автоколебания численностей видов [3]. Следовательно, при формулировке динамической модели система должна рассматриваться как автоколебательная.

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

1Зайцев Валерий Васильевич (zaitsev@ssu.samara.ru), Карлов-младший Александр Владимирович, Телегин Сергей Сергеевич, кафедра радиофизики и компьютерного моделирования радиосистем Самарского государственного университета, 443011, Россия, г. Самара, ул. Акад. Павлова, 1.

обратной связи кардинально изменять характер динамической системы [4]. Динамические модели системы ”хищник-жертва” с запаздыванием предложены и исследованы в работах Колесова [5] и Вангерски и Каннингема [6]. В математическом плане модели с запаздыванием представляют собой системы нелинейных дифференциальных уравнений с запаздывающим аргументом. Они не имеют точных аналитических решений, и при анализе моделей речь идет в основном о численных результатах. Ситуация еще более усложняется с введением в динамические модели флуктуаций параметров [7] и случайных внешних воздействий для описания реально существующих в системе стохастических процессов [8]. Для расчета вероятностных характеристик стохастических автоколебаний не удается воспользоваться методами теории марковских процессов, т. к. в данном случае из-за наличия в системе последействия предположение о марковости неправомерно.

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

1. Дифференциальная модель системы.

Стационарные режимы и их устойчивость

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

N = „N1 - ^NN2 - Рт (N1(1 - в)),

^ = -»2^ + 72N1 (г - т) N2 (г - т) .

Здесь V!, ^2, 71, 72 — положительные постоянные: V! — коэффициент прироста (репродуктивности) жертв, V2 — коэффициент вымирания хищников, 71 — коэффициент гибели жертв из-за встречи с хищниками, 72 — коэффициент размножения хищников. Запаздывание Т имеет смысл усредненного интервала времени между моментом гибели одной особи жертвы и моментом соответствующего увеличения числа взрослых хищников. Аналогично интервал в учитывает запаздывающую реакцию жертв на изменение условий поступления энергии из окружающей среды. В первом уравнении слагаемое Рт (N1^ - в)), где Рт — полином четной степени т, описывает внутривидовую конкуренцию у жертв, ограничивая рост их численности в

отсутствие хищников. Заметим, что при в = 0 и Рт(^) = ^2(^1) = в^і уравнения (1) описывают модель Вангерски — Каннингема [6], а при 71 =

= 0 и Рт(^і) = р2(^і) = в^і(£)Жі(і — в) первое уравнение системы (1) —

модель Ферхюльста-Перла [10] с запаздыванием.

В уравнениях (1) целесообразно перейти к нормированным численностям видов:

Пі(г) = 72Жі(і) , Н2(І) = 71 N2(і).

^2 ^і

После чего система (1) примет вид

= ^і«і — ^іНіП2 — Рт (пі(г — в)) , (2)

^ПГ = —»2П2 + ^2Пі (І — 1) П2 (і — Т) .

Здесь рт (пі (і — в)) = 72Рт (и2пі(і — в)/72) /у2- Дальнейшее обсуждение модели проведем для частного случая в = 0, полагая, что рт (пі(і — в)) =

= —пт (і).

Система уравнений (2) имеет стационарное решение

піо = 1, П20 = 1 — — = 1 — Ці, (3)

^і

описывающее состояние системы с неменяющимися во времени численностями видов.

Введем в рассмотрение новые переменные У і (І) и У2(г) — отклонения численностей от стационарных значений:

Пі(і) = піо + Уі(і), П2(і) = П20 + У2(і).

Для переменных у і (і) и У2(і) уравнения (2) можно записать в виде

^ = ——(т — 1)Уі — иіУ2 — иі Уі У2 — —Рт (Уl), (4)

= —»2У2 + ^2 (1 — —і) Уі(і — Т) + ^2У2(І — т) + ^2Уі(І — Т)У2(І — Т) ,

где рт(уі) = (1 + Уі)т — (1 + тУі) — нелинейная функция.

Устойчивость стационарного состояния (3) исследуется на основе анализа поведения решений вида Уі(і) = Аі ехр(ві), У2 (і) = А2 ехр(ві) линеаризованной системы уравнений (4). Этой системе соответствует характеристическое уравнение

в2 + (—(т — 1) + ^2) 8 + — (т — 1)^2 + (^2 — —т^2 — ^в) ехр(—вТ) = 0. (5)

Комплексные корни в = вг +]ві уравнения (5) удается найти лишь численно. Устойчивость обеспечивается выполнением условия вг ^ 0. Для отыскания границ областей устойчивости в пространстве параметров ^і,^2,— можно воспользоваться методом В-разбиений Неймарка [9].

На рис. 1 показано сечение поверхности вг = 0 в пространстве координат (иі,^2,—) плоскостями ^іТ = 0,25, ^іТ = 0, 5 и ^т = 0, 75. Ниже

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

8т

0.04

0.03

0.02

0.01

Стац. состояние У ¡1 =0.75

ур = ( ).5

У,Т = 0.25

Предельный цикл

0.1

0.2

0.3

0.4

0.5

V д

Рис. 1. Области устойчивого стационарного состояния и автоколебательного режима

Что же касается нулевого решения (що = 0,П20 = 0) системы (2), то его устойчивость характеризуется линеаризованными уравнениями

йДщ

= ^іДщ,

йДи2

= —^2 ДП2.

( (И

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

2. Уравнения движения системы в дискретном времени

В работах [10, 11] предложен способ перехода к дискретному времени в нелинейных динамических системах и формирования алгоритмов генерации ДВ-автоколебаний. Следуя [11], выделим в системе (4) устойчивую диссипативную подсистему, записав уравнения следующим образом:

^ = —-(т — 1) уі — VIу2 + fl(t,Уl,У2), (6)

^ = —^У2 + Ж^УъУ2).

При этом введены обозначения

Ь^^ЪУ2) = —^Уі(%2(^ — Ц<Рт (Уі^)) , (7)

f2(t,Уl,У2)= V2 (1 — -і) Уі(І — т) + V2У2(í — Т) + V2Уl(t — Т)У2^ — Т).

Формально считая выделенные линейные части уравнений (6) уравнениями состояния некой линейной системы с двумя независимыми входами и двумя выходами, введем в рассмотрение матричную системную функцию

в в

Н(в) = ( 5_0-51 5-51 1 3-32 \ (8)

\ 0 5-52 /

с полюсами ві = —ц,(ш — 1) и в2 = —^2; В = у1/ (ц,(ш — 1) — и2).

Для перехода к дискретному времени с интервалом дискретизации А предлагается применить принцип инвариантности импульсных характеристик, широко используемый в практике синтеза линейных цифровых фильтров [12]. В рамках такого подхода в элементах системной функции НВ-системы (8) проводится замена простых дробей вида

1А

в — ві 1 — ЄХр(віА)г-1’

В результате получается системная функция ДВ-системы

( 1 в________________В \

Н(г) = А ( 1-ехР(*1Д)^ 1 1—ехр(«1Д)-г 1 1 1—ехр(«2Д)-г 1 | (д)

\ 0 1—ехр(«2Д).г~1 у

Учитывая, что слагаемое с аргументом г-1 в системной функции соответствует временной задержке А в уравнении движения, по элементам матрицы (9) можно записать следующие уравнения движения ДВ-системы:

У(1)[п] = а111)у(1)[п - 1] + Ь011)/1[и,у1,у2],

У(2)[п] = а!12)У(2)[п - 1] + а212)у(2)[п - 2]+ &112)/2[п - 1,Уl,У2], (10)

У1[п] = У^Н + У(2)[п],

У2[п] = а122)У2[п - 1] + Ь022) /2 [п, У1, У2] •

Коэффициенты в уравнениях (10) связаны с параметрами НВ-системы соотношениями

,(11) _ехр(„ А) „ (12) =„хр(„ А) , „хр(„ А) „ (22)

аі =ехр(в1А), а1 = ехр(в1А) + ехр(в2А), а! = ехр(в2А),

а212) = — ехр(в1А + в2А),

Ь0П) = А, &112) = АВ (ехр(в1А) — ехр(в2А)), = А,

а нелинейные функции дискретного временного аргумента в соответствии с выражениями (7) имеют вид

¡1[п,У1,У2] = —^1У1[И]У2Н — №т (У1 М) ,

к[п,У1,У2] = ^ (1 — Ц1) У1[п — N] + У2У2[и — N + ^У^П — Щу2[п — Щ.

(11)

При этом предполагается, что интервал дискретизации составляет целую часть времени запаздывания: т = N А.

Система уравнений движения (10) содержит лишь одно нелинейное разностное уравнение (первое соотношение). К нему можно применить метод

1

последовательных приближений. Все остальные соотношения являются рекуррентными формулами, что обеспечивает высокую вычислительную эффективность дискретной модели (10).

Пример результатов, полученных в рамках модели (10), представлен на рис. 2: график процесса возбуждения автоколебаний (а) и их амплитудный спектр (б). При этом параметры модели имеют следующие значения: А = = т/10, v\r = 0, 5, V2T = 0,1, ßr = 0, 003, m = 4. Сплошная линия соответствует автоколебаниям жертв Y\[n] = yi[n] + nio, пунктирная — хищников Y2[n] = y2 [n] + П20. Частота на графике спектра измеряется в единицах частоты Q = 2п/т.

Рис. 2. Процесс возбуждения (а) и амплитудный спектр (б) автоколебаний

численностей видов

Как видно из графиков, в рассматриваемом примере в системе возбуждаются существенно нелинейные одночастотные автоколебания, характеризующиеся высоким уровнем гармоник основной частоты. Наличие постоянной составляющей в спектре автоколебаний означает, что закон сохранения средних, имеющий место в классической модели Вольтерра [1], в данном случае не выполняется.

3. Стохастическая модель со случайным запаздыванием

Динамическую модель (10) нетрудно дополнить шумовыми источниками, моделирующими случайные аддитивные и мультипликативные воздей-

ствия на систему. Стохастическая модель системы при наличии случайных изменений коэффициентов репродуктивности жертв VI и хищников ^2 приведена в [13]. Здесь мы рассмотрим автоколебания в системе со случайным запаздыванием. Отметим, что вопрос о наличии в системе ”хищник-жерт-ва” распределенного во времени запаздывания также обсуждался Вольтер-ра в монографии [1].

Будем считать, что дискретное время запаздывания N в уравнениях движения (10) является элементом временной последовательности случайных целых чисел с заданным распределением вероятностей Р(Ы) в интервале значений от Ыт[п до Ытах. Тогда движение системы будет представлять собой реализацию автоколебаний со случайным запаздыванием.

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

На рис. 3 представлен отрезок реализации случайного дискретного времени запаздывания, полученной путем выборки с Ыт[п = 2, Ытах = 18 из гауссова белого шума со средним значением < тп >= 10 и дисперсией = = 100 и последующей повышающей дискретизации (интерполяции) с коэффициентом I = 10. Последовательность Ып характеризуется выборочными средним Ыа = 10 и дисперсией = 21, 9. Нормированная корреляционная функция к(п) = К(п)/аN последовательности показана на рис. 4.

т

Рис. 3. Реализация случайного запаздывания

Рис. 4. Корреляционная функция случайного запаздывания

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

График процесса возбуждения автоколебаний представлен на рис. 5, а, усредненные амплитудные спектры автоколебаний У\(£) и у 2 (^) —

на рис. 5, б. Спектры рассчитаны методом Бартлетта по отрезкам реализаций установившихся автоколебаний: число отсчетов в реализации равно 16384, число реализаций в псевдоансамбле — 8. Фазовый портрет системы в координатах У\[п\ и У2[п] показан на рис. 6.

О 0.2 0.4 0.6 0.8 1

Рис. 5. Процесс возбуждения (а) и амплитудный спектр (б) автоколебаний численностей видов в модели со случайным запаздыванием

о--------1------1-------1------

О 0 5 1 15 . •

Рис. 6. Фазовый портрет автоколебаний численностей видов в модели со случайным запаздыванием

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

Заключение

Описанный здесь способ формирования дискретной временной модели нелинейной системы путем выделения линейной диссипативной подсистемы и использования принципа импульсной инвариантности применим к различным модификациям модели Вольтерра (А.Д. Базыкина [2], Ю.С. Колесова [5] и др.). Отметим также, что вольтерровские модели находят применение в формирующемся направлении современной экономики, ставшим известным в последнее десятилетие под названием ”Эконофизика” [15] или ”Динамическая экономика”.

Литература

[1] Вольтерра В. Математическая теория борьбы за существование. М.: Наука, 1976. 288 с.

[2] Базыкин А.Д. Нелинейная динамика взаимодействующих популяций. М.: Институт компьютерных исследований, 2003. 368 с.

[3] Мари Дж. Нелинейные дифференциальные уравнения в биологии. Лекции о моделях / пер. с англ. М.: Мир, 1983. 400 с.

[4] Рубаник В.П. Колебания квазилинейных систем с запаздыванием. М.: Наука, 1969. 288 с.

[5] Колесов Ю.С. Математические модели экологии // Исследования по устойчивости и теории колебаний. Ярославль: Изд-во ЯрГУ, 1979. С. 3-40.

[6] Бабский В.Г., Мышкис А.Д. Математические модели в биологии, связанные с учетом последействия // Нелинейные дифференциальные уравнения в биологии. М.: Мир, 1983. 400 с.

[7] Музычук О.В. Вероятностные характеристики системы ”хищник-жерт-ва” со случайно изменяющимися параметрами // Изв. вузов. Прикладная нелинейная динамика. 1997. Т. 5. № 2. С. 80-86.

[8] Гардинер К.В. Стохастические методы в естественных науках. М.: Мир, 1986.

[9] Эльсгольц Л.Э., Норкин С.Б. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. М.: Наука, 1971. 296 с.

[10] Зайцев В.В., Давыденко С.В., Зайцев О.В. Динамика автоколебаний дискретного осциллятора Ван дер Поля // Физика волновых процессов и радиотехнические системы. 2000. Т. 3. № 2. С. 64-67.

[11] ДВ-осцилляторы, порождаемые томсоновскими автоколебательными системами / В.В. Зайцев [и др.] // Физика волновых процессов и радиотехнические системы. 2008. Т. 11. № 4. С. 98-103.

[12] Оппенгейм А., Шафер Р. Цифровая обработка сигналов. 2-е изд. М.: Техносфера, 2006. 856 с.

[13] Зайцев В.В., Телегин С.С. Интегральная модель автоколебаний в системе ”хищник-жертва” // Физика волновых процессов и радиотехнические системы. 2009. Т. 12. № 2.

[14] Hall C.A.S. An assessment of several of the historically most influential theoretical models used in ecology end of the data provided in their support // Ecological modeling. 1988. V. 43. P. 5-31.

[15] Романовский М.Ю., Романовский Ю.М. Введение в эконофизику. Статистические и динамические модели. М.; Ижевск: Институт компьютерных исследований, 2007. 280 с.

Поступила в редакцию 17//V/2009;

в окончательном варианте — 17//V/2009.

THE DISCRETE TIME ”PREDATOR-PREY” MODEL

© 2009 V.V. Zaitsev, A.V. Karlov-junior, S.S. Telegin2

Method of impulse invariance, which is often used to design digital filters, was used to create the discrete time Volterra model with lagging.

The predator-prey system was studied as a discrete self-oscillating system.

The results of simulation modeling of a system with random lagging are listed.

Key words: Volterra model, lagging links, self-oscillations, discrete time.

Paper received 17/IV/2009. Paper accepted 17/IV/2009.

2Zaitzev Valeriy Vasilievich (zaitzev.ssu.samara.ru), Karlov-juniour Alexander Vladimirovich, Telegin Sergey Sergeevich, Dept. of Radiophysics and Computer Modelling of Radiosystems, Samara State University, Samara, 443011, Russia.