«В.Н. Разжевайкин Модели динамики популяций ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР ИМ. А.А. ДОРОДНИЦЫНА РАН Москва 2006 УДК 519.6 Ответственный редактор доктор физ.-матем. наук А.П. Абрамов В работе излагается ряд методов исследования ...»

РОССИЙСКАЯ АКАДЕМИЯ НАУК

ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР

ИМ. А.А. ДОРОДНИЦЫНА

В.Н. Разжевайкин

Модели динамики популяций

ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР

ИМ. А.А. ДОРОДНИЦЫНА РАН

Москва 2006

УДК 519.6

Ответственный редактор

доктор физ.-матем. наук А.П. Абрамов

В работе излагается ряд методов исследования моделей

динамики биологических популяций и сообществ, отражающий современное состояние ряда направлений в математической экологии Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (коды проектов 07-07-00104, 06-01-00244).

Рецензенты: Д.О. Логофет, Д.А. Саранча Научное издание c Вычислительный центр им. А.А.Дородницына Российской академии наук, 2006.

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

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

Учитывая огромный поток литературы по рассматриваемой здесь тематике, а также глубокую проработанность многих из затрагиваемых вопросов, я позволил себе распределить внимание не везде достаточно равномерно, отдавая предпочтение тем из них, которые либо представляются мне наиболее существенными с точки зрения принципиальных подходов к математическому моделированию в биологии, либо не находят должным образом адаптированного под прикладные задачи подхода к выбору моделей. Поэтому, например, первый раздел, посвященный классификации моделей, изложен достаточно конспективно, в то время как второй, затрагивающий сравнительно простые вопросы одномерных систем, оказывается довольно большим по объему. Но такое внимание к этим вопросам оказалось неслучайным. Многие современные монографии, например, пестрят бифуркационными диаграммами типа «вилки» как представляющими бифуркацию потери устойчивости положения равновесия, хотя, как это хорошо известно и как это изложено в предлагаемой читателю книге, такая картина не является грубой, т.е. требует выполнения ряда специальных условий типа равенства или симметрии. Подобные условия зачастую возникают в качестве артефакта при попытках прикладного использования излишне простых математических моделей, таких, как, например, моделей возникновения диссипативных структур в системах уравнений реакции – диффузии на одномерных пространственных областях. Надеюсь, что расставленные мной акценты позволят предостеречь читателя, дабы он не принял эти артефакты за имеющие какое-либо отношение к реальным свойствам моделируемых объектов.

Много внимания уделено в работе также и бифуркации Хопфа. Я позволил себе не ограничиваться иллюстративным изложением и передать материал таким образом, чтобы у читателя осталось ощущение понимания описываемых явлений. Поскольку такое ощущение может возникнуть исключительно в процессе самостоятельной работы, мною приводится расчетная формула первой ляпуновской величины для плоских задач и некоторые результаты ее применения. Отмечу здесь, что эту, как и еще более громоздкие формулы, удобнее в настоящее время использовать с помощью компьютера, поскольку имеется большое количество достаточно глубоко проработанных математических пакетов (таких, как, например, последние версии Maple), позволяющих легко выполнять алгебраические преобразования.

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

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

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

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

У меня нет сомнений, что предлагаемая вниманию читателя работа может быть во многих отношениях улучшена.

Поэтому прошу присылать ваши замечания и пожелания по адресу: [email protected] В.Н. Разжевайкин. Москва, декабрь 2006 г.

1. Структура математических Таксономическая структура математических моделей биологических систем определяется следующими дихотомиями.

1) Имитационные – аналитические модели Имитационные модели предназначены для решения задач прогнозирования динамики реальных биологических систем и управления ими.

Объектами моделирования являются:

– лабораторные популяции, – культивируемые сообщества, – эксплуатируемые популяции, – экосистемы, подверженные антропогенным нагрузкам, – население регионов, – эпидемии, эпизоотии, эпифитотии.

Методика исследования имитационных моделей основывается на следующей последовательной схеме:

— концептуальная модель, – математическая модель, – идентификация параметров модели, – численные эксперименты и анализ модели, – верификация модели, – наработка сценариев управления и прогнозов по ним.

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

Примеры имитационных моделей.

– Моделирование роста хлопчатника при наличии паразитов.

– Моделирование движения рыбы в Великих озерах Северной Америки.

– Моделирование распространения сине-зеленых водорослей при эвтрофикации пресноводных водоемов.

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

Аналитические модели предназначены для изучения причин наблюдаемых в природных системах характерных явлений. Объектами моделирования являются:

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

Методика исследования аналитических моделей основывается на сокращенной схеме:

– концептуальная модель, – математическая модель, – аналитический и численный анализ модели, – выводы.

Замечание. На стадии анализа допускается корректировка как математической, так и концептуальной моделей.

Принципиальные отличия от имитационных моделей следующие:

1. Отсутствует необходимость привязки к измеряемым параметрам.

2. На выходе получаются результаты качественного характера.

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

Примеры математических моделей и описываемых ими явлений.

– Системы «хищник – жертва» (колебательные режимы).

– Конкурентный отбор (занятие экологических ниш).

– Пандемии (порог эпидемии).

– Системы «реакция – диффузия» (скорость пространственного распространения, пространственные структуры).

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

2) Непрерывные – дискретные модели – Непрерывные модели имеют независимую переменную времени t R. При отсутствии других независимых переменных описываются системами обыкновенных дифференциальных уравнений. Наличие других независимых переменных (как правило, тоже непрерывных) приводит к уравнениям более сложного вида (уравнения с частными производными, уравнения с отклоняющимся аргументом, интегродифференциальные уравнения).

– Дискретные модели имеют независимую переменную времени t Z. Описываются разностными уравнениями.

Предпочтительны при наличии жесткой временной периодизации (сезонность, непересекающиеся поколения и т.п.) 3) Линейные – нелинейные модели – Линейные модели используются для описания стадий свободного (при отсутствии лимитирующих факторов) развития и в качестве аппроксимации нелинейных (например, в окрестности положения равновесия).

– Нелинейные модели используются при учете внутривидовых и межвидовых взаимодействий популяций.

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

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

Помимо перечисленных дихотомий математические модели различаются таксономически также следующими признаками.

5) Выбираемые фазовые переменные В качестве таковых принято выбирать:

1. Численности популяций (популяционная динамика),.

2. Составляющие вещества и энергию.

3. Частоты генов (математическая генетика).

6) Факторы структурированности По этому признаку модели различаются на следующие:

1. Точечные (отсутствие какой-либо структуры).

2. Модели с пространственной структурой.

3. Модели с возрастной структурой.

4. Модели с отклоняющимся аргументом.

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

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

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

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

2. Численность популяции достаточно велика для того, чтобы использовать детерминированные модели.

Простейшей из таких моделей является 1) Модель Мальтуса Пусть N(t) – численность популяции в момент времени t, b – ее рождаемость (темп воспроизводства в расчете на одну особь), d – смертность (темп вымирания), a = b d – мальтузианский параметр. Прирост популяции N за время t определяется как разница между числом родившихся и умерших особей: N = bNt dNt, так что динамика численности популяции задается решениями задачи Коши имеющими вид экспоненты N(t) = N0 eat.

2) Логистическая модель Включение внутривидовой конкуренции в виде линейно убывающего по численности коэффициента прироста приводит к уравнению вида Здесь a > 0, K > 0 – емкость среды, определяющая предельное значение для численности. Уравнение (2.2) интегрируется разделением переменных, и его решение с начальными условиями, как в (2.1), имеет вид Заметим, что для N = K N уравнение (2.2) может быть переписано в виде что отличается от (2.2) направлением времени. Отсюда следует:

1. Симметрия кривой N(t) относительно точки (t = 0, N = K/2) 2. Асимптотика является экспоненциальной по обоим направлениям времени.

Замечания.

1. Вид решения (2.3) не зависит от знака величин a и N0 K.

2. Вид решения (2.3) сохраняется также и в неавтономном случае a = a(t) Надо лишь заменить везде произведение at интегралом a(s)ds В случае расходимости этого интеграла по какому-либо направлению сохраняется и асимптотика по этому направлению.

3) Модель поиска партнера При возникновении трудностей в поисках брачного партнера модель может принимать вид В окрестности нулевого положения равновесия уравнение (2.4) аппроксимируется уравнением dN решения в этой окрестности N(t).

Замечания.

1. Малая скорость размножения при малых численностях может повлечь исчезновение вида при наличии дополнительных факторов (например, пространственной распределенности).

2. При a < 0 этот же эффект имеет обратный характер:

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

4) Модель поиска партнера при наличии смертности Учет естественной смертности в уравнении (2.4) приводит к уравнению где N1,2 = K 1 ± 1 dc, а dc = aK – критическое значеd ние для коэффициента смертности d, превышение которого приводит к потере нетривиальных положений равновесия.

5) Модель безымунной эпидемии Пусть N – постоянная численность популяции носителя, I – число инфицированных, a – коэффициент заражения, и d – коэффициент выздоровления. Исходное уравнение для доли инфицированных p = может быть переписано в виде говое значение численности для возникновения эпидемии).

При этом N может служить параметром, по которому при N = Nc наблюдается бифуркация смены устойчивой ветви (см. ниже). При прохождении этого значения тривиальное равновесие теряет устойчивость, которая передается нетривиальному положительному равновесному состоянию p = k.

Замечание. Учет естественной смертности в уравнении (2.2) также дает уравнение типа (2.6), т.е. снова сводится к уравнению вида (2.2).

3. Бифуркации положения 1) Грубый случай одномерной однопараметрической бифуркации Под грубостью системы (синоним – система общего положения) понимается ее устойчивость в смысле схожести фазовых портретов (наборов направленных фазовых траекторий) для нее и для возмущенной системы по отношению к малым возмущениям в пространстве непрерывно дифференцируемых функций, применяемым к правой части описывающих ее дифференциальных уравнений (см., например, [2]).

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

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

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

При последовательном увеличении числа параметров возникают все новые неустранимые шевелениями вырождения.

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

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

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

x R с параметром R. Нарушение неравенства требует выполнения двух условий типа равенства, что в грубом случае соответствует изолированным точкам на плоскости (x, ). При этом можно предполагать, что в этих точках правая часть (3.1) отлична от нуля, так что на множестве ее корней выполнено (3.2). Но последнее условие означает локальную применимость теоремы о неявной функции, по крайней мере, относительно одной из переменных на указанной плоскости, из чего в свою очередь вытекает возможность представления множества корней в виде одномерной кривой на этой плоскости. На невырожденных участках кривой, где выполнено неравенство fx = 0, эта кривая может быть задана некоторой гладкой функцией корней x = g(), причем устойчивым корням соответствует случай fx < 0, а неустойчивым – fx > 0. Компоненты кривой, задаваемой условием f (x, ) = 0, разделяют плоскость (x, ) на две области с противоположными (вертикальным и горизонтальным) направлениями фазовых траекторий, согласующимися со знаками указанных выше производных. Это позволяет чисто геометрически находить знак производной fx в каждой точке кривой только на основании информации об этом знаке для единственной из ее невырожденных точек. В тех точках кривой (xc, c ), где выполнено равенство fx = (точки складки), она может быть параметризована в виде:

= c +(xxc )2 +o ((x xc )2 ). (Из разложения f (x, ) в ряд Тейлора в точке (xc, c ) при f (xc, c ) = fx (xc, c ) = 0 получается указанное выше разложение с остаточным членом вида o ((x xc )2 + | c |). Последующее применение теоремы о неявной функции по отношению к переменной позволяет избавиться от последнего слагаемого в остаточном члене.) Случаю > 0 соответствует рождение пары стационарных решений уравнения (3.1) (устойчивого и неустойчивого) при увеличении параметра в точке = c. В случае противоположного знака наблюдается слияние таких корней и их исчезновение.

Те точки плоскости (x, ), где нарушено условие (3.2), являются в случае общего положения либо изолированными точками локального экстремума функции f (x, ), либо ее седловыми точками. Если такая точка оказывается еще и корнем этой функции, то малое шевеление вида f (x, ) + позволяет в первом случае либо избавиться от нее, либо получить вместо нее замкнутую кривую, близкую к эллипсу, а во втором – вместо двух пересекающихся под ненулевым углом кривых получить две непересекающиеся, но сходящиеся близко друг к другу в точках своей максимальной кривизны кривые, близкие в окрестности этих точек к гиперболам.

2) Ветвления в случае одномерной однопараметрической бифуркации В случае, когда правая часть уравнения (3.1) имеет вид и шевеления допускаются только в отношении каждого из сомножителей в (3.3), к описанным выше бифуркациям множества корней уравнения (3.1) добавляются ветвления, соответствующие точкам пересечения кривых i, задающих корни сомножителей fi (x, ). Поскольку при этом то на ветви i = {x : fi (x, ) = 0} все слагаемые в (3.4), за исключением i-го обратятся в нуль, так что для x i будет выполнено fx (x, ) = f1 (x, )... fi1 (x, )fix (x, )fi+1 (x, )... fn (x, ).

Поскольку при i = j множитель fj (x, ) присутствует в (3.5), то отсюда следует, что в точках пересечения i и j fx (x, ) = 0. То же верно, очевидно, и для производной по, так что во всех точках пересечения различных ветвей неравенство (3.2) нарушается. Как и в предыдущем случае, вне ветвей знак производной сохраняется и меняется на противоположный при пересечении каждой из них. Поскольку точки тройных пересечений, пересечения с касанием и «вилки»

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

Для построения полной картины устойчивости опять достаточно знать знак fx (x, ) лишь для одной нестационарной точки. Неустранимость такой бифуркации следует из сохранения (возможно с малым смещением) ветвей i при шевелении определяющих их сомножителей fi (x, ).

Заметим, что описанная здесь картина бифуркации пересечения различных ветвей корней уравнения (3.1) соответствует случаю общего положения в ситуации нарушения неравенства (3.2). Действительно, пусть для определенности имеется вырождение в начале координат, так что f (0, 0) = fx (0, 0) = f (0, 0) = 0. В этом случае при разложении в ряд Тейлора старшими окажутся квадратичные члены: f (x, ) = ax2 + 2bx + c2 + o(x2 + 2 ). В случае общего положения квадратичная форма невырождена, так что Эта форма посредством линейного ортогонального преобразования (x, ) (u, v) приводится к диагональному виду, так От изолированной особой точки, соответствующей случаю > 0, можно избавиться шевелением сомножителей функции f (x, ) точно так же, как это было описано в предыдущем параграфе. В случае < 0 множество корней функции f (u, v) представляет собой две гладкие кривые, пересекающиеся в начале координат. Действительно, переписывая разложение в ряд Тейлора в виде f (u, v) = u2( + U(u, v)) + v ( + V (u, v)), где U(0, 0) = V (0, 0) = 0, и полагая для определенности > 0, находим, что f (u, v) = g+ (u, v)g(u, v), где функции g± (u, v) = u + U(u, v) ± v V (u, v) удовлетворяют в окрестности начала координат условиям теоремы о неявной функции, обеспечивающей их гладкое некасательное пересечение.

Замечание. Если отказаться от условия ортогональности линейного преобразования формы, то можно считать || = || = 1.

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

Точнее (см., например, [2]), пусть f : Rn Rk Rn – достаточно гладкое отображение, и f (0, 0) = 0, так что начало координат является положением равновесия системы при значениях параметров = 0. Если n±, n0 – число собственных значений якобиана x f (x, 0) при x = 0, расположенных соответственно справа, слева и на мнимой оси, то существует такая непрерывно зависящая от параметров непрерывная замена фазовых переменных x (, y, z), где Rn0, y Rn, z Rn+, оставляющая на месте начало координат, что семейство систем (3.6) в его окрестности является надстройкой седла над семейством с фазовым пространством размерности n0, называемым нейтральным многообразием, т.е. в новых переменных решения системы (3.6) задаются системой уравнений вида:

Локальные инвариантные многообразия, образующие указанную надстройку, называются устойчивым (многообразие размерности n в окрестности положения равновесия x = системы (3.6), переходящее при указанном преобразовании в компоненту y системы (3.7)) и, соответственно, неустойчивым (размерности n+, переходит в компоненту z). При изменении знака времени эти многообразия меняются местами. Особый с прикладной точки зрения интерес представляет случай, когда n+ = 0, соответствующий устойчивости по компонентам настройки. В этом случае устойчивость положения равновесия полностью определяется его устойчивостью по отношению к возмущениям из нейтрального многообразия. В частности, все перестройки, связанные с изменением характера устойчивости при изменениях параметров, полностью определяются перестройками на этом многообразии, что является основной причиной, определяющей важность исследования бифуркаций для случаев малой размерности фазового пространства.

Заметим, что в этом случае надстройка уже не является седлом, а всего лишь узлом, но это не меняет терминологии. При этом следует отметить, что все устойчивые узлы и фокусы одинаковой размерности эквивалентны между собой в смысле формулировки приведенной выше теоремы (см.

[1]). Такая эквивалентность называется топологической эквивалентностью. В отличие от топологической орбитальной эквивалентности она при отображении траекторий не только оставляет на месте направление изменения времени, но и его значение. Это означает в частности, что она является более сильной в том смысле, что соответствующие ей классы эквивалентности являются более узкими. Так, например, две орбиты (фазовые траектории периодических решений), имеющие различные значения периодов, не являются топологически эквивалентными, хотя и становятся таковыми при ослаблении эквивалентности до орбитальной. Именно это обстоятельство и послужило главной причиной невозможности использования топологической эквивалентности при определении понятия грубости системы.

4) Пример системы «ресурс – потребитель»

с независимым ресурсом Пусть R(t) – количество ресурса в момент времени t, динамика которого описывается уравнением где функция f такая же, как и в одном из разобранных выше примеров. Пусть эффективная (т.е. учитывающая смертность) рождаемость потребителя x зависит только от величины доступного для него ресурса D = Rx ( – коэффициент использования ресурса потребителем), так что где V (D) > 0, V (0) = d < 0, V () = VM (0, ).

Пусть Ds > 0 – корень функции V (D), Rs = Rs () – корень функции f (R, ). Тогда, если Rs > Ds, то построенная система уравнений имеет положительное равновесное состояние (Rs, xs ), где xs = Rs Ds, устойчивость которого определяется исключительно устойчивостью его первой компоненты, т.е. знаком fR (Rs, ) (ибо якобиан является треугольным и его ненулевой диагональный член, соответствующий x, равен xs V (Ds ) < 0).

При слиянии двух корней первого уравнения (как, например, в системе 4) второго раздела) в условиях выполнения неравенства Rs () > Ds происходит бифуркация типа «седло – узел», являющаяся двумерной надстройкой над одномерной бифуркацией исчезновения (или появления) двух корней.

5) Одномерные бифуркации с большим В контексте теоремы сведения размерность n0, на которой происходят бифуркации положения равновесия, определяется числом собственных значений якобиана системы (3.6), расположенных на мнимой оси. Одномерность означает наличие единственного такого собственного значения в нуле, что в случае общего положения при наличии единственного параметра (k = 1) соответствует наличию неустранимой особенности типа «складка», соответствующей бифуркации рождения или исчезновения пары положений равновесия:

устойчивого и неустойчивого. Увеличение размерности пространства параметров приводит к появлению новых неустранимых особенностей. Так, оставаясь в рамках одной фазовой переменной, можно наблюдать появление особенностей типа «сборка» для k = 2 и типа «ласточкин хвост» для k = 3.

При k > 1 в случае общего положения неустранимым образом могут обнаруживаться изолированные точки(x, ) R Rk, в которых так что если считать такую точку началом координат, то правая часть уравнения (3.1), в котором уже = (1, · · ·, k ), приобретает при = 0 вид f (x, 0) = xk+1 (a + g(x)), где a = 0, g(0) = 0. Не меняющей знака заменой фазовой переменной y (x) = x k+1 |a + g(x)|, можно добиться того, чтобы в некоторой окрестности точки x = 0 функция f(y) такая, что f ((x)) = f (x, 0), имела вид f(y) = ±y k+1. Здесь знак совпадает со знаком signa. То, что левая часть в (3.1) при такой замене при = 0 изменится, не принципиально, поскольку при этом сохраняется ее знак (т.е. направление траектории) и неподвижная точка x = y = 0, что соответствует топологической эквивалентности уравнения (3.1) при = 0 и уравнения Заметим, что при нечетном k от знака «минус» в (3.9) можно избавиться, изменив знак y, поскольку в этом случае происходит только изменение знака левой части (3.9).

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

Это означает, в частности, что, вычислив все слабо неустранимые особенности заданной коразмерности, мы тем самым не пропустим и сильно неустранимые особенности той же коразмерности. В качестве узкого класса возмущений для стационарной точки y = 0 уравнения (3.9) принято выбирать возмущения в виде канонической формы Уитни (см., например, [3]) Pk (x, ) = 1 + 2 x + · · · + k xk1, так что уравнение (3.9) вкладывается при 1 = 2 = · · · = k = 0 в k – параметрическое семейство:

Оказывается, однако, что все особенности коразмерности k уравнения (3.9) не только являются подмножеством множества особенностей уравнения (3.10), но и в точности совпадают с ним в смысле топологической эквивалентности фазовых пространств, расширенных на пространство параметров одинаковых размерностей. При этом в процессе построения замены переменных, устанавливающей эту эквивалентность, удается посредством сдвига, определяющего выделение полной (k + 1)-й степени, избавиться от члена степени k, что, собственно, и обуславливает приведенный выше вид возмущения в виде канонической формы Уитни.

Заметим также, что если не интересоваться знаком устойчивости положений равновесия семейства уравнений (3.10), Рис. 1. Трехмерное изображение поверхности то для него можно всегда выбирать знак «плюс» независимо от четности k, поменяв в случае необходимости направление времени. В этом случае обычно в (3.10) ограничиваются поиском положений равновесия, т.е. считают левую часть равной нулю:

6) Бифуркация типа «сборки»

При k = 1 каноническая форма Уитни (3.11) задает уже знакомую нам «складку» W1 (, y) = y 2 + = 0, на которой происходит рождение (исчезновение) пары стационарных решений y = ±. При k = 2 получаем вместо (3.11):

При 1 = 2 = y = 0 для (3.12) наблюдается особенность типа «сборка». Геометрически она действительно напоминает сборку двумерной поверхности 2, задаваемой уравнением (3.12) в пространстве трех переменных (1, 2, y). Она устроена следующим образом (см. рис.1). В сечении 2 плоскостью 2 = 0 получается кубическая парабола 1 + y 3 = 0. При 2 = 0 она превращается либо в монотонную кривую (для 2 > 0), либо в кривую с тремя участками монотонности (для 2 < 0). В последнем случае точки изменения направления монотонности, т.е. экстремальные точки кривой – это точки складки, определяемые, как обычно, из условия равенства нулю первой производной по фазовой переменной:

Одновременное выполнение условий (3.12) и (3.13) задает трехмерную кривую S2 складок, описываемую параметрически как (1 = 2y 3, 2 = 3y 2 ). В проекции S2 на плоскость (1, 2 ) получается дискриминантная кривая 2, имеющая вид полукубической параболы 2 : 2 0, D(1, 2 ) = 2 2 = 0. Ее острие, т.е. точка возврата, получающаяся при 2 = 0, называется точкой сборки, а обе линии 2 – линиями складок. Внутренняя (меньшая) область плоскости параметров, ограниченная 2, т.е. определяемая двумя соотношениями 2 < 0, D(1, 2 ) > 0, соответствует наличию трех корней уравнения (3.12), а внешняя, получаемая при 2 > 0, или же при D(1, 2 ) < 0, – только одного корня. При пересечении 2 в случае непрерывного изменения параметров происходит обычное исчезновение (возникновение) двух таких корней. Внутренняя часть поверхности 2, ограниченная трехмерной кривой складок, состоит из устойчивых положений равновесия уравнения dy = 1 + 2 y + y 3, в то время как остальная ее часть – из неустойчивых. Изменение направления времени (т.е. знака левой части уравнения) меняет характер устойчивости на противоположный.

7) Бифуркация типа «ласточкин хвост»

При k = 3 уравнение (3.11) принимает вид:

Рис. 2. Трехмерное изображение поверхности 3.

для которого при 1 = 2 = 3 = y = 0 возникает особенность типа «ласточкин хвост». В качестве бифуркационной диаграммы этой особенности обычно рассматривается разбиение трехмерного пространства параметров (1, 2, 3 ) на области посредством двумерной дискриминантной поверхности 3 (см. рис. 2), являющейся проекцией поверхности складок S3 на это пространство. Поверхность складок S задается в пространстве (1, 2, 3, y) системой уравнений сительно (1, 2 ), получаем параметрическое представление для ее проекции: 3 : 2 = 4y 3 23 y, 1 = 3y 4 + 3 y 2.

Точки пространства параметров, не попадающие на саму поверхность 3, принадлежат одной из трех областей (см.

рис. 3), различающихся между собой количеством вещественных корней уравнения (3.14): два корня для самой большой Рис. 3. Сечение поверхности 3 при 3 = из них (куда попадают, например, точки = (1, 2 3 ) с 1 < 0, 3 > 0), отсутствие корней для другой неограниченной на сечениях 3 = const (сюда попадают, например, точки с 1 > 0, 3 > 0), и четыре корня для области, лежащей в полупространстве 3 < 0 и имеющей на указанных сечениях вид ограниченной треугольной области (к ней можно отнести, например, точки с 2 = 0, 1 > 0, 3 < 0). При этом в каждой из этих областей устойчивость корней по отношению к уравнению dy = 1 +2 y +3 y 2 +y 4 чередуется в порядке их возрастания таким образом, что наименьший из них является устойчивым. На самой поверхности 3 наблюдается бифуркация складки. На линии ее самопересечения эта бифуркация двойная, а на ее линиях возврата – бифуркация сборки. Эти линии составляют кривую сборок, задаваемую соотношениями W3 (, y) = W3 (,y) = W3 (,y) = 0, имеющими параметриy ческое решение (1 = 3y, 2 = 8y 3, 3 = 6y 2).

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

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

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

2) Теорема А.Н. Тихонова о сингулярно возмущенных системах уравнений Пусть задана система уравнений с непрерывно дифференцируемыми по всем своим аргументам во всей области рассмотрения правыми частями, включающая малый параметр 0, определяющий разницу двух временных масштабов:

относительно переменных x Rn, y Rm, называемых соответственно быстрой и медленной для системы (4.1). Таким образом, изначально предполагается известным разбиение всего набора фазовых переменных на медленные и быстрые в соответствии с характерными масштабами скоростей их изменений. При этом в системе (4.1) t считается быстрым временем, так что при подстановке в нее = 0 получаем так называемую быструю систему:

Здесь и далее f (x, y) = F (x, y, 0), g(x, y) = G(x, y, 0). Введение медленного времени = t позволяет переписать систему (4.1) в медленном виде:

в отличие от его быстрого вида (4.1). Подставляя = 0 в (4.3), получаем так называемую медленную систему:

Множество ее решений принадлежит медленной поверхности S Rn Rm, определяемой первым уравнением в (4.4). Она составляет множество положений равновесия быстрой системы. Пусть S S – та ее часть, на которой эти положения равновесия устойчивы в линейном приблиx,y) Более детальное рассмотрение позволяет уточнить результат и привести некоторые оценки сходимости. Пусть l() = | ln()|, l(0) = 0. Оказывается, что можно подобрать такие положительные не зависящие от константы Ci, что при достаточно малых > 0 решение системы (4.3) с начальным условием (x1, y1) U таким, что x1 () = (y2 ) и |y1,2 y0 | < C1 l() будет отличаться от решения медленной системы (4.4) с начальным условием ((y0 ), y0 ) не более чем на C2 l() на полуинтервале времени (C3 l(), T ]. Содержательный смысл функции l() определяется как время релаксации быстрой системы (4.2), необходимое для попадания ее решения с некоторым конечным начальным условием в -окрестность медленной поверхности, пересчитанное в медленном времени. В предположении линейной устойчивости рассматриваемых стационарных решений быстрой системы отклонения от них имеют экспоненциальный порядок убывания, что дает логарифмическую оценку в быстром времени.

Медленные переменные, изменяющиеся со скоростью порядка единицы, за время порядка l() успевают измениться на величину того же порядка.

Замечания.

1. В случае n = 1 медленная поверхность разбивает все фазовое пространство на области притяжения различных компонент S по отношению к быстрой системе (4.2), так что в определении сходимости по медленным переменным вместо последнего условия достаточно потребовать локализации начального значения возмущенной системы (xk, yk ) в той из областей, которая соответствует компоненте, содержащей ((y0 ), y0 ).

Те же соображения относятся и к вопросу о построении оценок. Последнее из условий, накладываемых на начальную точку (x1, y1 ) в случае выполнения остальных эквивалентно указанной локализации и выполнению неравенства |y1 y0 | < C1 l().

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

3) Релаксационные колебания В случае выполнения условий теоремы Тихонова решения системы (4.3) остаются вблизи одной из компонент S медленной поверхности. Однако в случае приближения к границе такой компоненты скорость медленных переменных продолжает оставаться конечной, что влечет попадание на нее за конечное время. Условия теоремы нарушаются, и траектория уже может покинуть малую окрестность поверхности. В случае одной быстрой переменной (n = 1) происходит это следующим образом. По мере приближения решеx,y) ния медленного уравнения к производная fx возрастает до нуля, что означает потерю устойчивости стационарного решения быстрого уравнения. Без ограничения общности можно считать, что траектория достигает границы в точке x = 0, y = 0. В случае общего положения эта точка не является стационарной для медленной системы, так что существует направление, по которому медленные переменные меняются с конечной скоростью. Если выбрать одну из медленных переменных по этому направлению, а остальные выразить через нее согласно уравнению медленной траектории, то без ограничения общности можно ограничиться одной медленной переменной, что позволяет представить медленную поверхность в окрестности точки x = 0, y = 0 как y = axk, a = 0 с k 2 (по поводу отсутствия остаточного члена см. рассуждения, применявшиеся для обоснования соотношения (3.9)). Поскольку dx = xak d, а d = 0, то при x 0 быстрая переменная бесконечно разгоняется в медленном времени dx, что соответствует смене медленного движения быстрым.

Такое явление называется срывом, а сама точка (x, y), на которой оно происходит, – точкой срыва. Для возмущенной системы с ненулевым значением малого параметра такой срыв означает изменение скорости быстрой переменной на величину этого параметра. В результате срыва траектория уходит из окрестности прежней компоненты S и устремляется к некоторой другой такой компоненте или уходит в бесконечность. Куда она устремится, зависит уже от того, в Рис. 4. Релаксационные колебания в системе Ван-дер-Поля области притяжения какой из компонент S оказалась точка срыва. После того, как срыв произошел, снова начинают выполняться условия теоремы Тихонова. Траектория снова достигает медленной поверхности и начинает изменяться в соответствии с ее утверждениями. Чередования медленных (квазистационарных) и быстрых ( ационных) движений называются ационными колебаниями. Это понятие включает как произвольные блуждания с перескоками с одной компоненты S на другую, включая попадание, в конце концов, в положение равновесия или уход в бесконечность, так и собственно колебания, заключающиеся в периодическом движении по замкнутой траектории.

Пример. Уравнение Ван дер Поля.

Релаксационные колебания можно продемонстрировать на примере системы Ван дер Поля Здесь медленная поверхность определяется уравнением y = x3 x, а устойчивость ее элементов в быстром времени – неравенством |x| > 3. Точки ± 3, 33 являются точками срыва, ограничивающими с одной из сторон две компоненты S. На каждой из этих компонент медленное движение направлено в сторону точки срыва, а поскольку их проекции на ось медленной переменной пересекаются, то в системе возникают релаксационные колебания (см. рис. 4).

4) Релаксационные колебания в системе Предположим, что между двумя биологическими видами имеется трофическое взаимодействие. Это означает, что один из видов служит кормом для другого. Первый из них назовем ресурсом, второй – потребителем. Соответствующие им зависимости биомассы от времени будем обозначать как R(t) и N(t). Предположим также, что характерное время изменения биомассы потребителя существенно больше характерного времени изменения биомассы ресурса. Таким образом, R(t) выступает в роли быстрой, а N(t) – медленной переменной.

При отсутствии потребителя будем считать динамику ресурс подчиняющейся уравнениям поиска партнера, но с возможным дополнительным внешним источником. Потребление ресурса осуществляется, как предполагается, с линейной трофической функцией. Такое взаимодействие может быть Рис. 5. Релаксационные колебания в системе «ресурс – описано следующей системой уравнений:

Здесь R выражено в единицах его максимальной численности при отсутствии внешних воздействий, Q – внешний приток ресурса, d – смертность потребителя, выраженная в единицах ресурса. Заметим, что первое уравнение в (4.5) не содержит коэффициентов перед слагаемыми, поскольку от них в случае необходимости можно избавиться перемасштабированием N и t. Обозначая его правую часть как fQ (R, N), находим медленную поверхность из условия fQ (R, N) = 0, точки срыва (складки) – из дополнительного условия QR = 0 и точку сборки из еще одного дополfQ (R,N ) нительного условия = 0. Рассматривая в последR нем случае Q как второй параметр (в дополнение к N), получаем, что сборка возникает при Q = Qc = 27. Для всех 0 < Q < Qc существуют по две точки срыва m = (Rm, Nm ) и M = (RM, NM ), так что для таких Q при Rm < d < RM и малых значениях параметра возникают релаксационные колебания. На рис. 5 по оси ординат отложено произведение NR, так что в этих координатах медленная поверхность совпадает с графиком функции fQ (R, 0). Релаксационные колебания будут происходить в контуре (SM, m, Sm, M ), где SM и Sm – точки пересечения лучей, касающихся графика функции fQ (R, 0) в точках M и m соответственно, с этим графиком.

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

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

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

2) Устойчивость в случае дискретных Пусть f (x) – гладкая вещественная функция вещественРис. 6. Диаграмма Ламерея. Сходимость дискретных ного аргумента. Рассмотрим последовательность точек {xk }, задаваемых по правилу На рис. 6 она изображена в виде проекции на одну из осей последовательности {(xk, xk+1 )}, получаемой построением диаграммы Ламерея, состоящей из направленной цепочки чередующихся вертикальных и горизонтальных отрезков, соединяющих биссектрису первого координатного угла и график функции f (x), причем вертикальные направлены от биссектрисы к графику. Точки их пересечения задают стационарные точки отображения (5.1). Устойчивость стационарной точки x = f (x ) понимается в смысле локализации в любой наперед заданной ее окрестности хвоста последовательности {xk }, стартующей достаточно близко от нее. Поскольку для отклонений zk = xk x выполнены соотношения zk+1 = Azk + o(zk ), где A = dx x=x, то стационарная точка x является заведомо устойчивой при |A| < 1 и заведомо неустойчивой при |A| > 1.

Аналогичное утверждение имеет место и для многомерных отображений (5.1) с xk = (x1,..., xn ). В этом случае в правой части соотношения в отклонениях в качестве линейного оператора A рассматривается якобиан правой части отображения (5.1), вычисленный в стационарной точке x. При этом для сравнения с единицей вместо абсолютной величины используется спектральный радиус r(A) = max(A) ||, где (A) – спектр, т.е. набор собственных значений оператора A.

3) Отображение Пуанкаре для цикла Пусть система уравнений имеет периодическую траекторию (цикл) S Rn+1, содержащую некоторую точку y0 S, причем такую, чтоf (y0) = 0.

В окрестности U(y0 ) этой точки рассмотрим содержащую ее n 1 мерную плоскость, пересекающуюся с S без касания. В качестве примера таковой можно взять = {y :

(y y0, f (y0)) = 0}, где (y1, y2 ) = n+1 y1 y2 – скалярное проll изведение векторов y1 = (y1,..., y1 ) и y2 = (y2,..., y2 ).

Выбрав на систему координат x = (x1,..., xn ) с центром в точке y0, можно построить отображение P : U(y0 ), переводящее точку плоскости x посредством движения по траекториям системы (5.2) за один оборот вблизи периодической траектории S в точку той же плоскости. Такое отображение называется отображением последования, или отображением Пуанкаре. Для него точка x = 0 является стациоРис. 7. Отображение Пуанкаре для фокуса нарной, а свойства ее устойчивости не зависят от конкретного выбора подходящей секущей плоскости. Эти свойства также совпадают со свойствами устойчивости периодической траектории S в смысле локализации положительных полутраекторий возмущенных решений системы (5.2) в окрестностях цикла S.

4) Отображение Пуанкаре для фокуса Пусть система уравнений (5.2) на плоскости, т.е. при значении n = 1, имеет положение равновесия типа «фокус» в начале координат. Тогда, как и в случае цикла, для каждой проходящей через него направленной прямой, образующей с положительным направлением оси абсцисс угол [0, 2), упорядоченной координатой x, такой, что y1 = x cos, y2 = x sin, локально (т.е. в малой окрестности нуля) определено отображение P (x) за один полный оборот вдоль траекторий системы. Это отображение, которое мы также буРис. 8. Иллюстрация доказательства теоремы дем называть отображением последования или отображением Пуанкаре, обладает рядом очевидных свойств: для любого [0, 2) Для его отклонения от тождественного q (x) = P (x) x имеет место следующая Теорема. Если для некоторого > 0 q0 (x) = 0 при x (0, ), то и для некоторого 1 > 0 q (x) = 0 при x (0, 1) и [0, 2).

Доказательство. Если это не так, то найдется угол (0, 2) такой, что для достаточно малого x > 0 P (x) = x. (см.

рис. 8). Рассматривая точку, задаваемую полярными координатами (x, ), в качестве начальной для исходной системы и двигаясь из нее назад и вперед по времени до пересечения траекторий с положительным направлением оси абсцисс, получим, соответственно, два малых положительных значеРис. 9. Отклонение отображения Пуанкаре от тождественного в теореме Ляпунова ния x1 и x2, причем, поскольку P0 (x1 ) = x2, то x1 = x2. Но продолжение второй траектории проходит в силу сделанного предположения через точку (x, ), так что через эту точку в направлении обратного времени проходят, по крайней мере, две различные траектории, что противоречит теореме единственности решений обыкновенных дифференциальных уравнений.

Следствие (теорема Ляпунова). (см., например, [5]) Разложение функции q0 (x) в ряд Тейлора начинается с членов нечетной степени: q0 (x) = l x2l+1 + o(x2l+1 ).

Доказательство. Из доказанной теоремы и приведенных выше свойств отображения Пуанкаре следует, что q0 (x) в окрестности нуля меняет знак, чего не могло бы быть в случае начального члена четной степени.

Коэффициент l = 0, задающий начало указанного разложения, называется lй ляпуновской (или фокусной)величиной.

При этом если l < 0, то фокус является асимптотически устойчивым, а если l > 0, то асимптотически неустойчивым. Таким образом, для определения характера устойчивости фокуса следует определить знак первой из ляпуновских величин l, необращающихся в нуль. Ее знак не зависит от выбираемого для x масштаба, хотя сама она, конечно, имеет размерность, равную [x]2l. Заметим, что случай, когда 0 = 0, не является интересным с точки зрения исследования устойчивости фокуса, поскольку в этом случае она определяется по линейной части системы (5.2), т.е. по ее якобиану, так что необходимость дальнейшего анализа отпадает.

5) Бифуркации при потере устойчивости Без ограничения общности будем считать, что при изменении параметра, обуславливающем потерю устойчивости положения равновесия типа «фокус», рассматриваемое положение равновесия всегда остается в начале координат. В противном случае можно было бы сделать подстраивающуюся под параметр замену фазовой переменной, сдвигая начало координат в положение равновесия. В таком случае рассматриваемое семейство уравнений принимает вид где x = (x1, x2 ), g(x, ) = o ( x ). Последнее соотношение подразумевает равномерность по из некоторой окрестности нуля. Будем считать, что квадратная матрица второго порядка A() при прохождении параметром нулевого значения слева направо теряет устойчивость, связанную с пересечением ее комплексно сопряженными корнями мнимой оси. Это означает, что 2 = det A(0) > 0, trA(0) = 0 и след является монотонно возрастающей функцией параметра. Снова, не ограничивая общности, мы можем считать, что trA() = 2, так что вещественные части корней матрицы A() пересекают мнимую ось с единичной скоростью = 1 d trd A() В конкретных задачах с заданным параметром, обеспечивающим ненулевое значение скорости, следует делать поправку на нее, рассматривая вместо старого параметра новый =.

При сделанных предположениях для малых || вместо (5.3) можно записать где матрица A(), такая, что A(0) = A(0), имеет чисто мнимые корни ±i () = ±i(1 + + o()). Выбирая в качестве единицы времени один период линейной системы T = 2 и () переходя к времени в периодах = Tt, получим умножением (5.4) на T систему d = T A()x + T (x + g(x, )). Матрица T A() подобна матрице A = T0 A(0), причем замена переменных x x, сопрягающая эти матрицы, гладко зависит от и совпадает с тождественной при = 0. В переменных x эта система примет вид dx = Ax + T (x + g (x, )), где g (x, ) = g(x, ) = g(x, )+o(g(x, )). Игнорируя индексы в и выделяя главные члены, получим Здесь при оценке остаточного члена предполагается сходимость || + x 0. Полагая r = T0 и подставляя в (5.5) x = er y, получим для В (5.6) остаточный член уже включает слагаемые, зависящие от времени (ибо e±r = (1 ± r + o(|r |))), так что оценки следует рассматривать на интервалах времени порядка периода.

Предположим, что при = 0 уравнение (5.6) (без остаточных членов) задает некоторое отображение Пуанкаре переменной Y = C1 x1 + C2 x2 :

Отметим, что при l 1 фокусная величина масштабным образом зависит от выбора Y, так что можно считать l = l (C1, C2 ). Учет малых при отображениях за период = 1 можно осуществить, сделав переход от Y к X = er Y для значений и X = Y для аргумента отображения Пуанкаре. При этом для экспоненты используется приведенная выше оценка, а учет влияния малых ненулевых значений параметра на выражение (5.7) осуществляется поправкой остаточного чле на.

В выражении (5.8) остаточный член включает слагаемые вида pq r p X q для случаев p = 0, q > 2l + 1; p 1, q 1, p + q 3. Применяя диаграммы Ньютона (см. приложение А) к уравнению Pr (X) X = 0 с двумя малыми параметрами X и r и двумя крайними точками 11 = 1, 2l+1,0 = l находим, что отображение Пуанкаре Pr (X) имеет нетривиl r альное положение равновесия X = Pr (X ) вида X.

Для определения его устойчивости вычислим производную dPr (X) ражения видно, что его устойчивость противоположна (по крайней мере, для малых значений r) устойчивости тривиального положения равновесия. Именно в докритической области (т.е. при r < 0) устойчивым является тривиальное положение равновесия (производная по абсолютной величине меньше единицы), а нетривиальное – неустойчивым. В закритической же области все наоборот. Однако само нетривиальное положение равновесия существует в докритической области только в случае l > 0, т.е. тогда и только тогда, когда оно неустойчиво.

Наоборот, случаю l < 0 соответствует существование устойчивого нетривиального положения равновесия в закритической области (т.е. при r > 0). В терминах исходной системы дифференциальных уравнений нетривиальные положения равновесия отображения Пуанкаре – это циклы. Первый из рассмотренных случаев соответствует так называемой жесткой потере устойчивости фокуса (которому соответствует тривиальное положение равновесия отображения Пуанкаре), второй – его мягкой потере устойчивости. При этом при r = 0 возникающий сложный фокус неустойчив в первом случае и устойчив во втором, так что для того, чтобы узнать, какая потеря устойчивости имеет место, достаточно определить характер устойчивости сложного фокуса, т. е.

знак первой ненулевой ляпуновской величины l.

Замечание о терминологии. Использование терминов «мягкая» и «жесткая» потеря устойчивости связана с тем, что в случае возникновения устойчивого цикла траектории, покинувшие малую окрестность ставшего неустойчивым фокуса, сходятся к образовавшемуся устойчивому циклу. Его размер Рис. 10. Бифуркация Хопфа для случаев мягкой (слева) и жесткой (справа) потери устойчивости фокуса т.е. убывает до нуля при 0. Поэтому при потере устойчивости фокуса в этом случае будут наблюдаться колебания с амплитудой X (), т.е. тем меньшей, чем меньше параметр отклонился от своего критического значения. Таким образом решение системы (5.3), стартующее недалеко от положения равновесия, по мере возрастания этого отклонения локально, т.е. как бы «мягко» начинает удаляться в колебательном режиме. Для наиболее часто встречающегося при этом случая 1 < 0 величина такого отклонения пропорциональна, что изображено на рис. 10 слева.

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

Этот неустойчивый цикл по мере приближения параметра к своему критическому значению как бы «наваливается» на пока еще устойчивый фокус, сужая до нуля размер области его притяжения (ее диаметр также характеризуется величиl системы под воздействием случайных флюктуаций вырывается за пределы этой области еще на подходе к критическому значению параметра и уходит далеко в сторону к какомуто другому притягивающему режиму. Потеря устойчивости имеет явно нелокальный характер и из-за возникающих резких перестроек в описываемых реальных системах называется «жесткой». Для случая 1 > 0 соответствующая картина бифуркаций изображена на рис. 10 справа.

В заключение приведем выражение для первой ляпуновской величины через коэффициенты тейлоровского разложения правых частей системы до третьего порядка включительно. Пусть система на плоскости (x1, x2 ) имеет вид 2 = bc a2 > 0 (производные здесь и далее вычисляются в положении равновесия (x1, x2 ) = (0, 0)), так что начало координат является положением равновесия типа «центр» по линейному приближению.

Тогда величина L1 = 1 (1, 0), соответствующая выбору Y = x1, может быть вычислена из выражения (см., например, [5]):

L1 = 4b3 {[ac (a3 + a11 b02 + a02 b11 ) + ab (b2 + b11 a20 + b20 a11 ) +c (a11 a02 + 2a02 b02 ) 2ac (b2 a20 a02 ) 2ab (a2 b20 b02 ) b2 (a20 b20 + b11 b20 ) + (bc 2a2 ) (b11 b20 a11 a20 )] (a2 + bc) [3 (cb03 ba03 ) + 2a (a21 + b12 ) + (ca12 b21 b)]}.

Приводимый ниже пример иллюстрирует использование этой формулы.

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

6) Система «хищник – жертва» с унимодальной плодовитостью жертвы Предположим, что в условиях отсутствия хищников внутренняя саморегуляция численности жертв x задается унимодальной функцией (x), достигающей невырожденного положительного максимума в некоторой точке M > 0. В качестве биологической интерпретации такого выбора можно рассматривать случаи наличия проблем с поиском брачного партнера при малых численностях с одновременным самолимитированием при больших численностях жертв. Взаимодействие с хищником y осуществляется в соответствии с гипотезой встреч, т.е. пропорционально численностям обоих видов, что соответствует линейной трофической функции этого взаимодействия. Система уравнений, описывающая такого рода систему, может быть представлена в виде Здесь равенство (с точностью до знака) коэффициентов при произведении xy в обоих уравнениях обеспечивается согласованием масштабов измерения их численностей (считается, что единица съеденной жертвы увеличивает на единицу численность хищника), а естественная смертность хищника m рассматривается в качестве параметра. Функция (x) задает разницу рождаемости и смертности жертвы.

Нетривиальное положение равновесия, в котором представлены оба вида x = m, y = (m), устойчиво при < 0, т.е. при m > M, и теряет устойчивость при пересечении параметром m сверху вниз точки максимума функции (m). В окрестности этой точки указанное положение равновесия является фокусом. При этом, поскольку след якобиана правой части системы (5.11) в положении равновесия равен m (m), то скорость пересечения корнями мнимой оси равна В качестве параметра бифуркации удобно выбрать параметр = (mM) = 1 (M)(mM), меняющийся в направлении потери устойчивости положения равновесия. Для определения характера бифуркации в критической точке m = M достаточно определить, как мы видели, характер устойчивости возникающего в этой точке сложного фокуса. Для этого вычислим первую фокусную величину L1 по формуле (5.10) для x1 = x, x2 = y. Для входящих в эту формулу элементов якобиана имеем: a = 0, b = M, c = (M), так что 2 = M(M).

Из остальных коэффициентов ненулевыми оказываются следующие a20 = 1 M (M), a11 = 1, a30 = 1 (M) + M (M), b11 = 1. В результате подстановки получаем Для амплитуды колебаний находим В частном случае (M) = 0, например для квадратичной функции (x), это дает амплитуду устойчивого предельного цикла X = 2 M(M m) при малых M m > 0.

6. Автоколебания в системах 1) Модель Вольтерра Простейшей моделью, описывающей взаимодействие двух видов по типу «хищник – жертва», является модель, предложенная В. Вольтерра (см., [6]). В основе этой модели лежат гипотезы о том, что динамика участвующих во взаимодействии видов при их изолированном рассмотрении происходит по линейным законам, в то время как само их взаимодействие осуществляется в соответствии с гипотезой встреч. Эта гипотеза, восходящая к законам химической кинетики, заключается в том, что число встреч, приводящих к пожиранию одного вида другим, пропорционально численности каждого из видов. При этом считается, что фиксированная доля биомассы съеденных жертв, пропорциональная их численности, переходит в биомассу хищника.

Пусть x и y обозначают численности жертвы и хищника соответственно. Описывающая их взаимодействие на основе указанных выше принципов система дифференциальных уравнений может быть представлена в виде Здесь и m – положительные коэффициенты, масштабы измерения x и y подобраны так, чтобы коэффициенты при нелинейных членах были единичными (см. комментарии к системе (5.11)). Отметим здесь, что система (6.1) может быть получена из (5.11) при = const. Ее нетривиальное положение равновесия (x, y ) = (m, ) является положением равновесия типа центр по линейным членам с частотой локальных колебаний = m. Подсчет ляпуновских величин дает для них нулевые значения (для L1 это следует из формулы для (5.11)). И это не случайно. Оказывается, что система (6.1) является консервативной, т.е. имеет первый интеграл.

Действительно, из (6.1) находим, что dx + dy = x my.

С другой стороны m d ln x + d dt y = my + x. Сложив эти равенства, можно разделить переменные m d ln x dx = dy d dt y, что дает интеграл вида C = m ln xx+ ln yy, вместо которого удобнее рассматривать постоянную на траекториях функцию Иногда вместо логарифмического выражения (6.2) пишут экспоненциальное yey xex = const. Интегрируя на одном периоде уравнения для логарифмов, можно получить для средних значений (, y ) = (x, y ). Заметим, что сам период возрастаx ет по мере удаления циклов от положение равновесия (x, y ) от T0 = 2 вблизи него до бесконечности (предельной траекторией служит цепочка сепаратрис нулевого седлового положения равновесия с входящей по оси x = 0 и выходящей по оси y = 0).

Наличие первого интеграла (6.2) означает, то положение равновесия (x, y ) системы (6.1) является положением равновесия типа центр не только по линейному приближению, но и для системы в целом, что означает, в частности, негрубость этой системы. Так, например, добавление к первому уравнению (6.1) слагаемого вида x2 для сколь угодно малого > 0, означающего с биологической точки зрения включение в рассмотрение внутренней стабилизации жертв, делает положение равновесия (x, y ) = (m, m) асимптотически устойчивым (неустойчивым при < 0).

Действительно, в этом случае функция (6.2) служит функцией Ляпунова, ибо V (x, y) = xx, yy, так что (x, y ) является ее единственной экстремальной точкой, а вдоль траекторий возмущенной системы она убывает почти всюду, ибо dV = (x x )2, причем для локально ограниченdt ной третьей производной в стационарных неэкстремальных точках функции V (x, y), т.е. при x = x, y = y выполнено d V = 2x2 (y y)2 < 0. В предположении отсутствия асимптотической устойчивости, т.е. в предположении существования ограниченной замкнутой области U = {(x, y) :

V (x, y) > > 0, x > 0, y > 0} ее представление в виде объединения областей U1 = U (x, y) : dV dt < 1 < и U2 = U \ U1 приводит к противоречию из-за невозможности распределения времени по условию пребывания траектории системы (6.1) в каждой из областей. (Локализация траектории в U требует конечного времени ее пребывания в U1, а значит и за пределами U2, где ограниченная третья производная может быть положительной. В таких условиях бесконечно долгое пребывание траектории в области U2 ведет к бесконечному убыванию второй производной, что противоречит предполагаемой ограниченности U.) 2) Трофические функции в системе Трофическим взаимодействием двух видов называется такое, при котором биомасса одного из видов служит ресурсом для выживания и развития другого. К числу таковых относятся взаимодействия по принципу «хищник – жертва», «ресурс – потребитель» в случае, когда ресурсом является живая биомасса какого-либо вида, как правило, растительного, или «паразит – хозяин». Такой тип взаимодействия упорядочивает виды в так называемую трофическую цепь (см., например, [7]). В звене такой цепи вид, служащий ресурсом, считается расположенным на нижнем трофическом уровне, в то время как на верхнем – вид, пользующийся этим ресурсом. Зависимость мальтузианского параметра последнего от численности V (x) первого называется трофической функцией. Так, например, в случае вольтерровской системы «хищник – жертва» (6.1) V (x) = x. Такая функция в силу естественных соображений должна удовлетворять ряду условий. 1) V (0) = 0, 2) V (x) 0, 3) 0 < V (x) < при 0 xmin < x <. В последнем условии xmin обозначает некоторый порог «недоступности» жертвы, ниже которого она находится вне поля зрения хищника. Очень часто в моделях считается xmin = 0. Замена линейной трофической функции в модели Вольтерра (6.1) на произвольную приводит к системе вида Эта система имеет тривиальное положение равновесия, являющейся седлом с входящей по оси y и выходящей по оси x сепаратрисами. Кроме того, в случае V () > m существует еще и нетривиальное положение равновесия (x, y ) = (V 1 (m), V 1 (m)/m). В грубом случае можно считать, что V (x ) > 0, а, значит, такое положение равновесия единственно и является узлом или фокусом (поскольку определитель якобиана системы (6.3), вычисленного в этом положении равновесия, равен x V (x ) > 0). Оно асимптотически устойчиво в случае отрицательности следа этого якобиана, т.е. при выполнении условия V x ) < V (x ). Последнее неравенство всегда выполнено в случае выпуклости V (x) > 0 на интервале (0, x ) (интегрирование по частям функции xV (x) на этом интервале), а также эквивалентно условию монотонного возрастания потребляемой доли (V (x)/x)x |x=x > 0. Неравенства противоположного знака соответствуют асимптотической неустойчивости положения равновесия.

Более того, в случае знакопостоянства функции (x) = (V (x)/x)x при выполнении дополнительного условия «возвращаемости» V () > m+, нетривиальное положение равновесия является глобальным аттрактором в области x > 0, y > 0 системы (6.3) в прямом времени (т.е. при t +) при (x) > 0 и в обратном времени (т.е. при t ) при (x) < 0. Это следует из возвращаемости траекторий в некоторую ограниченную область, а следовательно, и ее ограниченности при движении по одному из направлений времени, отсутствия в указанной области других положений равновесия и циклов и структуры фазовых траекторий автономных систем на плоскости (см. приложение Б). Вопрос о существовании циклов будет рассмотрен ниже.

По поводу возвращаемости отметим следующее. Пусть x таково, что > 1 (при достаточно больших x это неравенство выполнено в силу условия возвращаемости). Все фазовое пространство x > 0, y > 0 системы (6.3) делится ее нуль-изоклинами и нижним отрезком прямой x = x на пять областей: A = {0 < x < x, 0 < y < x/V (x)}, B = {x Рис. 11. Направления траекторий в модели «хищник – Нетрудно видеть (см., например, рис. 11 ), что направление векторного поля в каждой из областей таково, что траектория системы (6.3), проходящая через нее, с необходимостью переходит из нее в следующую в соответствии с порядком A B C D E A. Так, например, переход C D обуславливается тем, что в области C = y(V (x)m) > y(V (x)m) > x, так что d ln x > 1, в то вреd ln y мя как нуль-изоклина переменной x, задаваемая уравнением yV (x) = x, в логарифмическом масштабе при x > x будет близка к прямой, имеющей единичный наклон. Рассматривая отображение последования на границе областей A и B, остается только проверить, при каком направлении времени траектория оказывается ограниченной. Случаю возрастания величины y для этого отображения соответствует положительное направление.

3) Достаточные условия отсутствия циклов Достаточные условия отсутствия циклов системы (5.9) в односвязной области G на плоскости могут быть заданы условиями Бендиксона – Дюлака (см., например, [5]), заключающимися в существовании гладкой функции B(x1, x2 ) такой, что функция (x) = (x, BX) знакопостоянна в области G и может обращаться в нуль лишь на множествах нулевой двумерной лебеговой меры. Здесь X = (P, Q)T – вектор правых частей системы (5.9) x = (x1, x2 ) – градиент по фазовым переменным (·, ·) – скалярное произведение в R2, так что проверяемую функцию можно переписать в виде (x) = div(BX) = x1 (BP ) + x2 (BQ). Достаточность является следствием теоремы Гаусса, ибо в случае существования цикла интегрирование по внутренности цикла дает 0 = (x)dx = (BX, n(s))ds = 0, где n(s) – вектор внешней единичной нормали к в точке s, ортогональный траектории, проходящей по. Отсюда также следует, что в случае выполнения указанного условия исключаются не только периодические решения, но и любые замкнутые цепочки траекторий в указанной области. В частности, не может быть сепаратрисных петель или двух узлов, напрямую связанных траекторией (ибо рядом с такой траекторией должна существовать и другая такая же, а полученная пара траекторий образует цикл).

Заметим, наконец, что в случае B(x, y) 1 условия Бендиксона – Дюлака эквивалентны знакопостоянству следа якобиана системы. Так, например, отсюда следует отсутствие циклов у системы (6.3) в случае знакопостоянства функции (x) = (V (x)/x)x в области {x > 0}.

4) Периодические решения в системе В случае нарушения условия знакопостоянства функции (x) = (V (x)/x)x в системе (6.3) будут в случае общего положения наблюдаться циклы по крайней мере для некоторого интервала значений параметра m. Если для некоторой точки x > 0 имеет место равенство () = 0, причем в этой точке происходит смена знака функции, то в качестве такого интервала можно выбрать одну из полуокрестностей точки m = V (). Ситуация здесь полностью укладывается в описанную выше теорию бифуркации Хопфа. Направление убывания функции (x) в окрестности точки x определяет направление потери устойчивости положения равновесия (x, y ) = (V 1 (m), V 1 (m)/m) при прохождении параметром m значения m в том же направлении. Для величины m, близкой к значению m, положения равновесия будет фокусом, причем при m = m – сложным фокусом. Характер его устойчивости, а следом за ним и характер сопутствующих периодических решений, включающий характер их устойчивости и область их существования среди возможных значений параметра m, определяется знаком первой ненулевой ляпуновской величины. Подсчет по формуле (5.10) первой из проверяемых дает значение Заметим, что здесь V () > 0, в то время как характер второй производной остается неопределенным. При этом в случае выполнения неравенства V () = 0 знак и значение третьей производной также могут быть произвольными, так что знак в (6.4) в принципе может быть любым. Устойчивые периодические решения, соответствующие мягкой потере устойчивости положения равновесия (x, y ), будут наблюдаться в закритической области (т.е. при (V 1 (m)) < 0) для L1 < 0, а неустойчивые, соответствующие его жесткой потере устойчивости, – в докритической (т.е. при (V 1 (m)) > 0) для L1 > 0. В негрубом случае V () = 0 неравенство V () = оставляет функцию (x) знакопостоянной в окрестности точки x, что исключает согласно изложенному в предыдущем параграфе существование циклов в ее окрестности. Если же в этом случае, V () = 0, то тогда и L1 = 0, так что для изучения характера бифуркаций требуется вычисление ляпуновских величин более высокого порядка.

5) Модель Колмогорова В качестве альтернативы модели Вольтерра А.Н. Колмогоров предложил (см. [8]) свою модель системы «хищник – жертва», не обладающую указанным выше свойством негрубости. Его модель имеет вид системы Здесь трофическая функция V (x) обладает обычными свойствами V (0) = 0, V (x) > 0, мальтузианская функция жертвы (x) считается монотонно убывающей по ее численности x с единственным положительным корнем x, а мальтузианская функция хищника K(x) считается не зависящей от его численности y и монотонно возрастающей по численности жертвы с единственным положительным корнем x, так что (0) > 0, (x) < 0, () < 0, K(0) < 0, K (x) > 0, K() > 0. Заметим, что при отказе от требования монотонности (x) в случае K(x) = V (x) m эта модель совпадает с моделью (6.3), так что глобальная структура ее фазового портрета позволяет судить о такой же структуре и для рассмотренной выше модели.

Структура же эта такова. Система (6.5) имеет три положения равновесия. Тривиальное положение равновесия (0, 0) является седлом с входящей по оси y и выходящей по оси x сепаратрисами. Также седлом является и равновесное положение жертвы в отсутствии хищника (, 0). Здесь устойчиx вым направлением является ось x, в то время как исходящая из него траектория касается прямой (1 2 )(x x) = V ()y, сюда, в частности, следует, что при x > x исходящая сепаратриса наклонена влево. Третье положение равновесия (x, y ) с y = (x )x имеет обе положительные компоненты в случае x > x. Поскольку определитель якобиана системы (6.5) в этом положении равновесия положителен, то свойства его устойчивости определяются исключительно знаком следа, т.е. величиной a = (0 V0 + x1 V0 x0 V1 )/V0. Здесь и ниже x = x, fi = d f (x ), f (x) = (x), V (x), K(x). Случаю a = 0 соответствует сложный фокус, для определения характера устойчивости которого следует определить знак первой ляпуновской величины. Она равна L1 = 16 (K2 (2 xV0 (xV1 V0 ) + 1 V0 (2xV1 x2 V2 2V0 )) +K1 (22 (3V02 4xV1 V0 + xV12 ) + 21 (2V1 (xV1 V0 ) + x2 (V3 V V1 V2 )) 2x3 V0 (xV1 V0 )))/ xK1 V0 (xV1 V0 ) 1 V0 K Здесь мы исключили 0 из условия a = 0, так что знаменатель в (6.6) является отрицательным (в силу положительности выражения под знаком корня). Из (6.6) видно, что, не нарушая наложенных на задаваемые функции (x), V (x), K(x) ограничений, всегда можно добиться любого знака у L1, хотя бы только за счет подходящего выбора третьих производных.

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

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

раздел 4). Независимо от причины быстрые процессы установления, вызванные конкуренцией, можно рассматривать в рамках одних и тех же моделей, основной целью которых является отыскание условий, обеспечивающих устойчивость стационарных состояний, соответствующих тому или иному набору выживших видов. Основные характеристики таких условий обнаруживаются уже на модели двух конкурирующих видов. Такая модель в простейшем случае записывается в виде системы двух уравнений (см., например, [9]):

Здесь x = (x1, x2 ) – вектор численностей видов, участвующих в конкуренции, f (x) = (f1 (x1, x2 ), f2 (x1, x2 )) – вектор соответствующих им мальтузианских функций. Предполагается, что f (0) > 0 (покомпонентно) и что все элементы якобиана A = f = aij i = 1, 2, отрицательны. Последнее из этих условий подразумевает наряду с межвидовой конкуренцией также наличие и внутривидовой. Более того, обычно предполагается, что все эти элементы ограничены и отделены от нуля > |aij | > > 0, так что нуль-изоклины i = {(x1, x2 ) : fi (x1, x2 ) = 0} имеют конечный отрицательный ненулевой наклон что обеспечивает их однократное пересечение с каждой из положительных координатных полуосей. Пусть fi (i ) = 0, i = 1, 2, где x = (1, 0), x = (2, 0). Система (7.1) имеет, по меньшей мере, три положения равновесия: одно в начале координат x0 = (0, 0) – неустойчивый узел и два на осях координат x1 и x2, имеющих как минимум одно устойчивое направление вдоль своей оси. Точки пересечения нуль-изоклин в области положительных значений x > 0 дают набор положительных положений равновесия системы (7.1) S = {k = x Множество S предполагается конечным, хотя может оказаться и пустым. Устойчивость каждого из нетривиальных положений равновесия системы (7.1) определяется по ее линейным членам. Якобиан системы в каждой из точек множества S равен diag(k ) · A, где по определению В силу отрицательности следа устойчивости соответствует положительность определителя якобиана, т.е. положительность det A, задаваемая неравенством Производные, стоящие в левой части (7.3) имеют естественную биологическую интерпретацию интенсивностей внутривидовой регуляции, в то время как производные, стоящие в правой части – интенсивностей межвидовой регуляции. Таким образом, равновесное состояние, включающее оба вида с ненулевой численностью, оказывается устойчивым тогда и только тогда, когда произведение интенсивностей внутривидовой регуляции превышает произведение интенсивностей межвидовой регуляции.

Неравенство (7.3) можно трактовать и по-другому – в терминах наклонов нуль-изоклин в точке их пересечения. В силу (7.2) его можно переписать в виде:

что означает локализацию кривой 2 выше кривой 1 справа от рассматриваемой точки их пересечения в декартовой системе координат (x1, x2 ). Геометрическая интерпретация, Рис. 12. Нуль-изоклины и чередующиеся по устойчивости положения равновесия в модели конкуренции двух видов исключающая зависимость от порядка нумерации видов, такова: устойчивое положение равновесия расположено выше прямой, соединяющей точки пересечения касательных к кривым i в рассматриваемой точке с соответствующими им координатными осями xi. Поскольку в случае общего положения (т.е. при пересечениях без касания) при последовательном переходе от одной точки пересечения к другой происходит чередование взаимных расположений кривых, то и устойчивые положения равновесия будут чередоваться с неустойчивыми. Что же касается двух крайних положений равновесия xi,i = 1, 2, то устойчивость каждого из них определяется знаком соответствующего значения fj (i ), j = i, причем отx рицательной величине соответствует устойчивость. Этот результат является чрезвычайно важным с точки зрения теории эволюционного отбора, и мы отложим его интерпретацию до следующего параграфа, где он будет обобщен на многомерный случай. С точки зрения геометрической интерпретации устойчивость означает, что точка xi расположена дальше от начала координат вдоль оси xi, чем точка пересечения j, j = i с этой же осью. Отсюда следует, в частности, что ближайшая к xi точка пересечения обеих кривых окажется неустойчивой, так что обе точки xi, i = 1, 2 могут быть включены в построенную цепочку положений равновесия с сохранением характера чередования устойчивостей в качестве ее крайних точек (см. рис. 12).

Наконец, переходя в системе (7.1) к логарифмам, можно убедиться, что след якобиана системы в новых переменных будет всегда отрицательным, что позволяет сделать вывод об отсутствии у нее циклов, а значит и об отсутствии циклов у системы (7.1) в области положительных значений переменных. Поскольку область G = {(x1, x2 ) : 0 x1 X1, x2 X2 } при достаточно больших X1, X2 (они определяются из условий f1 (x1, x2 ) < 0, f2 (x1, x2 ) < 0 при x1 > X1, x2 > X2, так что, должны выполняться, по крайней мере, неравенства x1 < X1, x2 < X2 ) является инвариантной по отношению к системе (7.1), а за ее пределами в области положительных значений переменных отсутствуют положения равновесия, то любая траектория с положительными начальными условиями сходится к одному из описанных выше равновесных состояний (см. приложение Б).

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

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

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

С математической же точки зрения эти установившиеся состояния есть ни что иное, как устойчивые стационарные решения в подходящей модели. При этом устойчивость должна быть двоякая – как по отношению к возмущениям численностей присутствующих в равновесии видов, так и по отношению к внедрению отсутствующих, т.е. имеющих нулевую численность. Первый из этих двух типов устойчивости называется внутренней устойчивостью, второй – внешней. Что же касается вопроса об основных характеристиках видов, разделяющих их в рассматриваемом установившемся состоянии на присутствующих (т.е. живых) и отсутствующих, то принципиальный ответ на него был дан еще в дарвиновском постулате о выживании сильнейших. Отметим здесь, что сама формулировка постулата в форме превосходной степени подразумевает использования для установления отличительных характеристик некоторых экстремальных принципов. Что следует считать «силой» вида, по которой производится отбор, можно установить на основе применения математической модели (микроэволюционного) конкурентного отбора.

Пусть процесс конкуренции n видов с вектором плотностей x = (x1,..., xn ) описывается системой уравнений:

Предположим без ограничения общности, что в некотором устойчивом стационарном состоянии x этой системы присутствуют только первые m ее видов x = (1,..., xm, 0,..., 0); xi > 0, i = 1,..., m. Тогда, во-первых, а, во-вторых, якобиан системы (7.5) в этом положении равноAB весия имеет следующую блочную структуру, где странство численностей выживших видов, B = fi (), i = x 1,..., m, j = m + 1,..., n, 0 – нулевая m r-матрица (r = n m), D = diag (fi ()), i = m + 1,..., n. Спектр такого якобиана локализован в левой комплексной полуплоскости тогда и только тогда, когда там же локализованы спектры каждой из матриц A и D. Такая локализация спектра матрицы A означает устойчивость положения равновесия x по отношению к возмущениям ненулевых его численностей, т.е.

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

Таким образом, необходимым и достаточным условием устойчивости положения равновесия x является одновременное выполнение условий его внутренней и внешней устойчивости. Но матрица D, определяющая внешнюю устойчивость, является диагональной, так что последнее условие заключается в выполнении совокупности неравенств fi () < 0, i = m+1,..., n. Необходимыми для их выполнения с учетом (7.6) являются равенства где мы обозначили N = {i = 1,..., n} – полный набор учитываемых в модели видов, а M = {i = 1,..., m} – набор реально присутствующих, т.е. имеющих ненулевую численность в рассматриваемом положении равновесия. Равенство (7.7) имеет характер экстремального соотношения и называется принципом эволюционной оптимальности. Его биологический смысл заключается в том, что выжившие в том или ином состоянии виды обязаны иметь максимальные значения мальтузианских коэффициентов среди всех потенциально допустимых, вычисленных в этом состоянии. Именно эти коэффициенты характеризуют силу видов в ее дарвиновском понимании.

Заметим, что на функции fi (x), входящие в правые части системы (7.5), не накладывалось никаких специальных ограничений, учитывавшихся, например, в предыдущем параграфе. Более того, также как и в случаях других известных (например, в физике) экстремальных законов, проверка выполнимости условия (7.7) не требует даже их гладкости. Поэтому соотношения типа (7.7) относятся к категории «универсальных» принципов, которые можно использовать в качестве фундамента при построении любой биологической теории, так или иначе апеллирующей к учету факторов эволюционного отбора. Так, обобщая соотношения (7.7) на случай бесконечного числа потенциальных видов можно в качестве параметров оптимизации вводить векторные переменные с компонентами, соответствующими реальным возможным характеристикам видов (например, считать, что Rk – некоторое подмножество k-мерного пространства характеристик). В таком случае количество уравнений, задаваемых соотношениями (7.8), будет в грубом случае в точности равно размерности векторного пространства характеристик. Это означает, что решение экстремальной задачи (7.8) позволит выделить из континуального набора возможных решений некоторый дискретный (чаще всего – конечный) набор допустимых состояний, которые только и следует принимать во внимание при дальнейших рассмотрениях.

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