CC BY
68
Изв. вузов «ПНД», т. 17, № 1, 2009
УДК 519.245+517.9+517.39
ВЛИЯНИЕ ПЕРЕМЕШИВАНИЯ И ДИФФУЗИИ НА ПРОСТРАНСТВЕННО-ВРЕМЕННУЮ ДИНАМИКУ В СТОХАСТИЧЕСКОЙ СИСТЕМЕ ЛОТКИ-ВОЛЬТЕРРЫ С ДИСКРЕТНЫМ ФАЗОВЫМ ПРОСТРАНСТВОМ
А.В. Ефимов, А.В. Шабунин
В работе рассматривается влияние диффузии и перемешивания на динамику стохастической системы Лотки-Вольтерры. Моделирование осуществляется с помощью метода Монте-Карло. Показывается, что локальная диффузия сильно изменяет динамику модели, ускоряя процессы взаимодействий на решетке, а перемешивание приводит к появлению глобальных периодических колебаний. Выясняется, что рождение глобальных колебаний происходит благодаря явлению фазовой синхронизации. В работе подробно рассматриваются различные характеристики системы и их зависимости от параметров. Представленные в статье материалы служат основанием для дальнейших исследований, направленных на изучение возможности управления системами данного типа, а также демонстрируют одну из причин видового многообразия и устойчивости динамики популяций в экосистемах.
Ключевые слова: Пространственные структуры, синхронизация, самоорганизация.
Введение
На протяжении многих десятилетий ансамбли взаимодействующих частиц, а в частности, модели типа Лотки-Вольтерры не перестают интересовать ученых всего мира. На то есть несколько причин. Во-первых, это большая универсальность такого рода систем. Класс моделей Лотки-Вольтерры описывает различные процессы, такие как взаимодействия внутри популяций и в химических системах, пожары в лесах и другие явлениями, основу которых составляют взаимодействия автокаталитического характера. Во-вторых, динамика подобных систем оказывается очень богатой на различные нелинейные эффекты, которые сами по себе представляют немалый интерес для современных исследований. Это такие процессы, как синхронизация, образование пространственно-временных структур, индуцированные шумом переходы, различного рода бифуркации и т.д. [1-21].
Еще одной причиной пристального внимания ученых к системам типа Лотки-Вольтерры, особенно в последние годы [22-26], является расхождение между реаль-
57
ной эволюцией популяций и поведением, предсказанным динамическими моделями Лотки-Вольтерры. В рамках классической модели Лотки-Вольтерры в системе «хищник-жертва» существуют консервативные колебания численности видов. Как показали исследования [22], учет флуктуаций приводит к выживанию только одного вида, что соответствует экспериментальным данным, если в среде обитания есть сильное перемешивание. Однако чаще в реальных экологических и химических системах локальность взаимодействий сочетается с отсутствием сильного перемешивания. Это приводит к малым устойчивым колебаниям численностей видов вблизи состояния равновесия, то есть наблюдается межвидовое разнообразие, к которому мы все так привыкли. Кроме того, следует отметить, что такое поведение очень похоже на поведение диссипативных динамических систем.
В последнее время поиску объяснения подобного несоответствия было уделено достаточно большое внимание. Вопросы устойчивости колебаний и процессы самоорганизации в популяционных моделях подробно рассмотрены в работах [22,24-26] и на данный момент причины расхождений в поведении классических динамических моделей типа Лотки-Вольтерры и реальных систем связывают с эффектами конечной численности популяций, наличием флуктуаций среды и миграционными (диффузионными) процессами. Однако, по-прежнему, остается не совсем ясным влияние каждого из вышеперечисленных эффектов на динамику конкретных моделей, а также то, какими изменениями в пространственной организации системы сопровождается добавление того или иного типа диффузии. Например, хорошо известно, что в открытых нелинейных системах типа «реакция-диффузия», находящихся вблизи пространственно-однородного состояния равновесия, диффузия взаимодействующих элементов приводит к образованию стационарных пространственных структур (структур Тьюринга) [27-29]. В классическом случае для развития такого сценария необходимо, чтобы скорости диффузии отдельных взаимодействующих компонентов стали различными. В данной работе рассмотрен случай, когда независимо от типа диффузионного процесса (локальная диффузия или нелокальное перемешивание) все элементы ансамбля частиц диффундируют с одинаковой скоростью. Более того, как было показано в предыдущих работах [7-9], в нашем случае возникновение пространственных структур (нестационарных) обусловлено чисто реакционными механизмами и протекает в отсутствие диффузии.
В качестве рабочей была выбрана известная система класса Лотки-Вольтерры (Lattice Lotka-Volterra в англоязычных изданиях, или сокращенно LLV) [7]. Мы показываем, что процессы перемешивания играют роль связи между отдельными участками пространства системы, в которых происходят несинхронные колебания. Увеличение степени перемешивания приводит к бифуркационному рождению устойчивых глобальных колебаний в системе, благодаря явлению фазовой синхронизации между локальными осцилляторами на поверхности.
Статья построена следующим образом. В первом разделе описывается исследуемая система и методы моделирования ее динамики. Во втором и третьем разделах с помощью разновидности метода Монте-Карло, представляющего собой «виртуальный эксперимент», проводится исследование влияния диффузии и перемешивания на систему. В четвертом разделе проводится исследование явления фазовой синхронизации локальных колебаний в системе. Наконец, в Заключении обсуждаются основные результаты.
58
1. Система LLV. Динамическая модель среднего поля и метод Монте-Карло
Рассмотрим квадратную решетку, содержащую N = L х L идентичных ячеек. Каждая ячейка представляет собой элемент ансамбля, который может находиться в одном из трех состояний (фаз), обозначенных далее по тексту X, Y и S. Каждый элемент ансамбля может совершать циклические вероятностные переходы из одной фазы в другую под воздействием своих ближайших соседей. Если кинетические константы, характеризующие вероятность переходов, обозначить как ki, k2 и k3, то для исследуемой модели можно записать следующую схему превращений:
Первая строка схемы означает, что случайно выбранная ячейка X переходит в состояние Y с вероятностью, пропорциональной коэффициенту ki, если случайно выбранная соседняя ячейка находилась в фазе Y. Если соседняя ячейка находилась в любом другом состоянии, то ячейка X сохраняет свою фазу неизменной. Таким же образом осуществляются остальные переходы, описываемые второй и третьей строкой схемы.
Схема (1) принадлежит к классу вероятностных клеточных автоматов. Данная модель очень проста, однако она позволяет изучать ряд нетривиальных эффектов, возникающих в результате локальных взаимодействий в активных средах различной природы.
С динамической точки зрения подобные модели описываются с помощью уравнений среднего поля (так называемое среднеполевое приближение). Однако на сегодняшний день хорошо известно, что динамика систем, предсказанная в рамках этого приближения, плохо соответствует наблюдаемым на практике результатам. Расхождение результатов обычно связывают с конечным числом элементов в ансамбле, наличием пространственных ограничений в решеточных моделях и, как следствие, локальностью взаимодействий, которые не могут быть учтены методом среднего поля. Однако этот метод остается тем инструментом, с помощью которого можно легко узнать общую информацию о поведении системы в первом приближении.
Динамика системы (1) в приближении среднего поля подробно описана в работе [7], поэтому коротко упомянем лишь те результаты, которые будут необходимы для последующего рассмотрения. Уравнения среднего поля для исследуемой системы записываются в следующем виде:
здесь x и у - динамические переменные, обозначающие относительные концентрации соответствующих состояний, а ki, k2 и кз - кинетические константы фазовых переходов. При выводе уравнений (2) учтено свойство, отражающее постоянство
X + Y 2Y, Y + S ^ 2S, S + X ^ 2X.
(1)
X = -kixy + k3x(1 — x — у) y = kixy — k2y(1 — x — у),
(2)
59
общего числа ячеек, что позволило исключить из рассмотрения относительную концентрацию s: x + у + s = 1. Уравнения (2) относятся к уравнениям типа Лотки-Вольтерры. Переменные имеют физический смысл в области, ограниченной следующими прямыми на фазовой плоскости: x = 0, у = 0, x + у = 1. Поэтому будем рассматривать поведение системы внутри замкнутого контура, образованного этими прямыми.
Система дифференциальных уравнений (2) является консервативной и имеет четыре состояния равновесия: три седловых точки Pi(0; 0), P2(1; 0), P3(0; 1) и одну точку типа центр P4 (k2/(k1 + k2 + k3); k3/(k1 + k2 + k3)). Модель демонстрирует консервативные периодические колебания внутри рассматриваемой области. Фазовый портрет системы представляет собой континуум вложенных замкнутых траекторий вокруг центра P4. Вариации параметров системы ki не приводят к качественным изменениям ее динамики в смысле бифуркаций, но влияют на частоту колебаний и на их форму.
В силу указанных выше недостатков метода среднего поля будем использовать другой метод, основанный на алгоритмах Монте-Карло, для моделирования кинетики ансамблей частиц (Kinetic Monte Carlo, далее по тексту KMC). Данный метод представляет собой разновидность вероятностного клеточного автомата и позволяет провести моделирование процессов на микроуровне. При использовании этого метода вначале осуществляется инициализация решетки в соответствии с начальными условиями. В данной работе рассматриваются статистически пространственнооднородные начальные распределения состояний по решетке. После инициализации производится моделирование динамики решетки в дискретном времени. Каждый временной шаг метода Монте-Карло (Monte-Carlo Step, далее по тксту MCS) состоит из N = L2 микрошагов, где L - линейный размер решетки в ячейках. Каждый микрошаг алгоритма включает в себя три этапа, которые приведены ниже.
1. Случайный выбор ячейки решетки и одной из ее ближайших соседей. Ближайшими соседями будем считать ячейки, имеющие с выбранной смежные стороны.
2. Проверка состояний выбранных ячеек и их сопоставление с вариантами, описанными в исходной схеме (1).
3. Если ячейки удовлетворяют одной из строк схемы, то происходит соответствующий переход с вероятностью ki/max(ki, k2, k3). Если нет, то состояния ячеек остаются неизменными.
Считаем, что после N микрошагов алгоритма каждая ячейка была «посещена» в среднем один раз, и решетка переходит в новое, измененное состояние, а дискретное время увеличивается на At = 1 MCS.
Для получения результатов, представленных в данной работе, в качестве основного мы использовали этот же алгоритм, дополняя его необходимыми нововведениями, моделирующими либо локальную, либо нелокальную диффузию. Суть нововведений заключается в том, что после выполнения N микрошагов (то есть после каждого шага MCS), в алгоритм добавляется эффект тасования содержимого ячеек, соответствующий тому или иному типу диффузии. Эти части алгоритма будут подробно описаны в соответствующих разделах.
Результаты моделирования систем LLV без диффузии методом KMC описаны в статьях [7-9]. Прежде чем перейти к рассмотрению влияния различных типов диффузии на поведение системы, напомним известные результаты.
60
После инициализации решетки в соответствии с начальными условиями на поверхности начинается переходный процесс. В течение этого процесса одноименные состояния группируются во фрактальные структуры - кластеры. После окончания переходного процесса система переходит в стационарный режим. В силу особенностей исходной схемы (1) взаимодействуют только ячейки, принадлежащие границам кластеров. В ходе взаимодействий границы кластеров постоянно перемещаются, сами кластеры меняют форму и размеры, однако структурные изменения в системе не происходят. Если рассматривается решетка малых размеров (L < 30) или параметры ki сильно отличаются друг от друга (например к\ = 0.1, к2 = к3 = 1.0), то велика вероятность того, что система придет к одному из трех седловых состояний равновесия Pi, Р2 или Р3, предсказанных в рамках среднеполевого приближения (случай «отравления» решетки).
Подобная пространственная динамика приводит к соответствующему поведению системы во временном масштабе. В качестве переменных, описывающих состояние поверхности, удобно выбрать усредненные по пространству взаимодействия концентрации состояний ячеек. Обозначим их по аналогии с методом среднего поля как x и у. Такие переменные нельзя называть динамическими в строгом смысле, так как они неоднозначно описывают состояние системы. Одним и тем же значениям x и у (макроусловиям) могут соответствовать различные распределения состояний в пространстве (микроусловия). Однако, как показывают исследования, если разные микросостояния соответствуют одному и тому же макросостоянию, то и эволюция средних концентраций происходит одинаковым образом. То есть из каждой точки xo = x(0), уо = у(0) выходит единственная траектория x(t), y(t), аналогично фазовой траектории динамических систем. Соответственно, x и у ведут себя аналогично динамическим переменным.
Временная реализация системы (рис. 1) начинается с участка, соответствующего переходному процессу. Длительность участка и его форма зависят только от начальных макроусловий и параметров ki и не зависят от начальных микроусловий (в случае их однородности) и линейного размера решетки, если он достаточно велик (L > 30).
После переходного процесса начинается эргодический участок реализации.
Усредненные по пространству взаимодействия концентрации состояний совершают нерегулярные колебания вокруг точки, имеющей в фазовом пространстве (x, у) координаты, близкие к координатам точки Р4 из метода среднего поля. Средний период таких колебаний определяется кинетическими параметрами системы, однако они слабо влияют на интенсивность колебаний. Только в случае сильно отличающихся кинетических параметров, то есть при приближении к порогу «отравления» решетки, наблюдается увеличение дисперсии процесса по сравнению со случаем равенства параметров между собой. Значительно сильнее интенсивность
Рис. 1. Начальные участки временных реализаций x(t) для решеток с линейными размерами L = 128 (■) и L = 1024 (о) ячейки. Начальные концентрации xo = 0.6, у0 = 0.25. Кинетические параметры к1 = 0.7, к2 = к3 = 1.0
61
колебаний зависит от размеров решетки. При увеличении L интенсивность колебаний спадает по закону ox2 ~ L-2. При достаточно больших размерах решетки колебания затухают, и динамика системы представляет собой малые флуктуации переменных вокруг состояния равновесия, расположенного на месте точки P4. Однако, в отличие от среднеполевого приближения, где P4 - центр, состояние равновесия теперь проявляет себя как зашумленный устойчивый фокус.
Таким образом, можно заключить, что динамика системы в рамках модели Монте-Карло сходна с диссипативным поведением динамических систем, а поверхность решетки представляет собой набор пространственно-распределенных независимых либо слабо связанных стохастических осцилляторов.
2. Локальная диффузия в системе LLV
Как уже отмечалось выше, диффузионные процессы играют огромную роль в динамике активных сред. В химических системах диффузия может лимитировать скорости автокаталитических реакций или наоборот, ускорять их. Диффузионные процессы в сочетании с особенностями схемы реакции могут приводить к образованию на поверхности раздела сред автоволновых структур различного типа (реакция Белоусова-Жаботинского и др.) В экологических системах с диффузионными процессами обычно связывают разного рода миграции популяций. На протяжении длительного времени эти процессы тщательно исследуются учеными. Однако связь подобных перемещений популяций с динамикой их численностей не всегда очевидна. Немалую роль играет диффузия и в нанотехнологиях, где процессы образования и разрушения кластеров особенно важны.
Учет даже малой диффузионной составляющей в моделях вышеописанных систем может привести к сильно отличающимся результатам. Кроме того, независимо от природы рассматриваемой системы, фундаментальные аспекты влияния диффузии на динамику могут оказаться полезными для управления данными процессами.
Перед тем как перейти к обсуждению результатов, остановимся на описании алгоритма моделирования системы. Основу алгоритма составляет процедура, уже представленная нами в предыдущем разделе. Дополнением к ней служит алгоритм тасования содержимого ячеек. После каждого шага MCS случайно выбирается ячейка решетки. Далее случайно выбирается одна из ее ближайших соседей и выбранные ячейки обмениваются состояниями. Это приводит к тому, что состояние первой ячейки «перемещается» во вторую. Далее первое состояние «перемещается» в следующую выбранную соседнюю ячейку. Число подобных «перемещений» выбранного состояния определяется округлением случайной величины ^, распределенной по нормальному закону с истинным средним (%) = D и стандартным отклонением 0| = v, до ближайшего неотрицательного целого. Таким образом выбранное состояние совершает случайные блуждания вокруг своего первоначального положения. То же самое происходит с другими ячейками по всей поверхности, что достигается повторением вышеописанной процедуры для N случайно выбранных элементов ансамбля.
Таким образом, интенсивность и характер локальной диффузии на решетке определяется двумя параметрами: D и v. Первый из них представляет собой наибо-
62
лее вероятное число «перемещений», которые совершает содержимое каждой ячейки после очередного шага MCS, а второй - разброс числа «перемещений» относительно D. В данной работе не рассматривается зависимость поведения системы от параметра v, поскольку в ходе предварительных исследований не удалось выявить заметное влияние вариации этого параметра. Мы ограничились малым фиксированным значением v = 0.4, задав таким образом малую дисперсию для числа «перемещений». Как будет показано ниже, намного более заметный эффект на динамику системы оказывает вариация параметра D, отвечающего за интенсивность диффузии.
Рассмотрим влияние локальной диффузии на поведение исходной системы. Для начала обратимся к пространственной динамике. На рис. 2 представлены характерные пространственные распределения состояний ячеек поверхности в начальный момент времени и после переходного процесса для различных параметров диффузии. Состояния ячеек отображаются следующими цветами: X - серый цвет, Y - белый и S - черный. Начальное распределение состояний по поверхности решетки соответствует начальным макроусловиям, носит случайный характер и обладает свойством пространственной однородности (рис. 2, а, б). После переходного процесса (t & 2000 MCS) в отсутствие диффузии на поверхности образуются кластеры (рис. 2, в, г), фрактальная структура которых подробно исследовалась в работе [8]. В силу локального характера процессов на решетке без диффузии, ячейки соседних кластеров могут взаимодействовать только на их границах. Вследствие этого, количество ячеек, способных изменить свое состояние, после образования кластеров уменьшается. Чем больше размеры кластеров, тем меньше ячеек провзаимо-действуют на каждом шаге алгоритма Монте-Карло.
Теперь начнем постепенно увеличивать диффузию в системе. Как видно из следующей пары снимков поверхности (рис. 2, д, е), введение в систему даже слабой
Рис. 2. а - поверхность решетки (512 х 512), б - ее увеличенный участок (128 х 128)при t = 0. Та же решетка после переходного процесса к моменту времени t = 2000 MKS: в, г - в отсутствие диффузии (D = 0); д, е - случай слабой локальной диффузии (D = 2); ж, з - случай сильной диффузии (D = 10). Начальные макроусловия и кинетические параметры для всех случаев: x0 = у0 = 0.1, fci = 0.5, k2 = k3 = 1.0
63
локальной диффузии (D = 2) приводит к росту размеров кластеров и к изменению их пространственной структуры. Границы кластеров как бы «расплываются», становясь более изрезанными, а сами кластеры при приближении к границам теряют свою однородность, оказываясь испещренными ячейками с другими состояниями. Дальнейшее увеличение интенсивности диффузии (рис. 2, ж, з) приводит к усилению вышеописанных эффектов.
Еще одним характерным следствием диффузии является образование спиральных и концентрических автоволновых структур на поверхности решетки. Однако при взаимодействии волны с большим количеством других пространственных структур велика вероятность ее разрушения. Поэтому такие волны, возникая в той или иной части поверхности, при слабой диффузии быстро разрушаются, уступая место новым волновым фронтам и другим пространственным формациям. При увеличении параметра диффузии количество пространственных структур на решетке падает вследствие роста размера каждой из них, что приводит к увеличению времени жизни автоволн. Все вышесказанное соответствует общим представлениями о влиянии локальной диффузии на пространственную динамику рассматриваемых систем.
Изменение пространственной структуры решетки приводит к соответствующим изменениям во временной динамике модели. Рис. 3 демонстрирует участки временных реализаций x(t) после переходного процесса. Они соответствуют решеткам, «снимки» поверхности которых представлены на предыдущем рисунке. Даже без проведения дополнительных исследований видно, что степень регулярности колебаний увеличивается с ростом коэффициента диффузии. Также растут интенсивность и средняя частота колебаний. Кроме того, значение хср, около которого происходят колебания, становится ближе по значению к соответствующей координате точки
Для более детального изучения свойств колебаний рассмотрим спектральные и корреляционные характеристики процесса. На рис. 4 и рис. 5 изображены нормированные автокорреляционные функции и спектральные плотности для процесса x(t) при различных параметрах диффузии. При введении даже слабой диффузии в систему время корреляции заметно возрастает, а ширина спектрального пика уменьшается, что соответствует увеличению регулярности колебаний. Сам спектральный пик смещается в область больших частот, что соответствует уменьшению характерного временного масштаба колебаний, то есть увеличению скорости всех процессов, протекающих на решетке. Кроме того, при достаточно сильной диффузии (D = 20) в спектре колебаний, помимо основного пика (/о — 0.0514 MCS-1), появляется полочка, лежащая в интервале частот от / — 0.04 MCS-1 до
/2 — 0.047 MCS-1. Таким образом, в дополнение к основным колебаниям со сред-
равновесия P4 из динамической модели.
и различной интенсивности локальной диффузии: D = 2 (кривая 2), D = 10 (кривая 3). Размер решетки, начальные условия и остальные параметры те же, что и на рис. 2
64
ним периодом То ~ 20 MCS, в системе появляется новый временной масштаб в интервале от 88 MCS до 227 MCS.
Подобные изменения во временной динамике системы при введении локальной диффузии можно объяснить следующим образом. При увеличении параметра диффузии кластеры теряют свою однородность, а общая длинна границ кластеров резко возрастает (вследствие увеличения их изрезанности). Это, в свою очередь, означает рост числа взаимодействий на каждом шаге MCS, то есть уменьшение характерного временного масштаба системы. С дальнейшим усилением диффузии число взаимодействующих на каждом шаге ячеек должно увеличиваться за счет роста неоднородности кластеров и изрезанности их границ. Но этот эффект может в той или иной степени компенсироваться увеличением размеров пространственных структур, ведущим к снижению общей протяженности их границ. Таким образом, средняя частота колебаний fo может иметь довольно сложную зависимость от величины диффузии. Однако общая тенденция в этой зависимости - рост fo. Так, например, для представленного на рис. 5 диапазона параметра D средняя частота колебаний fo растет и достигает некоторого предельного значения, после чего практически не изменяется.
Образование на поверхности новых автоволновых структур, при значительной диффузии охватывающих большие области решетки, приводит к появлению нового временного масштаба, связанного с их временем жизни. Это время примерно на порядок больше периода То и сильно зависит от конкретной пространственной конфигурации. Поэтому в спектре колебаний полочка на низкой частоте не имеет четких очертаний пика. Пока на решетке существует в более или менее неизменном виде текущая глобальная структура, образованная интерференцией отдельных автоволн, средняя интенсивность колебаний высока. Когда структура начинает разрушаться, на поверхности происходит перестройка к новому режиму, сопровождающаяся наличием большего числа более мелких образований, что приводит к уменьшению средней интенсивности колебаний. Далее процесс формирования новой структуры заканчивается и все повторяется. Вышеописанные явления также ответственны за измене-
Рис. 4. АКФ процесса x(t) при разных значениях диффузии D: (a) - 0.0; (б) - 1.0; (в) - 5.0; (г) - 20.0. L = 256 ячеек, k1 = k2 = k3 = 1.0
Рис. 5. Спектры мощности колебаний x(t), построенные по АКФ, изображенным на рис. 4
65
ние характера спада автокорреляционной функции при последовательном увеличении параметра D (см. рис. 4).
Помимо прочего, образование глобальных структур на решетке и некоторая «перемешанность» состояний являются причиной роста регулярности и интенсивности колебаний. Остановимся более подробно на зависимости интенсивности колебаний от параметра диффузии D и линейных размеров решетки L. На рис. 6 представлены зависимости дисперсии процесса x(t) для различных размеров решетки. Из графиков видно, что в указанном диапазоне дисперсия растет линейно с ростом параметра D. Это означает, что средняя амплитуда колебаний растет как квадратный корень из значения параметра, то есть ох ~ л/D. Угол наклона характеристик на рис. 6 к оси абсцисс определяется линейным размером решетки. Запишем уравнение прямой, аппроксимирующей экспериментальные точки зависимости ох2 (D) для фиксированного значения L в виде
Ox2 (D) = kD + Ох2о, (3)
где ох20 - дисперсия процесса в отсутствие диффузии, а к - угловой коэффициент прямой. Коэффициент к и слагаемое ох2о зависят от L. Что касается второй зависимости, то ее вид качественно совпадает с аналогом для системы LLV большей размерности (см. [9]), а график хорошо аппроксимируется функцией f (x) = axa, где a =const и a =const. График зависимости к = к (L) изображен на рис. 7 в логарифмическом масштабе, и представляет собой прямую. Таким образом, если точки графика на рис. 7 принадлежат прямой f (x) = fix + c (где в и c - некоторые константы), то исходная зависимость имеет следующий вид:
к (L) = ecLe. (4)
Определив в и c из рис. 7, а также используя уже известные результаты для аппроксимации ох2о, на основании (3) и (4) получим итоговое оценочное выражение для значения дисперсии процесса:
ох2 (D, L) = dLeD + aLa, (5)
Рис. 6. Зависимость дисперсии колебаний x(t) от параметра D в решетках различного размера. Параметры ki были выбраны равными единице
Рис. 7. Зависимость угловых коэффициентов прямых, аппроксимирующих экспериментальные точки графиков рис. 6, от линейного размера решетки в логарифмическом масштабе
66
где c1 = ec = 32.83, a = 8.26, а a & в & —2. То есть ox ~ 1/\/N, что свидетельствует о суперпозиции множества несинхронных колебаний на решетке.
Однако зависимости на рис. 6 можно описать формулой (5) только в приведенном на рисунке интервале значений параметра D. При более высокой степени диффузии зависимость ox2 (D) становится нелинейной. При этом на поверхности решетки начинают формироваться пространственные структуры большого размера, охватывающие области порядка 100 х 100 ячеек и более (средний размер кластеров в системе без диффузии составляет примерно 20 х 20 ячеек). Отметим, что в силу случайного фактора размер, форма и тип структур не отличаются постоянством, потому приведенные выше цифры весьма условны. Кроме того, время жизни этих структур мало (в среднем не превышает 100 MCS), а зависимости их размеров от размеров всей решетки выявлено не было. При еще больших значениях диффузии происходит «отравление» решетки и колебания прекращаются. Кроме того, ход зависимости ox2 (D) определяется еще и кинетическими параметрами системы. Здесь стоит отметить, что, чем ближе значения параметров ki друг к другу, тем больше линейный участок вышеупомянутой зависимости и колебания в решетке продолжают существовать при больших значениях D.
3. Исследование влияния перемешивания на динамику решетки
Рассмотрим теперь систему с нелокальной диффузией. Как и при моделировании локальной диффузии, в этом случае дополнительным шагом к основному алгоритму также является обмен состояниями между парами случайно выбранных ячеек. Однако на этот раз выбираются не только соседние ячейки. Выбор любых двух ячеек теперь равновероятен и не зависит от их координат. Кроме того, при данном типе диффузии за один временной шаг выбранное состояние перемещается только один раз (если, конечно, случайный выбор снова не падет на него уже в новой ячейке). Интенсивность такого типа диффузии вводится как параметр p = n/N, где n - число пар ячеек, которые обмениваются состояниями на каждом шаге алгоритма MCS. Таким образом, если p = 0, то эффект нелокальной диффузии полностью отсутствует в системе, а при p = 1.0 в среднем каждая ячейка дважды обменяется своим состоянием с другой за At = 1 MCS.
В силу особенностей данной диффузии, связанных с равновероятным перемещением выбранного состояния в любую другую ячейку на поверхности, будем также называть этот тип диффузии перемешиванием. Модель системы с перемешиванием на первый взгляд может показаться несколько нереалистичной из-за нелокального характера такого типа диффузии. Однако это не совсем так. Например, в популяционных системах к подобному эффекту может привести деятельность человека, переселяющего представителей каких-нибудь малоподвижных видов с одного места на другое. Либо сами особи могут совершать дальние миграции с последующим пребыванием на одном месте (например, сезонные миграции птиц). В системах химической природы роль перемешивающего фактора могут играть сторонние частицы, осуществляющие транспорт реагентов, но сами не принимающие участия в реакциях.
Заметим, что перемешивание непосредственно не меняет значение средних концентраций состояний, однако, как и в случае с локальной диффузией, сильно влияет на пространственно-временные процессы в решетке. Изменения наблюдают-
67
ся как на микро-, так и на макроуровне. Рассмотрим сначала глобальное поведение. Пример временных реализаций x(t) для p = 0 и p = 0.03 приведен на рис. 8. Видно, что, как и в случае с локальной диффузией, даже очень малое перемешивание меняет колебания в системе. Теперь вместо слабых флуктуаций вокруг состояния равновесия наблюдаются почти периодические колебания большой амплитуды. Причем, как показывают исследования, амплитуда колебаний больше не зависит от размера решетки. Форма колебаний напоминает колебания в динамической системе, их средний период изменяется при вариации параметров hi так же, как и в среднеполевом приближении. Однако, если в модели среднего поля колебания являются консервативными, то в нашем случае налицо диссипативное поведение: как и прежде, траектория может быть разделена на стадии переходных и установившихся колебаний. Во время первой стадии траектория, стартуя от начальных условий, переходит на соответствующий «аттрактор» (кривые, обозначенные символами «о» и «Д» на рис. 9).
Шум при моделировании методом Монте-Карло является неотъемлемой частью динамики системы. Но в предположении бесконечного размера решетки за
Рис. 8. Временные реализации x(t) с начальными условиями х0 = 0.33, у0 = 0.33 для системы без перемешивания (пунктирная кривая) и для той же решетки с перемешиванием. Параметры hi = к2 = к3 = 1.0. Размер решетки L = 2048
Рис. 9. «Аттракторы» для разных значений р. Размер решетки L = 1024. Начальные условия и параметры одинаковы для всех случаев и имеют те же значения, что и на рис. 8
Рис. 10. Спектры мощности глобальных колебаний для двух различных размеров решетки L. Значение р = 0.02. Начальные условия и параметры имеют те же значения, что и на рис. 8
Рис. 11. Дисперсия процесса x(t) в зависимости от интенсивности перемешивания р для разных значений к1. Параметры к2 = к3 = 1.0
68
счет лучшего усреднения колебания на «аттракторах» можно считать периодическими. Это подтверждается тенденциями в изменениях спектров колебаний с ростом L (рис. 10). Спектры колебаний имеют характерные пики на кратных частотах, а шумовой пьедестал значительно снижается при увеличении линейного размера решетки всего в два раза.
Рассмотрим теперь динамику системы при изменении параметра p более детально. Будем постепенно увеличивать уровень перемешивания в системе и следить за дисперсией колебаний. Графики такой зависимости отображены на рис. 11, а несколько фазовых портретов приведены на рис. 9. Пока перемешивание остается достаточно малым, значительных изменений в глобальном поведении системы не обнаруживается. Затем, при достижении p определенного порогового значения в системе происходит рождение периодических колебаний. Сразу отметим, что наличие порогового характера в появлении глобальных колебаний и их близость к периодическому процессу - главные отличия во временной динамики системы с перемешиванием по сравнению с системой с локальной диффузией. Дальнейшее увеличение перемешивания влечет за собой рост амплитуды колебаний. Как и в случае локальной диффузии, дисперсия глобальных колебаний растет линейно. Это означает, что амплитуда колебаний растет как квадратный корень из разности между текущим значением p и его критическим значением ркр
Ox(p) ~ yjp - Ркр•
Таким образом, глобальные колебания рождаются так же, как предельный цикл при бифуркации Андронова-Хопфа в динамических системах. Рост амплитуды колебаний приводит к изменению их формы: они становятся релаксационными. При достаточно сильном перемешивании предельный цикл «влипает» в контур из инвариантных многообразий седловых точек равновесия, изображенный на рис. 9 стрелками. Это ведет к прекращению дальнейших колебаний на решетке и ее «отравлению» одним из состояний.
Детальный анализ глобальных колебаний показал, что их амплитуда и пороговое значение параметра ркр не зависят от размеров решетки, если он достаточно велик. Для небольших решеток картина иная. Это можно увидеть из рис. 12, где построена зависимость дисперсии колебаний от p для решеток разных размеров. Для малых решеток эта зависимость не имеет порогового характера: колебания начинаются практически сразу после «включения» перемешивания. Увеличение размера решетки приводит к изменению формы кривой: она становится ближе к кусочно-линейной зависимости, характерной для больших решеток. Для решеток с линейным размером более 300 элементов эта зависимость в дальнейшем почти не меняется.
возникновения глобальных колебаний для разных значений L: 60 (А); 100 (...); 200 (—); 300 (о); 500 (■). Параметры ki = k2 = k3 = 1.0
69
4. Синхронизация локальных колебаний на решетке
Из исследований глобальных колебаний видно, что малые и большие решетки ведут себя по-разному: если для больших решеток наблюдается порог возникновения колебаний, то на малых решетках колебания начинаются непосредственно с нулевого значения р. Чтобы объяснить это различие, можно предположить, что мы имеем дело с процессом синхронизации колебаний на разных частях большой решетки, результатом которого и является возникновение глобальных колебаний. Можно ли этот процесс синхронизации описать в терминах фазовой синхронизации, которая наблюдается для связанных динамических систем? Чтобы ответить на этот вопрос, нужно разделить большую решетку на малые части и рассмотреть колебания x и у на этих подсистемах.
Выберем две небольшие области на решетке, не имеющие общих границ. Обозначим средние концентрации частиц X на них как xi(t) и X2(t). Рис. 13 отображает временные зависимости Xi(t) и X2(t) для двух различных случаев: с перемешиванием и без. При отсутствии перемешивания колебания xi и X2 выглядят несинхронными (рис. 13, а). «Включение» малого перемешивания ведет как к увеличению их амплитуды, так и к фазовой синхронизации между ними (рис. 13, б). Чтобы продемонстрировать явление фазовой синхронизации более наглядно, можно построить временную зависимость разности мгновенных фаз, полученную из этих реализаций. На рис. 14 приведены графики для трех значений р. Кривая, соответствующая случаю р = 0, демонстрирует броуновское движение разности мгновенных фаз, что говорит о полном отсутствии фазовой синхронизации. Если перемешивание присутствует, но значение р ниже порога возникновения глобальных колебаний (р = 0.007), наблюдаются короткие интервалы фазового захвата, перемежающиеся интервалами броуновского движения. Наконец, при переходе через порог (р = 0.02) фазы захватываются и наблюдается полная фазовая синхронизация.
В наличии синхронизации среди локальных осцилляторов в решетке с перемешиванием можно убедиться, обратившись к пространственным «снимкам» системы. На рис. 15 представлены такие «снимки» в момент времени t = 1000 MCS для разных значений параметра р. В отсутствие нелокальной диффузии мы наблюдаем
Рис. 13. Временные реализации локальных колебаний в разных участках решетки размером 64 х 64 ячейки при p = 0 (а) и p = 0.02 (б). Линейный размер всей решетки L = 256 элементов
70
известную кластерную структуру (см. рис. 2). Если решетку условно разделить на локальные осцилляторы, то они будут иметь случайный разброс начальных фаз, а колебания концентраций в них будут асинхронны. По мере увеличения степени перемешивания отдельные части решетки начинают вести себя как одно целое, равномерно окрашиваясь соответствующими оттенками. При сильном перемешивании (рис. 15, б) вся решетка поочередно окрашивается цветами определенных состояний, что соответствует большой амплитуде колебаний. Кластерные структуры в этом случае практически отсутствуют, либо имеют малый размер.
Фазовая синхронизация наблюдается также в системе с расстройкой по параметрам. Введем расстройку по управляющему параметру к1, распределенную случайным образом по решетке
Рис. 14. Разности мгновенных фаз локальных колебаний для разных значений p
к1 = k(0) + mf (i),
где i - номер подсистемы, f представляет собой случайную функцию i с равномерным распределением на интервале (—к(0); 1 — к(0)), а m < 1 - параметр рас-стройки.Как известно из анализа среднего поля,значения ki влияют на собственную частоту колебаний. Поэтому в системе без перемешивания разность фаз почти монотонно меняется со временем (кривая 1 на рис. 16). Это означает, что средние периоды колебаний в разных областях решетки различаются. Данный случай соответствует отсутствию фазовой синхронизации. Перемешивание приводит к явлению частичной фазовой синхронизации в ее классическом смысле, когда горизонталь-
Рис. 15. а - поверхность решетки (512 х 512) и б - ее увеличенный участок (128 х 128) после переходного процесса к моменту времени t = 2000 MKS при различной степени перемешивания в системе p: а, б - 0.008; в, г - 0.02. Начальное состояние поверхности и пространственная динамика решетки без перемешивания приведены на рис. 2. Начальные макроусловия и параметры системы те же, что и на рис. 2
71
ные участки фазового захвата прерываются интервалами «проскальзывания» фазы (кривая 2 на рис. 16). Наконец, при достаточно большом р наступает полная фазовая синхронизация, когда разность фаз остается ограниченной все время.
Стоит отметить тот факт, что при сильном разбросе значений всех параметров ki по подсистемам режим полной синхронизации может и не наступить при рассматриваемых значениях параметра перемешивания. Однако, увеличивая, примерно, на порядок параметр р, можно добиться такого состояния решетки, когда колебания на ней обладают большей, чем без перемешивания, регулярностью, но амплитуда колебаний мала, то есть система находится вблизи равновесия. Поверхность решетки при этом напоминает начальное однородное случайное распределение состояний, что является следствием очень сильного перемешивания, а колебания в подсистемах оказываются почти полностью синхронизованными. Тем не менее, если две подсистемы значительно отличаются друг от друга, но разброс значений параметров в остальных подсистемах невелик, то полная фазовая синхронизация наступает и при малых значения p и сопровождается ростом интенсивности колебаний в подсистемах.
Рис. 16. Эволюция разности мгновенных фаз локальных колебаний на решетке с расстройкой по параметру hi для следующих значений р: 1 - 0.0, 2 - 0.025, 3 - 0.04
Заключение
В работе исследовано влияние локальной и нелокальной диффузии на динамику системы LLV при ее моделировании с помощью метода Монте-Карло. При отсутствии диффузии ячейки, содержащие одноименные состояния, группируются во фрактальные структуры - кластеры. Процессы перехода ячеек из одной фазы в другую протекают только на границах кластеров, а их характерный размер не зависит от размеров решетки. В результате этого колебания усредненных концентраций состояний во времени представляют собой случайный процесс, интенсивность которого спадает с ростом линейного размера решетки по закону о2 ~ 1/L2. Подобная зависимость характерна для ансамблей несвязанных либо слабосвязанных осцилляторов, совершающих несинхронные колебания.
Локальная диффузия увеличивает характерный пространственный масштаб системы и приводит к увеличению изрезанности границ кластеров. Вследствие этого растет число ячеек, состояния которых меняются на каждом шаге алгоритма. Это, в свою очередь, означает ускорение всех процессов и уменьшение среднего временного масштаба колебаний, а также приводит к росту регулярности и интенсивности колебаний. Зависимость дисперсии колебаний от параметра диффузии D носит линейный, возрастающий характер при малых значениях параметра (D < 20) и ста-
72
новится нелинейной при его увеличении. При достижении параметром диффузии некоторого порогового значения, зависящего от размера решетки и кинетических параметров исходной схемы, решетка «отравляется» одним из состояний, а колебания прекращаются. Угол наклона линейного участка зависимости ox2 (D) определяется набором кинетических констант и размером решетки. Спектр колебаний при D > 10, помимо основного, содержит дополнительный пик на меньшей частоте, соответствующей новому временному масштабу. Предположительно новый временной масштаб связан со средним временем жизни глобальной пространственной конфигурации, которая является результатом взаимодействия локальных структур и автоволн на поверхности решетки. Так как вид конкретной пространственной конфигурации, а значит и ее время жизни, определяется случайным фактором, то дополнительный пик в спектре не имеет четких очертаний.
Добавление нелокального перемешивания приводит к бифуркационному рождению глобальных колебаний в системе при превышении параметром перемешивания критического значения ркр. Колебания носят периодический характер в предельном случае решетки бесконечного размера. Как и в случае локальной диффузии, увеличение перемешивания влечет за собой линейный рост дисперсии. В целом, механизм появления колебаний в системе очень похож на бифуркацию Андронова-Хопфа в динамических системах. При дальнейшем росте параметра перемешивания «предельный» цикл, соответствующий колебаниям усредненных концентраций, в пространстве x-y «влипает» в контур, образованный инвариантными многообразиями седловых точек Pi, P2 и P3 из модели среднего поля. После этого решетка «отравляется».
Как показали исследования, появление близких к периодическим колебаний в системе с перемешиванием происходит благодаря явлению фазовой синхронизации между отдельными участками решетки. При последовательном увеличении перемешивания разность мгновенных фаз двух локальных осцилляторов сначала совершает броуновское движение во времени, затем появляются интервалы фазового захвата, чередующиеся с интервалами «проскальзывания» фазы, и, наконец, разность фаз оказывается ограниченной на протяжении всего процесса. Также показано, что фазовая синхронизация под действием перемешивания возможна и в системе с расстройкой по параметрам при условии, что разброс значений параметра неидентичности (в нашем случае ki) по подсистемам не слишком велик.
Таким образом, можно заключить, что учет шума, связанного со случайным характером взаимодействий в среде и конечным числом элементов ансамбля, в моделях Монте-Карло для систем типа LLV нарушает их консервативное поведение, предсказанное динамическим описанием среднего поля. Однако введение в систему перемешивания с заданной интенсивностью может стабилизировать колебания на предельном цикле, параметры которого не зависят от флуктуаций среды, а определяются лишь заранее заданной степенью перемешивания.
Авторы выражают свою признательность профессору Астеро Провате (Институт физической химии, научный центр «Demokritos», Афины, Греция) за проведенные обсуждения и консультации по материалам данной работы.
73
Библиографический список
1. Ziff R.M., Gulari E., Barshad Y. Kinetic phase transitions in
irreversible surface-reaction model // Phys. Rev. Lett. 1986. Vol. 56. P. 2553.
2. Albano E.V. and Marro J. Monte Carlo study of the СО-poisoning dynamics in a model for the catalytic oxidation of СО // J. Chem. Phys. 2000. Vol. 113. P. 10279.
3. Tammaro M. and Evans J.W. Chemical diffosivity and wave
propagation in surface reactions: lattice-gas model mimicking СО-oxidation with high СО-mobility // J. Chem. Phys. 1998. Vol. 108. P. 762.
4. Liu D.J. and Evans J.W. Symmetry-breaking and percolation
transitions in a surface reaction model with superlattice ordering // Phys. Rev. Lett. 2000. Vol. 84. P. 955.
5. De Decker Y., Baras F, Kruse N. and Nicolis G. Modeling the NO + H2 reaction on a Pt field emitter tip: Mean-field analysis and Monte-Carlo simulations // J. Chem. Phys. 2002. Vol. 117. P. 10244.
6. Zhdanov V.P. Surface restructuring, kinetic oscillations and chaos in heterogeneous catalytic reactions // Phys. Rev. E. 1999. Vol. 59. P. 6292.
7. Provata A., Nicolis G. and Baras F. Oscillatory dynamics in low-
dimensional supports: A lattice Lotka-Volterra model // J. Chem. Phys. 1999. Vol. 110. P. 8361.
8. Tsekouras G.A. and Provata A. Fractal properties of the lattice Lotka-Volterra model // Phys. Rev. E. 2002. Vol. 65. art. no 016204.
9. Shabunin A.V., Efimov A.V., Tsekouras G.A. and Provata A. Scalling, cluster dynamics and complex oscillations in a multispecies lattice Lotka-Volterra model // Physica
A. 2005. Vol. 347. P. 117.
10. Monetti R., Rozenfeld A., Albano E. Study of interacting particle systems: The transition to the oscillatory behavior of a prey-predator model // Physica A. 2000. Vol. 283. P. 52.
11. Antal T, Droz M., Lipowski A. and Odor G. Critical behavior of a lattice prey-predator model // Phys. Rev. E. 2001. Vol. 64. P. 036118.
12. Droz M. and Pekalski A. Different strategies of evolution in a predator-prey system // Physica A. 2001. Vol. 298. P. 545.
13. Satulovsky J.E. and Tome T. Spatial instabilities and local oscillations in a lattice gas Lotka-Volterra model // J. Math. Biology. 1997. Vol. 35. P. 344.
14. Spagnolo B., Cirone M., La Barbera A. and De Pasquale F. Noise-induced effects in population dynamics // J. Phys.: Condensed Matter. 2002. Vol. 14. P. 2247.
15. Ertl G. Oscillatory kinetics and spatiotemporal selforganization in reactions at solid surfaces // Science. 1991. Vol. 254. P. 1750.
16. Ertl G., Norton P.R. and Rustig J.Kinetic oscillations in the platinum-catalyzed oxidation of СО // Phys. Rev. Lett. 1982. Vol. 49. P. 177.
17. Voss C. and Kruse N. Chemical wave propagation and rate oscillations during the ^2/Н2 reaction over Pt // Ultramicroscopy. 1998. Vol. 73. P. 211.
18. Theraulaz G., Bonabeau E., Nicolis S.C., Sole R.V., Fourcassie V., Blanco S., Fournier
74
R., Jolly J.L., Fernandez P., Grimal A., Dalle P. and Deneubourg J.L. Spatial patterns in ant colonies // Proceedings of National Academy of Sciences USA. 2002. Vol. 99, № 15. P. 9645.
19. Ben-Jacob E., Shochet O., TenenBaum A., Cohen I., CzirokA. and VicsekT. Generic modelling of cooperative growth patterns in bacterial colonies // Nature. 1994. Vol. 368. P. 46.
20. Deneubourg J.L., Lioni A. and Detrain C. Dynamics of aggregation and emergence of cooperation // Biological Bulletin. 2002. Vol. 202, № 3. P. 262.
21. Saffre F and Deneubourg J.L. Swarming strategies for cooperative species //
J. Theoretical Biology. 2002. Vol. 214, № 3. P. 441.
22. Reichenbach T., Mobilia M. and Frey E. Coexistence versus extinction in the stochastic cyclic Lotka-Volterra model. // Phys. Rev. E. 2006. Vol. 74. P. 051907.
23. Tokita K. Statistical mechanics of relative species abundance // Ecological Informatics. 2006. Vol. 1, № 3. P. 315.
24. Washenberger M.J., Mobilia M. and Tauber U.C. Influence of local carrying capacity restrictions on stochastic predator-prey models // J. Phys.: Condensed Matter. 2007. Vol. 19. P. 065139.
25. Valenti D., Schimansky-Geier L., Sailer X., Spagnolo B. and Iacomi M. Moment equations in a Lotka-Volterra extended system with time correlated noise // Acta Physica Polonica B. 2007. Vol. 38, № 5. P. 1961.
26. Refael A., Schiffer M. and Shnerb N.M. Amplitude-dependent frequency, desynchronization, and stabilization in noisy metapopulation dynamics // Phys. Rev. L. 2007. Vol. 98. P. 098104.
27. Turing A.M. The chemical basis of morphogenesis // Philosophical Transactions of the Royal Society of London. Series B. Biological Sciences. 1952. Vol. 237, № 641.
P. 37.
28. Murray J.D. A pre-pattern formation mechanism for animal coat marking // Journal of Theoretical Biology. 1981. Vol. 88, № 1. P. 161.
29. Murray J.D. and Maini P.K. A new approach то the generation of pattern and form in embryology // Science Progress. 1986. Vol. 70, № 280. Part 4. P. 539.
Поступила в редакцию 23.09.2008
После доработки 14.01.2009
MIXING AND DIFFUSION EFFECT ON SPATIAL-TEMPORAL DYNAMICS IN STOCHASTIC LOTKA-VOLTERRA SYSTEM WITH DISCRETE PHASE SPACE
A.V. Efimov, A.V. Shabunin
The influence of two types of diffusion on dynamics of stochastic lattice Lotka-Volterra model is considered in this work. The simulations were carried out by means of Kinetic Monte-Carlo algorithm. It is shown that the local diffusion considerably changes
75
the dynamics of the model and accelerates the interaction processes on the lattice, while the mixing results in appearance of global periodic oscillations. The global oscillations appear due to phenomenon of phase synchronization. Various characteristics of the system and their dependence on parameters have been considered in detail. Submitted results form the basis for further researches of the control possibilities for systems with competitive dynamics. They also demonstrate one of the plausible reasons of species diversity and stability of population dynamics in ecosystems.
Keywords: Spatial structures, synchronization, self-organization.
Ефимов Антон Викторович - окончил физический факультет СГУ по специальности радиофизика и электроника (2004). Работает инженером кафедры радиофизики и нелинейной динамики физического факультета СГУ, с осени 2005 года - аспирант кафедры. Автор нескольких научных публикаций. Область научных интересов: нелинейная динамика, хаос, образование структур, хаотическая синхронизация, ее количественный анализ, клеточные автоматы. E-mail: efimov@chaos.ssu.runnet.ru
Саратовский государственный университет им. Н.Г. Чернышевского 410012, Саратов, ул. Астраханская, 83
Шабунин Алексей Владимирович - окончил Саратовский государственный университет (1990). Доцент кафедры радиофизики и нелинейной динамики физического факультета СГУ, кандидат физико-математических наук (1998). Научные интересы - нелинейная динамика, теория колебаний, синхронизация и управление хаосом. Автор более 40 научных публикаций.
E-mail: alexey@chaos.ssu.runnet.ru
Саратовский государственный университет им. Н.Г. Чернышевского 410012, Саратов, ул. Астраханская, 83
76
