Мультистабильность и стохастическая возбудимость структурообразования в модели Брюсселятора Е.Д. Екатеринчук Л.Б. Ряшко ek.ekaterinchuk@gmail.com Lev.Ryashko@urfu.ru УрФУ (Екатеринбург) Аннотация В данной работе рассматривается стохастическая пространственно- временная модель Брюсселятора. В детерминированном случае исследована устойчивость пространственно однородного (гомоген- ного) стационарного состояния, экспериментально получены воз- можные пространственно-временные структуры в зоне неустойчи- вости Тьюринга. Обнаружены зоны мультистабильности, в кото- рых пространственно-временные структуры зависят от начальных условий. В стохастическом случае изучена степень неоднородно- сти структуры в зависимости от интенсивности случайных возму- щений. Исследована индуцированная шумом генерация структур в зоне устойчивого гомогенного состояния. 1 Введение Процессы самоорганизации являются удивительными явлениями, примеры которых встречаются в раз- личных сферах науки [1, 2, 3], и их изучение в теории нелинейной динамики считается наиболее активно развивающимся направлением. В 1952 году Алан Тьюринг в работе «Химические основы морфогенеза» [4] опубликовал исследования, посвященные математической теории возникновения структур в изначально однородной системе типа «реакция-диффузия». Согласно Тьюрингу, в системе под влиянием диффузии может возникнуть неустойчивость, разрушающая исходное гомогенное состояние и приводящая к новой, периодической в пространстве и стационарной во времени диссипативной структуре [5]. Теория реакционно-диффузионных систем берет свое начало в работах [6, 7]. Одной из наиболее простых моделей, описывающих процессы самоорганизации, является предложенная в 1968 году И. Пригожиным и Р. Лефевром модель [8], получившая название «Брюсселятор». В данной работе для этой модели исследуются механизмы структурообразования и стохастической возбудимости. Раздел 2 посвящен исследованию детерминированной модели. В секции 2.1 проводит- ся анализ устойчивости стационарного гомогенного состояния детерминированной модели, представле- на бифуркационная диаграмма. Секция 2.2 посвящена численному моделированию и обзору возможных пространственно-временных структур. В области неустойчивости Тьюринга обнаружены зоны, в которых тип пространственно-неоднородной структуры зависит от начальных данных (явление мультистабильно- сти). В разделе 3 изучаются индуцированные шумом структорообразования в зоне устойчивого гомогенного стационарного состояния. В качестве количественной характеристики при описании процесса структуро- Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org 237 образования используется дисперсия. Исследуется степень неоднородности структуры в зависимости от интенсивности случайных возмущений. 2 Исследование детерминированной модели Пространственно-распределенная модель Брюсселятора задается системой дифференциальных уравне- ний в частных производных: ∂u = a − (b + 1)u + u2 v + Du ∇2 u ∂t (1) ∂v = bu − u2 v + Dv ∇2 v, ∂t 2 2 ∂ ∂ ∂u где ∇2 = ∂x 2 + ∂y 2 , (x, y) ∈ Ω = [0, L] × [0, L]. Краевым условием является непроницаемость границ: ∂n = 0, ∂v ∂n = 0, где n – вектор нормали к границе ∂Ω. Параметры a, b – постоянные концентрации поступающих реагентов. Вследствие своей простоты, модель часто становится предметом аналитических и численных исследований. Пространственно-временная динамика данной модели изучалась в работах [9, 10, 11]. 2.1 Анализ устойчивости стационарного гомогенного состояния Как и в случае точечных моделей, первым этапом изучения детерминированной распределенной систе- мы является исследование устойчивости ее однородного стационарного состояния. Пространственно од- нородным (гомогенным) стационарным состоянием называется состояние системы, при котором значения переменных не зависят от времени и одинаковы в каждой точке пространства. В нашей системе, гомо- генным стационарным состоянием является: (a, ab ). Пространственно однородное стационарное состояние является устойчивым, если малое возмущение, действующее на систему (в том числе и распределенное в пространстве), вызывает малое отклонение ее решения. Исследование устойчивости будет проводиться на основе анализа линеаризованной системы уравнений (1). Пусть u1 и v1 – малые отклонения от пространственно однородных решений u и v, u1 = Aept ei(kx x+ky y) , v1 = Bept ei(kx x+ky y) . Множитель ept характеризует поведение отклонения от стационарного состояния во времени. Множитель ei(kx x+ky y) характеризует отклонение величин переменных от однородного стацио- нарного состояния в точке с координатой (x, y) для собственных функций, соответствующих волновому числу k [12]. Подставляя u = ū+u1 , v = v̄ +v1 в систему (1) и линеаризуя диффузионные члены с помощью разложения Тейлора вокруг равновесия, мы получаем характеристическое уравнение для линеаризованной версии пространственной модели (Juv − pI)(u1 , v1 )T = 0, (2) где a11 − Du k 2   a12 Juv = . a21 a22 − Dv k 2 q Здесь a11 = b − 1, a12 = a2 , a21 = −b, a22 = −a2 , (kx , ky ) – волновой вектор, k = kx2 + ky2 – волновое число, I – единичная матрица 2 × 2. Запишем дисперсионное уравнение в виде: p2 − tr(Juv )p + det(Juv ) = 0, где tr(Juv ) = tr(J) − (Du + Dv )k 2 , det(Juv ) = k 4 Du Dv − k 2 (a11 Dv + a22 Du ) + det(J), a2   b−1 J= . −b −a2 Если для корней дисперсионного уравнения p1,2 выполняется неравенство Re p1,2 < 0, то стационарное состояние является устойчивым. Условием возникновения диффузионной неустойчивости или неустойчиво- сти Тьюринга является устойчивость по отношению к малым однородным возмущениям (соответствующая точечная система устойчива) и неустойчивость по отношению к малым пространственно неоднородным воз- мущениям [12]. Условие устойчивости точечной системы: tr(J) < 0, det(J) > 0. Условие неустойчивости в распределенной системе tr(Juv ) > 0 или det(Juv ) < 0. При положительных коэффициентах диффузии 238 выражение tr(Juv ) является всегда отрицательным, в силу устойчивости равновесия в точечной системе, следовательно для возникновения неустойчивости необходимо выполнения неравенства det(Juv ) < 0. Рассмотрим функцию det(Juv ) в зависимости от k 2 . Минимум данной функции по k достигается при 2 a22 Du + a11 Dv kcr = . 2Du Dv 2 Dv Из неравенств tr(J) = a11 + a22 < 0, a22 < 0 и kcr > 0 следует что a11 > 0 и Du > 1. Из условия неустойчивости det(Juv ) < 0 получаем неравенство r !2 Du b> 1+a . Dv Из этого неравенства вытекает явное выражение для бифуркационного значения √ ? Dv ( b − 1)2 Du = . a2 a) б) Рис. 1: а) Бифуркационная диаграмма при a = 3, Dv = 10; б) корни дисперсионного уравнения для Du = 4.5 (фиолетовая кривая), Du = 4 (красная кривая), Du = 3 (зеленая кривая) в зависимости от волнового числа k. На рис. 1a представлена бифуркационная диаграмма, на которой изображены бифуркационные кривые Тьюринга (пунктир) и Андронова–Хопфа (сплошная кривая) при a = 3, Dv = 10. В дальнейшем анализе рассматривается сечение b = 9 (точки), при этом бифуркационное значение Du? = 4.44. На рис. 1б показаны графики функций max Re pi (k), где pi (k) – корни дисперсионного уравнения для Du = 4.5 (фиолетовая i кривая), Du = 4 (красная кривая), Du = 3 (зеленая кривая). На бифуркационной диаграмме данным значениям параметра Du соответствуют звездочки того же цвета. При Du = 4.5 график лежит ниже нуля, что указывает на устойчивость однородного стационарного состояния и, как следствие, отсутствие структурообразования. При значении Du = 4 в некотором диапазоне k значения функции становятся положительными. Данные значения волновых чисел отвечают гармоникам, не затухающим с течением времени. При уменьшении коэффициента диффузии (Du = 3), диапазон волновых чисел увеличивается, что приводит к увеличению числа незатухающих гармоник. На процесс структурообразония большое влияние оказывает конкуренция гармоник [10]. Каждому значению q волнового числа k соответствует набор мод, для которых длина волновых векторов равна k = kx2 + ky2 . Гармоники, соответствующие заданному k, с течением времени конкурируя между собой, приводят к новой пространственной структуре, зависящей от конечного числа доминирующих режимов. Распространенный способ определения доминирующей гармоники связан с изучением распространения соответствующего ей волнового фронта с помощью метода амплитудных уравнений. 2.2 Численное моделирование и структуры Тьюринга Численное моделирование системы (1) проводилось с помощью использования явной схемы: двухто- чечный шаблон для производной по времени и трехточечный шаблон для моделирования диффузионной 239 компоненты. В качестве пространственной области был выбран квадрат со стороной L = 60, с простран- ственным шагом ∆x = ∆y = 1 и шагом по времени ∆t = 0.001, для обеспечения устойчивости данной схемы в рассматриваемом диапазоне параметров. В качестве граничного условия задается непроницаемость границ. В качестве начальных условий использовались малые возмущения гомогенного стационарного состо- яния: u(x, y, 0) = ū + εξxy , v(x, y, 0) = v̄ + εηxy , где ε = 0.01, ξxy , ηxy – пространственный гауссовский шум. Типичный пример возможных структур (Тьюринга) представлен на рис. 2. Здесь изображены структуры типа «ячейки-озера» (Du = 4), «лабиринты» (Du = 3.1), «ячейки-острова» (Du = 2), а также переходные структуры типа «ячейки + полосы» (Du = 3.8) и «полосы» (Du = 2.58). В дальнейшем анализе будут представлены структуры только для переменной u. а) б) в) г) д) Рис. 2: Характерный набор структур: а) структура типа «ячейки-озера» при Du = 4; б) переходная струк- тура типа «ячейки + полосы» (Du = 3.8); в) структура типа «лабиринты» при Du = 3.1; г) переходная структура типа «полосы» (Du = 2.58); д) структура типа «ячейки-острова» при Du = 2. Рассмотрим, как будет меняться динамика системы, если в качестве начальных условий взять не за- шумленное однородное стационарное состояние, а структуры, представленные на рис. 3. Анализ мульти- стабильности проводился в диапазоне значений параметра 3.8 6 Du < Du? . При значениях параметра 4 6 Du < Du? , независимо от выбора начальных значений, решения сходятся к структурам типа «ячейки-озера». В зоне 3.8 6 Du < 4, в зависимости от выбора начальных значений, решения сходятся к различным структурам, что указывает на мультистабильность. Продемонстрируем данное явление для фиксированного значения параметра Du = 3.8. На рис. 4 представлены аттракторы, к которым сходятся решения системы, если в качестве начальных условий взять соответствующие струк- туры из рис. 3. На рис. 4 видно, что одновременно сосуществует несколько аттракторов – разнообразных пространственно-временных структур. В этих структурах меняются соотношения «ячеек» и «полос»: в одних преобладают «ячейки» , в других – «полосы». Мульстистабильность может сыграть важную роль в возможных феноменах, вызванных случайными возмущениями. В данной работе рассматривается влияние шума на устойчивое гомогенное стационарное состояние. 240 а) б) в) г) д) Рис. 3: Структуры, принимаемые в качестве начальных условий при численном моделировании. а) б) в) г) д) Рис. 4: Сосуществующие аттракторы системы для значения параметра Du = 3.8. 3 Индуцированные шумом структурообразования Рассмотрим пространственно-распределенную модель Брюсселятора со случайными возмущениями ∂u = a − (b + 1)u + u2 v + Du ∇2 u + ε1 ξ ∂t (3) ∂v = bu − u2 v + Dv ∇2 v + ε2 η, ∂t где ξt , ηt – случайные стандартные гауссовские процессы, величина ε1,2 характеризует интенсивность воз- мущений. В данной работе ε = ε1 = ε2 . Моделирование стохастически возмущенного решения проводилось следующим образом: к значениям 241 а) б) Рис. 5: а) Устойчивое гомогенное стационарное состояние при Du = 4.5, ε = 0; б) индуцированная шумом структура типа «ячейки-озера» для Du = 4.5, ε = 0.005. udet (xi , yi ), vdet (xi , yi ), полученным при моделировании динамики детерминированной системы, добавля- лись случайные возмущения: u(xi , yi ) = udet (xi , yi ) + εξij , v(xi , yi ) = vdet (xi , yi ) + εηij . Напомним, что у детерминированной системы в зоне устойчивого однородного стационарного состояния структурообразование не наблюдается (см. рис. 5a, Du = 4.5). В присутствии случайных возмущений при определенных значениях интенсивности возникают индуцированное шумом структурообразование. На рис. 5б изображен аттрактор стохастической системы при Du = 4.5, ε = 0.005. Как видно из рис. 5б, решение системы сошлось к структуре типа «ячейки-озера». а) б) Рис. 6: а) Дисперсия решений u (фиолетовая линия) и v (синия линия) детерминированной системы (1) в момент времени t = 1000 в зависимости от Du ; б) дисперсия решения стохастической системы (3) в момент времени t = 1000 для трех значений параметров Du = 4.45 (зеленая линия), Du = 4.47 (красная линия), Du = 4.5 (синяя линия) в зависимости от интенсивности шума ε. Для детального количественного описания структур использовалась дисперсия 1 X D= (ui,j − ū)2 , n2 где u – гомогенное равновесие, n – количество узлов сетки в одном направлении. На рис. 6a представлен график дисперсии D решения детерминированной системы (1) в момент вре- мени t = 1000 в зависимости от параметра Du . В зоне устойчивого однородного стационарного состояния (Du > Du? ) дисперсия равна нулю. При уменьшении параметра Du дисперсия начинает увеличиваться, из чего следует, что амплитуда колебаний в структуре типа «ячейки-острова» больше, чем амплитуда в структуре типа «ячейки-озера». В стохастическом случае дисперсия использовалась для количественного анализа индуцированных шу- мом структурообразований. На рис. 6б представлены графики дисперсии в момент времени t = 1000, для трех значений параметров Du = 4.45 (зеленая линия), Du = 4.47 (красная линия), Du = 4.5 (синяя ли- ния). При малой интенсивности случайных возмущений дисперсия близка к нулю и структурообразования не происходит. Когда интенсивность шума достигает критического значения, происходит структурообра- зование и, как следствие, дисперсия резко возрастает и становится D ≈ 0.4. При увеличении параметра 242 системы Du , критическая интенсивность, при которой возникают индуцированные шумом структуры, уве- личивается. Таким образом, критический уровень шума, порождающий структуры, зависит от близости к точке бифуркации. Чем ближе к точке бифуркации, тем меньший шум генерирует структуры. 4 Заключение В данной статье рассматривалась пространственно распределенная модель Брюсселятор. Для данной модели проведен анализ устойчивости стационарного гомогенного состояния, построена бифуркационная диаграмма и представлен обзор типичных структур. В работе исследовано явление мультистабильности. Для значений 4 6 Du < Du? , независимо от начальных значений, решение системы сходится к структуре типа «ячейки-озера». Большое разнообразие структур наблюдается в зоне Du < 4. Для значения Du = 3.8 описан спектр возможных сосуществующих пространственно распределенных структур. В работе изучено явление индуцированного шумом структурообразования в зоне устойчивого стационарного гомогенного состояния. Структурообразование происходит в достаточно узком диапазоне интенсивности шумов, где пространственная дисперсия решений резко возрастает. В работе показано, как такой скачок дисперсии может быть использован для нахождения соответствующего критического значения интенсивности шума. Благодарности Исследование выполнено за счет гранта Российского научного фонда (проект №16-11-10098) Список литературы [1] J. Gollub, J. Langer. Pattern formation in nonequilibrium physics. Rev. Mod. Phys., 71(2):s396–s403, March 1999. [2] M. Wu, G. Ahlers, D. Cannell. Thermally-induced fluctuations below the onset of Rayleigh–Benard convection. Phys. Rev. Lett., 75, 1743 (1995). [3] P. B. Umbanhowar, F. Melo, H. L. Swinney. Localized excitations in a vertically vibrated layer. Nature, 382:793–796, 1996. [4] A. M. Turing. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. B. Biol. Sci., 237:37–72, 1952. [5] P. Glansdorff, I. Prigogine. Thermodynamic Theory of Structure, Stability and Fluctuations. Wiley, New York, 1971. [6] А. N. Kolmogorov, I. G. Petrovskii, N. S. Piskunov. A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem. Chapter in Selected Works of A.N. Kolmogorov: Volume I, editor V.M. Tikhomirov, 242–270, 1991. [7] R. A. Fisher. The Wave of Advantageous Genes. Ann. Eugenics., 7:355–369, 1937. [8] I. Prigogine, R. Lefever. Symmetry Breaking Instabilities in Dissipative Systems. J. Chem. Phys., 48:1695– 1700, 1968. [9] G. Nicolis, I. Prigogine. Self-Organization in Nonequilibrium Systems. Wiley, New York, 1977. [10] A. de Wit. Spatial patterns and spatiotemporal dynamics in chemical systems. Advances in Chemical Physics, 109:435–513, 1999. [11] B. Pena, C. Perez-Garcia. Stability of Turing patterns in the Brusselator model. Phys. Rev. E, 64:056213, 2001. [12] G.Y. Riznichenko. Lectures on mathematical models in biology. RCD, Moscow-Izhevsk, 2010 (in Russian). = Г.Ю. Ризниченко. Лекции по математическим моделям в биологии. РХД, Москва-Ижевск, 2010. 243 Multistability and stochastic excitability of pattern formation in the Brusselator model Ekaterina D. Ekaterinchuk, Lev B. Ryashko Ural Federal University (Yekaterinburg, Russia) Keywords: Brusselator, Turing instability, stochastic excitability, multistability. In this paper, we consider a stochastic spatio-temporal model of Brusselator. In the deterministic case, the stability of spatially homogeneous stationary state is investigated, and spatiotemporal structures in the zone of Turing instability are studied. For the quantitative analysis of these structures, a dispersion is used. Multistability zones in which spatio-temporal patterns depend on initial conditions are found. In the stochastic case, the degree of heterogeneity of the structure depending on the intensity of random perturbations is studied. 244