CC BY
35
СИНХРОНИЗАЦИЯ 2-ЦИКЛОВ В СИСТЕМЕ СИММЕТРИЧНО СВЯЗАННЫХ ПОПУЛЯЦИЙ, ЗАПАС-ПОПОЛНЕНИЕ В КОТОРЫХ ОПИСЫВАЕТСЯ ФУНКЦИЕЙ РИКЕРА
М.П. Кулаков, Е.Я. Фрисман
Исследуются системы связанных отображений, построенных на основе модели Рике-ра и описывающих пространственную динамику неоднородных популяций, представленных двумя связанными группировками особей, между которыми возможно небольшое миграционное взаимодействие. Рассматриваются бифуркационные механизмы синфазной и противофазной синхронизации мультистабильных режимов подобных систем. Для идентификации достигаемых режимов синхронизации вводится количественный показатель синхронизации.
Ключевые слова: Популяция, миграция, система связанных отображений, бифуркация, синхронизация, мультистабильность, бассейн притяжения.
Введение
Для большинства видов животных характерно неравномерное распределение особей по ареалу обитания. Такая неоднородность, будь она вызвана неоднородным распределением трофического (пищевого) ресурса, свойствами ареала (например, неоднородным ландшафтом), сложным типом нелинейной динамики (например, эффектом Олли) или пространственным перераспределением особей по однородному ареалу (для популяций, разные особи которых значительно различаются по своей подвижности), приводит к мозаичности пространственной структуры популяции. Изучению вопросов возникновения, сохранения и количественной динамики такой мозаичной структуры посвящено множество работ. Традиционно динамика подобных структур рассматривается по аналогии с физико-химическими структурами. Как правило, популяция рассматривается как диссипативная, непрерывная (как во времени, так и в пространстве) система, в которой из относительно однородного распределения вследствие мобильности и ограниченного радиуса активности особей образуются конгломераты - скопления особей - и устанавливается
некоторая мозаичность с диффузионной связью между близлежащими скоплениями особей. Вследствие неравновесности системы и перераспределения особей установившаяся мозаичность будет непостоянной как по форме, так и по плотностям скоплений [1].
Зачастую пространственное представление о реальных биологических популяциях формируется на основе достаточно большого числа «точечных» данных о плотностях, отмеченных на определенных учетных территориях, которые по каким-либо критериям рассматриваются как эталонные. При достаточной плотности на учетных территориях можно считать их очагами скоплений особей, а особей, находящихся вне скоплений, считать особями, совершающими сезонное миграционное движение между подобными очагами [2]. В такой ситуации обычно предполагается, что все очаги мозаичности определены и постоянны в пространстве и могут быть представлены как некоторые цепочки (линейные или пространственные) с глобальными связями на основе миграционных потоков [3-12]. Для описания динамики таких систем можно использовать набор логистических отображений, в которых каждый элемент связан с некоторыми другими элементами. Такая связь, независимо от вида отображения, может быть представлена комбинацией двух типов связи: инерциальной и диссипативной. Подобные системы хорошо изучены, обладают богатой феноменологией: здесь выявлены эффекты мультистабильности, синхронизации, хаоса, перемежаемости и др., а также явление кластеризации - образование в процессе развития групп элементов, состояния которых полностью совпадают [3,8,12].
С биологической точки зрения вопрос обнаружения, изучения механизмов и количественной идентификации пространственной синхронизации представляется наиболее интересным. Так, две слабо связанные популяции, динамика которых демонстрирует синфазное поведение, могут быть приняты за одну; слабое же изменение параметров связи может привести к возникновению противофазных колебаний и выявить принципиальное различие в динамическом поведении популящий, не отражаясь при этом на динамике их суммарной численности.
1. Уравнения динамики численности двух точечных популяций
Рассмотрим динамику простейшего из возможных случаев, когда имеющаяся пространственная неоднородность привела к формированию двух местообитаний, расселение между которыми возможно только один раз за сезон, что характерно для субпопуляций, сильно разнесенных пространственно. Предположив, что значимыми для развития популяции являются процессы размножения и «зимовки», можно утверждать, что расселение (миграция) возможно либо до периода размножения, либо после него.
Обозначим численность в первом местообитании (первая субпопуляция) в начале п-го сезона размножения через Хп, а во втором местообитании (вторая субпопуляция) - через Уп. В случае, когда местообитания полностью изолированы и между ними отсутствуют какие-либо миграционные обмены, численность двух этих субпопуляций может быть описана следующими рекуррентными уравнениями:
Хп+1 = /((Ч, Хп), Уп+1 = /(а2,Уп), (1)
где /- некоторая функция воспроизводства, равная сумме числа родившихся с уче-
том младенческой смертности, и выживших особей родительской части популяции; a\ и a2 репродуктивные потенциалы соответствующих субпопуляций. При таком задании процесса воспроизводства явно не учитывается возрастная структура и считается, что по истечении одного сезона после размножения особи не различаются по репродуктивным функциям. Подобный сезонный характер процессов воспроизводства наблюдается, например, у ряда насекомых, некоторых видов рыб, большинства видов мелких млекопитающих.
Перейдем теперь к описанию динамики миграционно связанных популяций. Из всех возможных вариантов связи между субпопуляциями рассмотрим миграции, происходящие непосредственно после размножения. Если процессы воспроизводства в местообитаниях определяются уравнениями (1), а доли особей, мигрирующих из местообитаний, равны mi(Xn) и m2(Yn), то динамика численностей субпопуляций может быть описана известным отображением с симметричной диффузионной связью
Xn+1 = f (ai,Xn)(l - mi(Xn))+ m2(Yn)f (a2,Yn),
(2)
Yn+i = f (a2,Yn)(1 - m2(Yn)) + mi(Xn)f (ai,Xn), симметричность в которой означает, что введенная связь не меняет общей численности, то есть сумма численностей в любой период размножения не зависит от доли мигрантов.
Аналогично может быть получена система уравнений динамики субпопуляций в случае, когда миграция происходит перед процессом размножения,
Xn+i = f (ai,Xn (1 - mi(Xn)) + m2(Yn)Yn), Yn+i = f (a2,Yn (1 - m2(Yn)) + mi(Xn)Xn).
ЗаменаXn = Xn (1 - mi(Xn))+Ynm2(Yn) и Yn = Y,n (1 - m2(Yn))+Xnmi(Xn), которая фактически меняет местами моделируемые фазовые переменные (численность после процессов воспроизводства и численность, наблюдаемую после сезонного расселения), показывает, что эта система тождественна системе (2).
2. Случай постоянной доли мигрантов
В большинстве работ [3,7,8,13-15], так или иначе связанных с изучением системы (2), для описания динамики локальных систем используются однопараметри-ческие отображения (квадратичное Фейгенбаума, кубическое) с областью значений [0; 1] или [-1; 1]. Несложно убедиться, что при подобном способе описания связи, если начальная точка лежит в этой области, а коэффициент связи m не отрицателен и не превосходит m = 1, то значения динамических переменных и число мигрантов для всех последующих моментов дискретного времени окажутся также лежащими в этой области.
Наиболее изучен случай равенства репродуктивных потенциалов (ai = a2 = a) и долей мигрирующих особей (mi = m2 = m = const)
Xn+i = f (Xn) + m (f (yn) - f (xn)),
yn+i = f (yn) + m (f (Xn) - f (yn)),
демонстрирующий богатую феноменологию - переход к хаосу через бифуркации удвоения периода, эффекты мультистабильности, синхронизации, хаос, перемежаемость и прочее. Для этой системы установлены соотношения подобия на основе законов подобия Фейгенбаума, описаны бифуркационные механизмы синхронизации динамических переменных и др. [7,8,13].
Здесь мы рассмотрим особенности динамики системы (3) при допустимых вариациях репродуктивного потенциала a и доли мигрантов m для несимметричного отображения f, основанного на модели запас-пополнение Рикера Xn+i = = aXne-bXn, которое при соответствующей нормировке x = bX сводится к од-нопараметрическому отображению
Xn+i = axne-xn. (4)
При f (x) = axe-x система (3) имеет вид
xn+1 = (1 — m)axne-Xn + mayn e-yn,
(5)
yn+i = (1 — m)ayne yn + maxne Xn
и описывает динамику симметрично связанных локальных популяций, воспроизводство которых происходит по модели Рикера. Как и одиночное отображение (4), система (5) имеет единственную нетривиальную неподвижную точку x = y = ln(a) = h.
Будем в дальнейшем называть синфазными реализациями (синфазным аттрактором) системы (5) асимптотически равные ее динамические состояния (итерации), то есть xn = yn при n — с. Соответственно, если динамические состояния асимптотически не равны одновременно, то есть xn = yn при n — со, будем говорить о несинфазных реализациях (5) (несинфазных аттракторах). Как известно, для систем вида (3) с унимодальной функцией воспроизводства f (x) полностью синфазные ее реализации могут быть описаны одиночным отображением xn+i = f (xn), а все динамические режимы реализуются на биссектрисе первого координатного угла вокруг нетривиальной неподвижной точки, что очевидно выполняется и для системы (5). Все несинфазные же режимы (5) располагаются симметрично вокруг отрезка y = 2h — x, x G [0;2h].
Учитывая, что для модели Рикера f '(h) = 1 — h, несложно получить характеристическое уравнение системы (5) в окрестности нетривиального стационара
А2 — 2А(1 — h)(1 — m) — (1 — h)2(1 — 2m) = 0
с корнями 1 — h и (1 — h)(1 — 2m). Откуда понятно, что нетривиальная неподвижная точка системы (5) устойчива там же, где и аналогичная точка отображения (4), то есть при 0 < h < 2, а дальнейшая эволюция синфазных режимов проходит по классическому сценарию перехода к хаосу через удвоение периода.
Численные расчеты показывают, что несинфазные режимы (xn = yn при n — — с) реализуются вне области, определяемой неравенством |(1 — h)(1 — 2m) | < 1. Причем на линиях h = 2(1 — m)/(1 — 2m) (при 0 < m < 0.5) и h = 2m/(2m — 1) (при 0.5 < m < 1) мультипликатор системы (5) равен (—1). При переходе параметров через эти линии в системе помимо синфазных могут достигаться различные
несинфазные и противофазные аттракторы, стремление к которым определяется начальной точкой итерирования. Это свойство называют мультистабильностью. Оно заключается в существовании в системе при одном наборе параметров различных динамических режимов - различных аттракторов с разной степенью синхронизации (в частности, циклов разной длины), реализуемых в зависимости от начального состояния системы.
Для определенности далее придерживаемся классической классификации колебаний системы (5) [13]. Циклы, обозначенные Мп - это такие колебательные режимы длины N, которые могут быть реализованы ровно N способами. Друг от друга они отличаются временным сдвигом колебательных реализаций системы (3) (решений) на количество итераций, равное п (п = 0,1, 2, ...N — 1). Так, циклы с п = 0 соответствуют синфазным колебаниям системы, циклы с п = 0 - несинфазным.
Известно, что каждый из этих Nn возможных мультистабильных режимов систем (3) и (5) возникает при определенной начальной численности, поэтому в фазовой плоскости можно выделить, по крайней мере, N подмножеств (бассейнов), соответствующих тому или иному аттрактору (если реализуются циклы только длины N). Все эти аттракторы будут различаться между собой реализуемыми на них колебательными режимами, а именно, возникающими на них режимами синхронизации. Для того чтобы отличать один аттрактор от другого, мы предлагаем вычислять для каждого из этих аттракторов количественный показатель, характеризующий достигаемый режим синхронизации. На основе этого показателя синхронизации разработана методика выделения на фазовой плоскости бассейнов аттракторов, соответствующих разным режимам синхронизации.
Для количественной характеристики и регистрации синхронизации используют различные показатели. Популярны коэффициенты, отражающие стабильность разности фаз колебаний связанных осцилляторов [16-18]. Для системы (3) при квадратичной или кубической функции воспроизводства f вопрос количественного описания синхронизации сводится к вопросу устойчивости полностью синхронного режима хп = уп, определяемой трансверсальным Х^ и тангенциальным Хц показателями Ляпунова
вычисляемыми для достаточно большого числа временных реализаций на аттракторе (более 10 тысяч), необходимых для получения малой погрешности [3,12,15].
Среди всех возможных показателей для нас особый интерес представляют такие показатели, которые позволяют «быстро» оценить меру синхронизации двух связанных осцилляторов по существующей динамике. Так, можно вычислять нормированную разницу между полной и безусловной энтропией, подсчитанной по реализа-
3. Количественная характеристика синхронизации
циям любых двух осцилляторов
Бх Б,
х\у
Бх
где Бх - плотность вероятности, подсчитанная по реализациям х, Бх\у - плотность вероятности, подсчитанная по реализациям хп, при условии, что состояние второго осциллятора равно у [16]. Чтобы избежать утомительных оценок плотности вероятности, можно воспользоваться суммарной ошибкой синхронизации
1
N
8 = х — у*\>
i=0
которая стремится к нулю для полностью синхронной траектории или аттрактора [10].
Нормировав эту величину по Xi и вычисляя ее только для (Ж — п) точек аттрактора, получаем показатель
1
N - п
N
Е
(6)
который равен нулю для полностью синхронного аттрактора, отличен от нуля для несинфазного аттрактора, близок к нулю в случаях частичной синхронизации. Более того, у двух разных аттракторов с разной степенью синхронизации значения величины (6) будут различны. Также не сложно убедиться, что величина (6) тем больше, чем сложнее устроен несинфазный аттрактор.
Варьируя начальную точку хо, уо, можно получить некоторые области фазовой плоскости, соответствующие различным начальным условиям системы (5), которые приводят к конкретным режимам синхронизации, то есть построить бассейн притяжения аттракторов (рис. 1 и 2). Причем очевидно, что бассейнам синфазных аттракторов будут соответствовать области, в которых кх = 0, а бассейнам несинфазных -
Рис. 1. Линии уровней функции кх(хо,уо) реализаций (5) при а = 13.0, т = 0.025; соответствующие режимы синхронизации: а, б - синфазные; в, г - противофазные, но в разных фазах
О
к
х
кх > 0. В свою очередь, бассейн несинфазных режимов будет состоять из N — 2) областей, каждой из которых соответствует свой режим синхронизации. Однако это не гарантирует того, что для двух «близких» аттракторов (со схожими режимами) величины (6) не окажутся равными (с учетом численного округления), поэтому в дальнейшем ограничимся применением показателя синхронизации (6) для количественной идентификации и анализа циклов небольшой длины.
На рис. 1 приведена «карта» уровней величины Кх в зависимости от начальной точки итерирования системы (5). Темным участкам соответствуют области фазового пространства, из которых траектории системы (5) всегда приходят к синфазному режиму, а белым участкам - области, траектории из которых приходят к несинфазному режиму. На рис. 1 выделены два принципиально разных колебательных состояния: синфазный 4-цикл (4о) и противофазный 2-цикл (2х); так же ярко выделяются те же циклы, но с иной фазой колебания, полученные одновременным сдвигом итераций в обоих осцилляторах на одну итерацию для 2-цикла и одну, две и три итерации для 4-цикла. Причем для последних оказалось, что величины Кх будут незначительно отличаться друг от друга. Подобное свойство величины (6) происходит за счет того, что при ее вычислении используются точки как аттрактора, так и предшествующей ему траектории, которые будут немного иными для аттракторов с разными фазами, так как величины N и п фиксируются для всей фазовой плоскости при получении среза (6). Различие фаз синфазных режимов на рис. 1 отражается в небольшом различии цветов в синфазном бассейне.
На рис. 2 приведены некоторые бассейны притяжений, иллюстрирующие трансформацию бассейна притяжения полностью синфазного режима. На рис. 2, а, б свет-
Рис. 2. Трансформация бассейнов притяжений аттракторов
лые участки, симметричные биссектрисе первого координатного угла, соответствуют бассейну синфазных колебаний 2о, черные - полностью противофазных - 2\. После второй бифуркации синфазного 20-цикла, его разрушения и формирования устойчивого 4о-цикла, в системе (5) постепенно формируются сложно организованная иерархия циклов длины 4 и 2 разных фаз, отмеченных оттенками серого. Это отражается в постепенном дроблении синфазного бассейна, сначала как показано на рис. 2, в, а затем, при последующем удвоении периода и формировании 8о-цикла, как на рис. 2, г. Дальнейшее изменение мультистабильных состояний и усложнение бассейнов притяжений происходит после разрушения противофазного цикла 21, формирования квазипериодической динамики и ее последующей эволюции до несинфазных циклов разных длин, как показано на рис. 2, д. На последнем рисунке приведен случай полностью хаотической динамики, в котором полностью отсутствуют регулярные синфазные режимы, однако фиксируется сложная фрактальная структура иерархии несинфазных хаотических режимов динамики системы (5).
Используя введенный показатель синхронизации, рассмотрим некоторые би-фуркационые механизмы формирования и разрушения синфазной динамики.
4. Особенности формирования синфазных и противофазных 2-циклов
Очевидно, что эволюция аттракторов системы (5) качественно подобна эволюции системы (3), основанной на квадратичной зависимости и описанной в цитируемых работах, поэтому не будем подробно на ней останавливаться. Наиболее интересными для нас являются механизмы формирования синфазных и противофазных циклических режимов. Рассмотрим их на примере циклов длины 2.
Для определения элементов 2-цикла дважды применим отображения (3) для точки Хп и Уп
Хп+2 = F (Г (Хп, Уп) , С (Хп, Уп)) = Р2(Хп, Уп),
(7)
Уп+2 = С (Г (Хп, Уп) , С (Хп, Уп)) = С2(Хп, Уп),
где Г(х,у) = I(Хп) + т (I(Уп) - I(Хп)), С(х,у) = I(уп) + т (I(Хп) - I(Уп)). Тогда элементы 2-цикла х = хп+2 = хп и у = уп+2 = Уп можно вычислить из системы
у = 2тт 111 (х)+т (1 (у) -1 (х))] - тт-гx,
(8)
х = тт 11 (у)+т (1 (х) -1 (у))\ - тт-~ у-
Решая численно каждое из уравнений (8) относительно одной из переменных (например, перебирая х, находим соответствующие значения у, удовлетворяющие первому уравнению, и симметричные пары переменных, удовлетворяющие второму уравнению), получаем две кривые на фазовой плоскости, симметричные относительно биссектрисы первого координатного угла, пересечением которых являются искомые элементы 2-циклы (рис. 3).
Рис. 3. Эволюция ненулевых решений (8) при уменьшении коэффициента связи m; f (x) = axe x, a = 8; элементы 2-цикла соответствуют точкам пересечения кривых
С другой стороны, решение системы, соответствующее полностью противофазному циклу 2i, удовлетворяет соотношениям xn = yn+i, Уп = xn+1 и xn+2 = = xn = x, yn+2 = yn = У, откуда с учетом (3) и (7) элементы противофазного цикла можно быстро найти из системы
1 - 2m
x = --f (У) +
m
1-m 1 — 2m
y
f (x) +
1 — m' m
(9)
-x.
1 - т 1 - т
Легко проверить, что системе (9) будут удовлетворять тривиальное х = у = 0 и нетривиальное равновесие х = у = Н, а также значения х и у, соответствующие как синфазным, так и несинфазным элементам цикла 2 системы (3), а при I(х) = ахе-х - и системы (5). Несложно показать, что (8) и (9) имеют в качестве своего решения несинфазные элементы систем (3) и (5) тогда и только тогда, когда Н > 2(1 — т)/(1 — 2т) (т < 0.5), то есть вне области устойчивости.
Проследим, как происходит изменение решений (8) и (9) при I (х) = ахе-х. В случае а < е2, то есть при устойчивом поведении системы кривые, получае-
мые из (8) и (9), пересекаются в точках, соответствующих стационарным решениям х = у = 0 и х = у = Н. С ростом репродуктивного потенциала а спектр решений (8) существенно расширяется. На рис. 3 изображены кривые, полученные из (8) после исключения нулевого решения.
В результате первой бифуркации (переход через Н = 2) в спектре решений (8) имеются три ненулевых точки - ставшая неустойчивой неподвижная точка X = у = 1п(а) = Н и элементы устойчивого синфазного 2-цикла (рис. 3, а). При уменьшении коэффициента связи т, по мере того как параметры попадают в область Н > 2(1 — т)/(1 — 2т), в множестве корней (8) появляется еще пара решений, соответствующая элементам неустойчивого противофазного цикла длины два (рис. 3, б). При дальнейшем уменьшении параметра т, этот 21-цикл становится устойчивым. Бифуркация сопровождается появлением еще двух пар решений (рис. 3, в-г), соответствующих неустойчивым несинфазным 2-циклам с частичной противофазной синхронизацией (без захвата амплитуды) Появление этих циклов приводит к возможности противофазной синхронизации 21, как показано на рис. 4, где помимо кривых, полученных из (8) и (9), показана траектория системы (5) (рис. 4, а, в) и увеличенный фрагмент фазовой плоскости (рис. 4, б-г), включающий элемент цикла 21 и неподвижную точку. Рис. 4 иллюстрирует процесс синхронизации и противофазной синхронизации динамики системы (5) соответственно случаям, когда цикл 21 существует, но неустойчив (см. рис. 4, а-б), и когда он становится устойчивым (см. рис. 4, в-г).
Рис. 4. Формирование устойчивого противофазного 2-цикла (21), при а = 10: а, б - неустойчивый цикл 2ь в, г - устойчивый цикл 21
Для определения параметрических границ устойчивости 2-циклов достаточно рассмотреть устойчивость противофазных режимов 21, так как эволюция синфазных состояний для системы (5) аналогична сценарию переходу к хаосу в одиночном уравнении (4).
Уравнения 1 + Б (а, т) + 3(а, т) = 0, 1 — Б (а, т) + 3(а, т) = 0 и 3 (а, т) = 1, построенные на элементах 2-цикла (синфазных и противофазных) системы (5), будут соответствовать бифуркациям их удвоения, седло-узловой бифуркации и бифуркации рождения тора, при которых мультипликаторы системы (7) соответственно равны А* = —1, А* = +1 и комплексно-сопряженным А* = а + гЬ (|А*| = 1). Коэффициенты Б и 3 характеристического полинома (А*)2 — Б А* + 3 определяются как след и детерминант якобиана отображения (7)
Б=
Г + С дх ду
и3
х* ,у*
дГ2 дС2 дГ2 дС2
дх ду ду дх
х*,у*
вычисленные для каждого элемента (х*,у*) исследуемого 2-цикла. Сами же значения динамических переменных, соответствующих 2-циклу, могут быть найдены численно как решения системы (8) или (9), а также могут быть получены в пределе прямым итерированием исходных уравнений (5).
Путем совмещения определенных таким образом границ устойчивости циклов 2о и 21 с картой динамических режимов, использованной для идентификации других циклов, была получена диаграмма существования и качественных переходов некоторых мультистабильных состояний системы (5), изображенная на рис. 5, а. Она представляет собой многолистную поверхность в пространстве параметров, где каждому листу соответствует свой режим синхронизации. На диаграмме отмечены также значения мультипликаторов А* на границах областей устойчивости цикла 21. При построении карты динамических режимов, каждой паре (а, т) параметрической плоскости задавалась начальная численность х0,у0, принадлежащая бассейну несинфазных режимов (если он существует). Для простоты начальная точка выбиралась на отрезке прямой у = 2Н — х (х £ [0; 2Н]) таким образом, чтобы величина кх была больше наперед заданного малого е (чтобы исключить бассейны с частичной синхронизацией).
Для большей наглядности на рис. 5 приведены бифуркационные диаграммы предельных режимов системы (5), полученные при движении одного из параметров по противофазным листам области существования мультистабильных режимов. Наблюдаемые разрывы и скачки предельных состояний на этих диаграммах соответствуют смене одного режима синхронизации на другой. Например, первый разрыв в области существования 2-цикла при движении в сторону роста репродуктивного потенциала а, при фиксированном т, соответствует пересечению границы устойчивости противофазных циклов длины 2, для которой А* = +1, и появлению в спектре (8) и (9) решений, соответствующих устойчивым несинфазным циклам длины 2, несимметричных к полностью синфазному циклу той же длины систем (3) и (5). Далее, после пересечения границы |А* | = 1 динамический режим определяется областью, обозначенной Q и ограниченной двумя противофазными предельными кривыми, сформировавшимися вокруг 21-цикла. Выход из этой области после пересечения
Рис. 5. а - Области существования синфазных и несинфазных циклов в фазовой плоскости; б-д -бифуркационные диаграммы смены режимов системы (3), наблюдаемых при изменении параметров и переходах через границы существования несинфазных циклов
X* = —1 сопровождается рождением трех типов противофазных циклов длины 4 (41, 42, 4з) или синфазного 4о цикла (рис. 5, д).
Особенности численного алгоритма нахождения бассейнов противофазных режимов и их сильно раздробленный характер при больших а, а также отсутствие механизма точной идентификации противофазных режимов приводят к тому, что движение по противофазным листам сопровождается частыми «перескоками» с одного листа синхронизации на другой. Подобные разрывы можно наблюдать при движении в направлении изменения коэффициента связи т. Так, с уменьшением т при пересечении границы, соответствующей X* = +1, устойчивый 21-цикл формируется из хаотической динамики (из области, обозначенной С на рис. 5, а-б) или из синфазного 8-цикла (рис. 5, в).
Рассмотрим теперь некоторые особенности формирования циклов 2о и 21. На рис. 6 показано взаимное расположение бассейнов притяжений, кривых (8) и траектории системы (5). Серым цветом обозначен бассейн синфазных синхронных режимов, белым - противофазных. После второй бифуркации в параметрической подобласти Н > 2(1 — т)/(1 — 2т) (при т < 0.5) устанавливается один из режимов синхронизации: полностью синфазный или противофазный 2-цикл. Прежде чем достигнуть одного из этих возможных устойчивых циклических состояний, траектория
Рис. 6. а-г - Фазовая плоскость, фазовые кривые и бассейны притяжений системы (5) при а = 8.0, т = 0.017294; достижение синфазного (а, в) и противофазного (б, г) режимов. Динамика (5) при х0 = 2.6, у0 = 1.9 до неустойчивой точки (д) и вблизи неустойчивой точки (г)
системы может проходить в окрестности неустойчивой точки (состояния с частичной синхронизацией), причем часть траектории до этой точки, как и сама эта точка, лежит на границе двух бассейнов синхронизации.
На рис. 6 показаны два варианта достижения синхронного и противофазного режима синхронизации при движении траектории вблизи неустойчивой пары точек. При таком расположении траектории, вернее, вследствие близости начального состояния к границе бассейна, процесс синхронизации оказывается длительным. В окрестности неустойчивой точки происходит выравнивание фаз колебаний двух переменных и один из осцилляторов на три итерации меняет характер своей динамики с периодического на монотонный, что отражается на временных реализациях, показанных на рис. 6, д, е.
Детальный анализ областей нерегулярной динамики, отмеченных на рис. 5 (С и Q), показал существование в них сложно организованных подобластей циклов конечной длины - окон периодичности (рис. 7, а,б): несинфазных, представленных вытянутыми вертикальными полосами сужающимися книзу («ласточки») в области С и Q, и синфазных, представленных горизонтальными полосами разной высоты («языки») в области С. На рис. 7 цифрами с индексом ноль отмечены синфазные циклы соответствующей длины, для несинфазных режимов индексы опущены. На рис. 7, б, приведен в увеличенном масштабе участок параметрической плоскости, на котором детализированы окна периодичности в области квазипериодической динамики.
Для более наглядного представления о выделенных областях приведены бассейны притяжений аттракторов (рис. 7, в-з). Области чисто хаотической динамики соответствует бассейн, содержащий нормально распределенный случайный набор точечных областей, по форме напоминающий белый шум (см. рис. 7, г). В области квазипериодической динамики зафиксирована градиентная структура бассейнов вдоль прямой у = 2Н — х (см. рис. 7, ж). Среди бассейнов циклов конечной длины выделяется бассейн, соответствующий области, в которой реализуются только
Рис. 7. а - Карта динамических режимов несинфазных режимов системы (3) (опущены индексы для несинфазных циклов), б - увеличенный фрагмент параметрической плоскости и в-з - бассейны притяжения этих режимов
синфазные 3-, 6-, 12-циклы, представляющий собой черный квадрат с редкими точками, для которых показатель синхронизации кх из-за ошибок округления случайно отличен от нуля (см. рис. 7, в). Причем данные синфазные циклы существуют при любом коэффициенте связи т, за исключением, может быть, области Q, где при тех же значениях репродуктивного потенциала, помимо всего, существуют несинфазные циклы конечной длины (см. рис. 7, е-з) со сложно устроенными бассейнами. Исключением также может быть случай очень малого т (т < 0.01), когда фиксируются несинфазные 3-, 6-, 12-циклы (см. рис. 7, д).
Заключение
Итак, в ходе данного исследования было показано, что система симметрично связанных популяций, динамика каждой из которых локально описывается логистическим отображением Рикера, обладает сложно устроенной иерархией мультиста-бильных состояний, особенно при малом коэффициенте связи. Для количественной идентификации мультистабильных состояний предложен показатель синхронизации, на основе которого на фазовой плоскости локализованы бассейны притяжений некоторых характерных аттракторов.
Вычисляемый показатель имеет различные значения для циклов длины N изменения переменных в которых различаются лишь начальной фазой. Это позволяет отличать друг от друга не только N их видов, различающихся временным сдвигом динамических переменных относительно друг друга, но плюс еще N их разных фаз, получаемых одновременным сдвигом по времени от начальной точки. Интересно заметить, что фиксированное число временных реализаций, по которым вычисляется этот показатель, позволяет успешно отличать друг от друга синфазные циклы разных фаз.
На основе предложенного количественного показателя синхронизации колебаний удалось локализовать бассейны притяжений на фазовой плоскости и выделить листы мультистабильных состояний. Было показано, что вследствие несимметричности отображения Рикера, переход с одного листа синхронизации на другой сопровождается значительной перестройкой как характера предельной динамики, так и предшествующей ей траектории.
Несимметричный характер модели Рикера с введенной симметричной миграцией приводит к тому, что противофазной синхронизации циклов длины 2 предшествуют три бифуркации. Первая соответствует потери устойчивости ненулевой неподвижной точки и появлению устойчивого синфазного 2-цикла. При второй наблюдается «рождение» вокруг неустойчивой неподвижной точки неустойчивого противофазного 2-цикла. Третья соответствует появлению вокруг неустойчивого противофазного 2-цикла двух пар неустойчивых циклов с частичной синхронизацией и переходу противофазного цикла в устойчивое состояние. Появление последних двух неустойчивых режимов соответствует первому дроблению бассейна притяжений синфазных состояний.
Показано, что при движении траектории вблизи границ бассейнов притяжений процесс синхронизации, оказываясь достаточно длительным во времени, сопровождается эффектами перемежаемости и временной синхронизацией колебаний. На границах бассейнов притяжений 2-циклов при установлении синфазного или противофазного режима происходит выравнивание фаз колебаний, а траектория на фазовой плоскости резко меняет свое направление.
Детальный анализ динамических режимов на параметрической плоскости показал, что помимо ветвящейся структуры листов мультистабильных состояний -явления, типичного для подобных систем, - области нерегулярной динамики характеризуются наличием «окон периодичности», в которых возможно появление как синфазных, так и несинфазных циклов, бассейны которых имеют разнообразную форму.
Описанные бассейны притяжений имеют зачастую сильно раздробленный и даже фрактальный характер. Это указывает на то, что пространственная синхронизация в большой степени оказывается достаточно случайным событием. Малые флуктуации популяционных параметров или долей мигрантов могут быстро изменить форму бассейнов и вызвать перескок популяционной динамики с одного листа синхронизации на другой.
Исследования проведены при финансовой поддержке РФФИ (проект №09-04-00146-а) и конкурсных проектов ДВО РАН (проекты №209-^15-01, №09-I-ОБН-12).
Библиографический список
1. Свирижев Ю.М.Нелинейные волны, диссипативные структуры и катастрофы в экологии. М.: Наука, 1987. 368 с.
2. Кисинг Ф., Остфельд Р. Связи популяционной динамики и расселения мелких млекопитающих с изменением сообществ в мозаичной среде: современное состояние и перспективы // Сибирский экологический журнал. 1999. № 1. С. 15.
3. Kaneko K. Clustering, coding, switching, hierarchical, ordering, and control in network of chaotic elements // Physica D. 1990. Vol. 41, № 2. P. 137.
4. Фрисман Е.Я. О механизме сохранения неравномерности в пространственном распределении особей // Математическое моделирование в экологии. 1979. С. 145.
5. Frisman E.Ya. Differences in densities of individuals in populations with uniform ranges // Ecological Modelling. 1980. Vol. 8. P. 345.
6. Фрисман Е.Я. Первичная генетическая дивергенция (теоретический анализ и моделирование). Владивосток: ИАПУ, ДВНЦ АН СССР, 1986. 160 с.
7. Кузнецов С.П. Универсальность и подобие связанных систем Фейгенбаума // Изв. вузов. Радиофизика. 1985. Т. 27, № 8. С. 991.
8. Кузнецов С.П.Переход от симметричного к несимметричному режиму хаотической динамики в системе диссипативно связанных рекуррентных отображений // Изв. вузов. Радиофизика. 1989. Т. 32. № 1. С. 49.
9. de Castro Manuela L., Silva Jacques A.L., Justo Dagoberto A.R. Stability in an age-structured metapopulation model // J. Math. Biol. 2006. № 52. P. 183.
10. Silva Jacques A.L., Giordani Flavia T. Density-dependent migration and synchronism in metapulsations // Bulletin of Mathematical Biology. 2006. № 68. P. 451.
11. Wysham Derin B., Hastings Alan. Sudden shifts in ecological systems: Intermittency and transients in the coupled Ricker population model // Bulletin of Mathematical Biology. 2008. № 70. P. 1013.
12. Omelchenko I., Maistrenko Y., Mosekilde E. Synchronization in ensembles of coupled maps with a major element // Discrete Dynamics in Nature and Society. 2005. № 3. P. 239.
13. Безручко Б.П., Прохоров М.Д. Селезнев Е.П. Виды колебаний, мультистабиль-ность и бассейны притяжений аттракторов симметрично связанных систем с удвоением периода // Изв. вузов. ПНД. 2002. Т. 10. № 10. С. 47.
14. Астахов В.В., Шабунин А.В., Стальмахов П.А. Противофазная синхронизация и формирование мультистабильности в симметрично связанных бистабильных системах // Изв. вузов. ПНД. 2006. Т. 14. № 6. P. 112.
15. Астахов В.В., Шабунин А.В., Стальмахов П.А. Бифуркационные механизмы разрушения противофазной синхронизации хаоса в связанных системах с дискретным временем // Изв. вузов. ПНД. 2006. Т. 14. № 6. С. 100.
16. Шабунин А.В., Демидов В.В., Астахов В.В., Анищенко В.С. Количество информации как мера синхронизации хаоса // Письма в ЖТФ. 2001. Т. 27, вып. 11. С. 78.
17. Kreuz T., Mormann F. Andrzejak R., Kraskov A. Lehnertz K., Grassberger P. Measuring synchronization in coupled model systems: A comparison of different approaches // Physica D. 2007. Vol. 225. P. 29.
18. Смирнов Д.А., Сидак Е.В., Безручко Б.П. Статистические свойства оценки коэффициента фазовой синхронизации // Изв. вузов. ПНД. 2008. Т. 16. № 2. С 111.
Институт комплексного анализа Поступила в редакцию 2.04.2010
региональных проблем ДВО РАН
SYNCHRONIZING THE PERIOD-2 CYCLE IN THE SYSTEM OF SYMMETRICAL COUPLED POPULATIONS WITH STOCK-RECRUITMENT BASED ON THE RICKER POPULATION MODEL
M.P. Kulakov, E.Ya. Frisman
We investigated coupled map lattices based on the Ricker model that describes the spatial dynamics of heterogeneous populations represented by two connected groups of individuals with a migration interaction between them. Bifurcation mechanisms inphase and antiphase synchronization of multistability regimes were considered in such systems. To identify a synchronization mode we introduced the quantitative measure of synchronization.
Keywords: Population dynamics, migration, coupled map lattices, bifurcation, synchronization, multistability, basin of attraction.
Кулаков Матвей Павлович - родился в 1982 году, окончил факультет математики и информатики в Биробиджанском государственном педагогическом институте (ныне Дальневосточная государственная социально-гуманитарная академия). Работает младшим научным сотрудником в лаборатории математического моделирования динамики региональных систем ИКАРП ДВО РАН. Основные научные интересы сосредоточены в области математического моделирования экологических систем, нелинейной динамики, синхронизации, теории колебаний.
679016 Еврейская автономная область, Биробиджан, ул. Шолом-Алейхема, 4 Институт комплексного анализа региональных проблем ДВО РАН E-mail: [email protected]
Фрисман Ефим Яковлевич - родился в 1948 году в Сталинабаде (ныне Душанбе), окончил Новосибирский государственный университет в 1971 году. После окончания НГУ работает в Институте автоматики и процессов управления ДВО РАН. С 1988 года зав. лабораторией математического моделирования экологических систем ИАПУ ДВО РАН. С 2002 года - директор Института комплексного анализа региональных проблем ДВО РАН (Биробиджан). Защитил диссертацию на соискание ученой степени кандидата биологических наук (1982) в Институте цитологии и генетики СО РАН (Новосибирск) и доктора биологических наук (1989) в Институте биофизики ДВО РАН (Красноярск) в области математической теории эволюции. Автор 4-х монографий (две вышли в издательстве «Наука» и две - в ДВО АН СССР). Опубликовал более 120 научных работ по общим вопросам математического моделирования и по конкретным моделям экологических и популяционных систем. Соросовский профессор.
679016 Еврейская автономная область, Биробиджан, ул. Шолом-Алейхема, 4 Институт комплексного анализа региональных проблем ДВО РАН E-mail: [email protected]
