Мультистабильность для математической модели динамики хищников и жертв на неоднородном ареале

Полный текст

ОГЛАВЛЕНИЕ Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 509 1. Модель динамики популяций хищников и жертв . . . . . . . . . . . . . . . . . . . . . . 510 2. Косимметрия . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 510 3. Равновесия локальной системы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 512 4. Вычисление мультистабильности . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514 5. Разрушение семейства стационарных решений . . . . . . . . . . . . . . . . . . . . . . . 515 6. Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 517 Приложение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 517 Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 518 Введение В нелинейных эволюционных задачах математической физики нередко наблюдается мультистабильность - реализация различных устойчивых решений (аттракторов) в зависимости от начальных данных [14]. Помимо изолированных аттракторов, возможно возникновение континуальных семейств решений вследствие непрерывных симметрий [11, 18] или косимметрий [9]. Теория косимметрии, развитая для объяснения существования однопараметрического семейства стационарных режимов в задаче фильтрационной конвекции [9], позволила исследовать ряд задач популяционной динамики. В [2, 12, 15, 16] были вычислены семейства равновесий (стационарных распределений) и установлена индивидуальность спектра устойчивости, что является характерным отличием косимметрии от симметрии [11, 22]. В [6] для системы обыкновенных дифференциальных уравнений, описывающей динамику двух хищников и двух жертв на однородном ареале, было аналитически найдено двупараметрическое семейство равновесий и показана связь этой мультистабильности с мультикосимметрией задачи [3]. В работах [7, 8] был проведен анализ уравнений реакции-диффузии-адвекции для моделирования динамики хищников и жертв в случае неоднородного ареала. Было показано, что функция © Российский университет дружбы народов, 2022 This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4.0/legalcode 509 ресурса жертвы должна учитываться при формулировании соотношений функционального отклика, описывающих локальное взаимодействие хищника и жертвы. В настоящей работе рассматривается случай неоднородного ареала, на котором взаимодействуют два хищника и два родственных вида жертв. Для уравнений с модифицированным функциональным откликом Холлинга второго рода [17] найдены условия на параметры системы, при которых у нее имеется мультикосимметрия. Выведены формулы для равновесий, составляющих двупараметрическое семейство стационарных решений в бездиффузионном приближении. В численном эксперименте показано, что и в случае неоднородного ареала реализуется мультистабильность в виде двупараметрического семейства стационарных распределений с индивидуальным спектром устойчивости. На основе полученных решений рассмотрена инвазия второго хищника в экологическую систему, состоящую из хищника и двух родственных видов жертв. 1. Модель динамики популяций хищников и жертв Рассматривается модель динамики двух хищников и двух жертв в случае неоднородного ареала на основе системы уравнений [2, 12]. Распределение видов в момент времени t на одномерном ареале x ∈ [0, 1] дается функциями плотностей uj(x,t) для жертв (j = 1,2) и vj(x,t) для хищников (j = 1,2). Уравнения баланса видов записываются через миграционные потоки qj и функции локального взаимодействия fj: (1.1) , (1.2) где kj - коэффициенты диффузии, а направленная миграция (таксис) для жертв и хищников определяется формулами , (1.3) Здесь, αj, βij - миграционные коэффициенты, p(x) - функция ресурса. Для реакций fj используются представления на основе функционального отклика Холлинга второго рода [7, 8, 17]: ; (1.4) (1.5) Функция f определяет гиперболический закон роста жертв [5]. Локальное взаимодействие регулируется положительными коэффициентами μj (прирост жертвы), lij (убыль жертвы из-за хищника), lj (смертность хищника), μij (прирост хищника за счет жертв), а также cj. Рассматривается кольцевой ареал, система (1.1)-(1.3) дополняется условиями периодичности uj (0,t) = uj (1,t), qj (0,t) = qj (1,t), j = 1,2, vj (0,t) = vj (1,t), qj+2 (0,t) = qj+2 (1,t), j = 1,2, и начальными распределениями плотностей популяций (1.6) uj (x,0) = u0j (x), vj (x,0) = vj0 (x), j = 1,2. (1.7) 2. Косимметрия При дополнительных условиях на параметры система (1.1)-(1.7) относится к классу косимметричных динамических систем [9], для которых возможно возникновение непрерывных семейств решений. Согласно [9] для уравнения u˙ = F(u) косимметрия L представляет собой векторное поле, ортогональное F в каждой точке фазового пространства. Следующая лемма устанавливает существование косимметрии для рассматриваемой системы при выполнении дополнительных соотношений между регулирующими динамику жертв параметрами. Лемма 2.1. Косимметрией системы (1.1)-(1.7) является векторное поле - , (2.1) если выполнены следующие условия на параметры модели: , (2.2) где γ1 -константа. Доказательство. По определению косимметрии векторное поле L должно быть ортогонально правой части системы (1.1)-(1.7) для любых функций ui(x,t), vj(x,t), т. е. , где Подстановка (2.1) в I1 и перегруппировка слагаемых позволяет записать интеграл I1 следующим образом: Из условий (2.2) следует, что u1f2 = γ1u2f1, поэтому I1 = 0. После интегрирования по частям интеграл I2 представляется в виде суммы Внеинтегральное слагаемое пропадает в силу условий периодичности (1.6). Из условий (2.2) получаем и . Аналогично выводится . В результате Векторное поле L1 = (ζ1,ζ2,0,0) из (2.1)-(2.2) ортогонально правой части системы (1.1)-(1.7) и является косимметрией системы (1.1)-(1.7). В следующей лемме формулируются условия существования косимметрии в системе (1.1)-(1.7) при дополнительных соотношениях на параметры, регулирующие динамику хищников. Лемма 2.2. Косимметрией системы (1.1)-(1.7) является векторное поле - , (2.3) если выполнены условия на параметры модели: , (2.4) где γ2 -константа. Доказательство. Аналогично проведенному для леммы 2.1. Леммы 2.1 и 2.2 позволяют сформулировать теорему о мультикосимметрии системы (1.1)-(1.7). Теорема 2.1. При выполнении условий (2.4) и (2.2) система (1.1)-(1.7) имеет мультикосимметрию L = κL1 + (1 - κ)L2, κ ∈ [0, 1]. 3. Равновесия локальной системы В бездиффузионном приближении координата x является дополнительным параметром. В этом случае система (1.1)-(1.7) допускает однородные по x решения, которые находятся из системы обыкновенных дифференциальных уравнений , (3.1) (3.2) Для любых значений параметров система (3.2) имеет неустойчивое нулевое равновесие u1 = u2 = v1 = v2 = 0 и семейство равновесий для различных соотношений жертв (без хищников) . (3.3) Для равновесия из семейства (3.3) характеристическое уравнение имеет вид σ (σ - σ2)(σ - σ3)(σ - σ4) = 0, σ2 = p(x)(μ1θ - θμ2 - μ1), σ3 = -θ(c1l2 - c2l2 - μ21 + μ22) - c1l2 - l2 + μ21 , c1θ - c2θ - c1 - 1 σ4 = -θ(c1l1 - c2l1 - μ11 + μ12) - c1l1 - l1 + μ11 . c1θ - c2θ - c1 - 1 Нулевое значение отвечает нейтральному направлению вдоль семейства. Видно, что σ2 < 0 для всех θ ∈ [0, 1] и любых x, семейство (3.3) устойчиво при условиях θ(c1lj - c2lj - μj1 + μj2) - c1lj - lj + μj1 < 0, j = 1,2. (3.4) При наличии дополнительных условий на параметры система может иметь равновесия, отвечающие сосуществованию одного хищника и одной жертвы: ; (3.5) ; (3.6) ; (3.7) Равновесия, отвечающие сосуществованию всех видов, получаются из системы алгебраических уравнений (3.9) . (3.10) Аналогично лемме 2.1 получается, что при выполнении условий на параметры (3.11) система (3.1)-(3.2) имеет косимметрию L1 = (γ1u2,-u1,0,0)T . При условиях система (3.1)-(3.2) имеет косимметрию T (3.12) L2 = (0,0,γ2v2,-v1) . При совместном выполнении условий (3.11) и (3.12) векторное поле задачи (3.1) и (3.2) ортогонально линейной комбинации L1 и L2: L ≡ (1 - κ)L1 + κL2, κ ∈ [0, 1]. (3.13) Это соответствует мультикосимметрии. С учетом (3.11) и (3.12) уравнения (3.9)-(3.10) примут вид . Двупараметрическое семейство равновесий получается в явной форме: , (3.14) где (3.15) Для положительности решений параметры должны удовлетворять условиям μ11 > c1l1, μ12 > c2l1, 0 < ξ < l1, 0 < ψ < ψ∗. Характеристическое уравнение для равновесия с континуальным номером (ξ,ψ) имеет вид , Лемма 3.1. Устойчивость равновесий двупараметрического семейства (3.14) не зависит от пространственной координаты x. Доказательство. Из (3.14) следует, что выражения для жертв и хищников соответственно пропорциональны p(x) и p2(x). Отсюда получается, что коэффициент A есть произведение p(x) и выражения, не содержащего x. Аналогично, коэффициент B обратно пропорционален p(x). Следовательно, устойчивость равновесий семейства зависит от параметров системы и не зависит от координаты x. В [6] проанализирована устойчивость равновесий двупараметрического семейства при однородной функции ресурса p(x) = 1 и cj = 0. В случае c1 = c2 = c, γ1 < 1 и при μ11 > μ12 > (c+1)l1 получается следующее условие устойчивости равновесия из (3.14): . (3.16) 4. Вычисление мультистабильности В случае ненулевых коэффициентов диффузии и таксиса задача (1.1)-(1.7) решается численно. Используется метод прямых и схема смещенных сеток, см. приложение и [12]. Для функции ресурса p(x) = 1 + 0,4sin2πx и фиксированных значений коэффициентов μ1 = l13 = l14 = 1, γ1 = 0,5, μ11 = 1,5, μ12 = 0,875, l1 = 0,5, γ2 = 1,5, c = 0,2 в бездиффузионном приближении получается двупараметрическое семейство равновесий , где . Здесь координата x играет роль дополнительного параметра. Семейство равновесий устойчиво при . В численном эксперименте показано, что данное семейство трансформируется при ненулевых коэффициентах диффузии, которые удовлетворяют соотношениям (3.11)-(3.12). При этом происходит сглаживание распределений видов по координате x. На рис. 1 приведены результаты установления при k1 = 0,005, k3 = 0,0025 из начальных данных, отвечающих равновесию в бездиффузионном приближении. Рис. 1. Установление стационарного распределения жертв и хищников из равновесия W(0,4,0,15) (4), k1 = 0,005, k3 = 0,0025. Fig. 1. Establishment of a stationary distribution of prey and predators from the equilibrium W(0,4,0,15) (4), k1 = 0,005, k3 = 0,0025. Далее для полученных стационарных распределений численно определялся спектр устойчивости, см. табл. 1. Видно, что первые два значения практически нулевые (10-11, 10-9). Это соответствует двум нейтральным направлением для семейства стационарных распределений. Следующие отрицательные значения порядка единицы характеризуют устойчивость в трансверсальных к семейству направлениях. Эти числа различны для решений, что подтверждает косимметричный характер полученного семейства. В эксперименте в качестве начальных данных брались равновесия при ξ = 0,35, ψ = 0,1 (I), ξ = 0,4, ψ = 0,1 (II), ξ = 0,4, ψ = 0,15 (III). Получившиеся стационарные профили распределений жертв (верх) и хищников (низ) при k1 = 0,005, k3 = 0,0025 представлены на рис. 2. Видно, что плотности хищников и жертв соответствуют характеру изменения функция ресурса p(x), но возможны различные комбинации видов. Имеются также неустойчивые стационарные решения, от которых ответвляются колебательные режимы (при 0 < ξ < 0,288). Рис. 2. Три стационарных распределения из двупараметрического семейства. Fig. 2. Three stationary distributions from a two-parameter family. σ1 σ2 σ3 σ4 I II III Таблица 1. Спектр стационарных распределений видов, приведенных на рис. 2. Table 1. The spectrum of stationary distributions of the species shown in Fig. 2. 5. Разрушение семейства стационарных решений Возникновение в системе двупараметрического семейства равновесий является следствием дополнительных условий, которым должны удовлетворять параметры системы. При нарушении некоторых из этих соотношений теряется косимметрия и происходит разрушение семейства. В [10] развит метод исследования нарушения косимметрии, основанный на вычислении косимметричного дефекта и анализе селективной функции. Этот метод позволяет указать на реализующиеся в результате разрушения семейства решения. В [16] для задачи о кинетике трех видов с учетом диффузии и линейной адвекции численно были найдены случаи трансформации семейства в предельный цикл и распада семейства на изолированные равновесия. В [15] показано, что при помощи метода селективной функции можно анализировать задачи инвазии - внедрения нового вида в экологическую систему. Далее рассматривается задача для системы (1.1)-(1.7) и считается, что вид v2 является инвазивным. Для локальной системы (3.1)-(3.2) при выполнении условий (3.11)-(3.12) имеется мультикосимметрия, и инвазия является успешной для ненулевых начальных распределений v2. В результате установления формируется решение, соответствующее сосуществованию двух хищников и двух жертв. Систему (3.1)-(3.2) можно переписать, предусмотрев возмущение параметров μ21, μ22, которые вместе с l14, l24 участвуют в описании взаимодействия хищника v2 с жертвами. В результате получается система с добавившимися параметрами ω1, ω2, изменения касаются двух уравнений из четырех (5.1) (5.2) (5.3) . (5.4) Вычислим косимметрический дефект системы (5.1)-(5.4), следуя [10]: (5.5) Здесь . Селективная функция S получается подстановкой формул для равновесий семейства (3.14) в косимметрический дефект D. Поскольку знаменатель в D (см. (5.5)) положителен, то его можно опустить, имея в виду анализ нулей селективной функции. В результате получается следующая функция трех параметров где . Нулями селективной функции, принадлежащими области определения, являются ψ = 0, (5.6) (5.7) . (5.8) - - Случай ψ = 0 отвечает решению исходной системы при v2 = 0, это означает неуспех инвазии. Из (5.7) получается решение с одним хищником v2, т. е. новый хищник «вытеснил» старого (v1 = 0). Для 0 < ψ < ψ∗ из (5.8) получается ξ = ξ∗, , (5.9) и однопараметрическое семейство равновесий W(ξ∗,ψ). k1 = 2k3 σ1 σ2 σ3 σ4 0,001 0,005 0,01 - × - × - ± Таблица 2. Спектр стационарных распределений видов W(ξ∗,ψ), ω1 = 0,08, ω2 = 0,01. Table 2. Spectrum of stationary distributions of the species W(ξ∗,ψ), ω1 = 0,08, ω2 = 0,01. k1 = 2k3 σ1 σ2 σ3 σ4 0,001 0,005 0,01 - × - × - ± Таблица 3. Спектр стационарных распределений видов W(ξ∗,ψ), ω1 = 0,45, ω2 = 0,05. Table 3. Spectrum of stationary distributions of the species W(ξ∗,ψ), ω1 = 0,45, ω2 = 0,05. В таблицах 2 и 3 представлены результаты вычисления спектра устойчивости для случаев ω1 = 0,08, ω2 = 0,01 и ω1 = 0,4, ω2 = 0,05 (ξ∗ = 0,408), соответственно. Для установления стационарных распределений расчеты проводились на промежутке времени t ∈ [0, 1000]. По сравнению с табл. 1, в каждой строке табл. 2 имеется только одно практически нулевое значение 10-9 (нейтральный спектр). Это является следствием нарушения мультикосимметрии, в результате которого остается косимметрия L1 и однопараметрическое семейство решений. В эксперименте обнаружено, что установление при больших значениях ωj требует большего временного промежутка. Наличие в спектре устойчивости величин порядка 10-4 может быть связано с «памятью» системы об исчезнувших равновесиях. 6. Заключение Математические модели популяционной динамики являются необходимым инструментом при прогнозировании экологических и биологических процессов [4, 20]. Для описания пространственно-временных миграций видов применяются различные варианты уравнений реакции-диффузии-адвекции [13]. Возникающие при этом задачи содержат большое количество экспериментально определяемых и при необходимости назначаемых параметров [7, 21]. Подход на основе теории косимметрии [3, 9, 12, 15, 22] позволяет найти соотношения на параметры, при которых имеется мультистабильность, а затем анализировать возмущения, нарушающие косимметрию. Возникающая при этом динамика напоминает череду состояний исчезнувшего семейства. Представленный в настоящей работе пример двупараметрического семейства стационарных распределений является основой для последующего анализа динамики системы при различных диапазонах значений параметров. Приложение Для дискретизации системы (1.1)-(1.7) по переменной x вводится равномерная сетка: . (6.1) Плотность распределения популяций uj, vj в узле xr далее обозначается через uj,r, vj,r. При вычислении потоков используется вспомогательная сетка: Для аппроксимации системы уравнений (1.1)-(1.5) по пространственной координате применяется метод баланса: уравнение (1.1) интегрируется по отрезку , а для потоков qj интегрирование проводится по отрезку [xr,xr+1]. Далее используются операторы разностной производной и вычисления среднего , (6.2) а также условия периодичности . (6.3) В результате получается система обыкновенных дифференциальных уравнений для uj,r, vj,r. Из (1.1), (1.2) следует u˙j,r = [-dqj + fj]r , j = 1,2, r = 1,...,N, (6.4) v˙j,r = [-dqj+2 + fj+2]r , j = 1,2, r = 1,...,N. Для локальных членов получается (6.5) (6.6) , (6.7) (6.8) где дискретный аналог функции ресурса определяется следующей формулой: (6.9) Выражения для потоков жертв получаются из (1.1) и (1.4): , (6.10) а из (1.2) и (1.4) находятся потоки для хищников: . (6.11) Построенная конечномерная модель (6.4)-(6.11) дополняется дискретными аналогами условий периодичности , (6.12) и может быть записана в виде W˙ = Φ(W), W (0) = W0. Здесь W - вектор значений переменных в узлах сетки: W = (u1,1,...,u1,N,...,u2,N,v1,1,...,v1,N,...,v2,N). Начальные данные для системы (6.3)-(6.12) получаются из (1.7): (6.13) . (6.14) Для интегрирования системы (6.13) по времени используется метод Рунге-Кутты 4-го порядка.

Об авторах

Т. Д. Ха

Южный федеральный университет; Вьетнамско-Венгерский индустриальный университет

Email: toanhd.viu@gmail.com
Ростов-на-Дону, Россия; Ханой, Вьетнам

В. Г. Цибулин

Южный федеральный университет

Автор, ответственный за переписку.
Email: vgcibulin@sfedu.ru
Ростов-на-Дону, Россия

Список литературы