CC BY
172
Вестник СамГУ — Естественнонаучная серия. 2009. № 6(72) 139
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 517.96: 574.34
ДВ-МОДЕЛЬ СИСТЕМЫ ’’ХИЩНИК-ЖЕРТВА”
© 2009 В.В. Зайцев, А.В. Карлов-младший, С.С. Телегин1
Принцип импульсной инвариантности теории синтеза цифровых фильтров использован для перехода к дискретному времени в модели Вольтерра с запаздыванием. Система ”хищник-жертва” представлена как дискретная автоколебательная система. Приведены результаты имитационного моделирования автоколебаний в системе со случайным запаздыванием.
Ключевые слова: модель Вольтерра, запаздывающие связи, автоколебания, дискретное время.
Введение
Математическое моделирование процессов в био- и экосистемах является одним из основных методов их теоретического исследования. Знаменитая модель Вольтерра [1], описывающая систему ”хищник-жертва”, известна уже более семидесяти лет. Тем не менее, исследования в рамках различных модификаций этой классической модели активно проводились и в последнее время [2]. Известны недостатки модели, главный из которых — ее консервативность, не позволяющие проводить количественных оценок характеристик процессов в системе ”хищник-жертва”. Согласно сложившимся в последние десятилетия представлениям о динамике системы ее основным движением являются устойчивые автоколебания численностей видов [3]. Следовательно, при формулировке динамической модели система должна рассматриваться как автоколебательная.
Одним из возможных способов преобразования модели Вольтерра к автоколебательной системе является учет естественных запаздываний в скоростях изменения видов по отношению к факторам, вызывающим данные изменения. При этом принимается во внимание способность запаздывающей
1Зайцев Валерий Васильевич ([email protected]), Карлов-младший Александр Владимирович, Телегин Сергей Сергеевич, кафедра радиофизики и компьютерного моделирования радиосистем Самарского государственного университета, 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.
