CC BY
44
Прикладные задачи
^^^^^^^^^^»нелинейной теории колебаний и вслн
Изв. вузов «ПНД», т. 13, № 5-6, 2005 УДК 531.36
МЕТОДЫ ИССЛЕДОВАНИЯ ДЕТЕРМИНИРОВАННОЙ И СТОХАСТИЧЕСКОЙ УСТОЙЧИВОСТИ ГЛИКОЛИТИЧЕСКОГО ОСЦИЛЛЯТОРА
Л.Б. Ряшко, А.В. Смирнов
В работе рассматриваются методы анализа чувствительности циклов к детерминированным и стохастическим возмущениям. В качестве исследуемой модели взята система уравнений Хиггинса, описывающая гликолитические автоколебания. Исследуется возможность описания локальной чувствительности цикла с помощью собственных чисел системы первого приближения. Результаты, наиболее согласующиеся с эмпирическими, были получены при использовании функции стохастической чувствительности. Выявлены наиболее чувствительные к случайным внешним возмущениям участки цикла. Отмечается уменьшение стохастической устойчивости отдельных участков цикла с ростом жесткости системы, в то время как детерминированная устойчивость цикла в целом увеличивается.
Введение
Предельный цикл в системе определяет поведение многих физических, биологических и химических процессов. Примерами таких систем являются модели Ван дер Поля, Тьюринга - Пригожина, реакции Белоусова - Жаботинского (орегонатор), брюсселятора и другие. В биологической кинетике важную роль играют релаксационные автоколебания. В работе рассматривается система дифференциальных уравнений, впервые предложенная Хиггинсом (Higgins) для описания гликолитических автоколебаний. Уравнения такого типа моделируют процессы в первичном жизненном цикле, имеющем немаловажное значение в теории эволюции.
Одни из первых моделей, представленные в работах Хиггинса, давали упрощенное общее описание процессов, имеющих отношение к гликолизу [1-2]. В дальнейшем исследовании процесса гликолиза использовались более сложные модели гликолитических осцилляций, например, модель Селькова -Хиггинса [3]. Изучалось
поведение моделей под влиянием большой внешней периодической силы (что может отражать реальные метаболические процессы в клетке). Вопросам чувствительности гликолитического осциллятора также уделялось внимание в работе Ю.М. Романовского [4]. Объемный обзор моделей колебательных реакций в гликолизе приведен в работе Гарела [5].
Существенный интерес представляет изучение устойчивости предельных циклов (как цикла в целом, так и отдельных его участков) к случайным внешним возмущениям.
Первые работы по исследованию поведения нелинейных осцилляторов под воздействием внешних периодических сил были проведены Ван дер Полем (van der Pol). Нерегулярность разброса возмущенных траекторий вдоль предельных циклов нелинейных систем была рассмотрена Томита (К. Tomita) и Дайдо (Н. Daido) [6], Куррер (С. Kurrer) и Шультен (К. Schulten) [7]. В последней работе авторы исследовали воздействие случайных шумов, действующих вдоль цикла и по нормали к нему на примере модели нейрона, задаваемой уравнением Бонхоффера - Ван дер Поля.
В данной работе рассматривается модель Хиггинса при малых случайных возмущениях. Цель работы - сравнение методов анализа чувствительности циклов к детерминированным и стохастическим возмущениям.
В п. 1 рассматриваются методы исследования детерминированной устойчивости автоколебательных нелинейных систем. Подход, основанный на мультипликаторах, не позволяет исследовать чувствительность отдельных участков цикла. Приводятся способы проведения такого анализа. Рассматривается метод исследования локальной устойчивости детерминированной системы, связанный с переходом в новую систему координат [8].
П. 2 посвящен исследованию стохастической устойчивости предельных циклов при случайных возмущениях. Полное описание поведения возмущенной системы дает уравнение Колмогорова - Фоккера - Планка. Рассматривается аппроксимация плотности распределения с использованием специальной функции - квазипотенциала. Квазипотенциал удовлетворяет уравнению Гамильтона - Якоби, аналитическое решение которого является очень сложной задачей.
В малой окрестности цикла первым приближением квазипотенциала является орбитальная квадратичная форма. Матрица этой квадратичной формы, заданная на детерминированной орбите, дает простое количественное описание поведения нелинейной системы при наличии малых случайных возмущений и позволяет оценить чувствительность отдельных участков цикла на внешние возмущения. Эта матричная функция рассматривается как функция чувствительности цикла.
В п. 3 с помощью мультипликаторов проведен подробный анализ асимптотической устойчивости циклов детерминированной системы Хиггинса. Рассматривается чувствительность соответствующей стохастически возмущенной модели. Методом прямого численного моделирования найдена экспериментальная дисперсия, характеризующая разброс случайных траекторий вдоль цикла. Выявлены чувствительные участки цикла. Отмечается уменьшение стохастической устойчивости отдельных участков цикла с ростом жесткости системы, в то время как детерминированная устойчивость цикла в целом увеличивается.
Исследуется возможность описания локальной чувствительности цикла с помощью собственных чисел системы первого приближения. Метод, использующий
подходящую замену системы координат, позволяет локализовать наиболее чувствительные участки цикла, но лишь приближенно.
Результаты, наиболее согласующиеся с экспериментальной дисперсией, были получены при использовании функции стохастической чувствительности.
В приложении к статье приводится алгоритм вычисления функции чувствительности.
1. Методы исследования детерминированной устойчивости
Математическая модель многих нелинейных автоколебательных детерминированных систем описывается в виде
x = f (x), (1)
где x - n-вектор, f (x) - достаточно гладкая вектор-функция. Система имеет T-периодическое решение x = ~%(t) (на фазовой плоскости - предельный цикл).
Исследование устойчивости периодического решения системы (1) основано на рассмотрении линейной системы первого приближения
z = F(t)z(t), F(t) = f №)), (2)
задающей динамику малых отклонений z(t) = x(t) — ~%(t).
Одной из стандартных характеристик устойчивости цикла являются мультипликаторы pi - собственные числа матрицы монодромии Z(T), где Z(t) - фундаментальная матрица решения системы (2). Аналогом мультипликаторов являются характеристические показатели Ляпунова.
Один из мультипликаторов (для определенности, pi) всегда равен 1. Условие |pi| < 1 для V i = 2, n является необходимым и достаточным для асимптотической орбитальной устойчивости цикла в целом.
Мультипликатор, характеризуя чувствительность цикла к начальным детерминированным возмущениям, является простым описанием потока около цикла за период, но эта интегральная характеристика не дает возможность исследовать устойчивость отдельных участков цикла.
Одну из естественных мер локальной устойчивости цикла дают собственные числа Xi(t) матрицы F(t). Цикл локально устойчив вблизи точки t, если Re(Xi(t)) < 0 для всех i. В противном случае - локально неустойчив.
В целом устойчивый предельный цикл может одновременно иметь участки локальной устойчивости (на которых отклонения траекторий уменьшаются) и локальной неустойчивости (на которых отклонения траекторий увеличиваются). На неустойчивых участках при воздействии внешних возмущений отклонения траекторий от детерминированной орбиты может значительно увеличиваться, даже если в целом (глобально) цикл очень устойчив.
Детальное исследование локальной устойчивости различных нелинейных систем и методы ее количественной оценки описаны в работе Али и Мензингера (Ali F. and Menzinger M.) [8]. Авторы приводят анализ различных сценариев, по которым предельный цикл теряет устойчивость при различных типах бифуркаций. Были рассмотрены примеры, в которых предельный цикл теряет устойчивость в целом ли-
бо «плавно» (в случаях бифуркации Хопфа), либо «резко», когда строго локально устойчивые участки сосуществуют с участками строгой локальной неустойчивости (в случае удвоения периода). Авторами был предложен новый количественный метод оценки локальной устойчивости. Метод заключается в переходе в новую систему координат, базис которой в каждой точке цикла состоит из касательного вектора и векторов нормалей. В новых переменных г = иг, где и(Ь) - матрица ортогонального преобразования, система первого приближения (2) выглядит следующим образом:
|ад = вщ,
где
dU
B(t) = UF(t)U-1 + dfU-1 - матрица устойчивости (stability matrix). (3)
2. Методы исследования стохастической устойчивости
Рассмотрим динамическую систему со случайными внешними возмущениями, которые описываются стохастическим дифференциальным уравнением Ито
X = / (х) + ео(х)и), (4)
где w(t) - винеровский процесс; е - интенсивность возмущений; а(х) - достаточно гладкая п х п - функция, определяющая зависимость помех от состояния системы.
В результате действия случайных возмущений траектории системы сходят с детерминированного цикла и формируют вокруг него некоторый пучок. При воздействии малых шумов большие отклонения от предельного цикла маловероятны. Уравнение Колмогорова - Фоккера - Планка дает наиболее детальное описание вероятностного поведения системы
| = - £ £/ш + е2 £ дХщ(ч(х)р),
где ач = [а(х)аТ(х)}^.
Функция р описывает плотность распределения пучка возмущенных траекторий вокруг детерминированного цикла. Так как уравнения такого типа в общем случае трудно решаются, то асимптотическое поведение характеристик системы при малых возмущениях представляет несомненный интерес.
Первые результаты по исследованию устойчивости нелинейных систем с использованием уравнения Колмогорова - Фоккера - Планка и различных методов преобразования и приближений были получены Рытовым в 1950-х годах [9-10]. Позднее работы в этом направлении велись Дыкманом (М. Бусктап) [11]. Анализ воздействия случайных возмущений на циклы дискретных систем представлен в [12]. Реакция хаотических аттракторов на случайные возмущения исследовалась в [13].
В настоящее время активно разрабатывается новый подход, связанный с использованием специальной функции - квазипотенциала [14]. Квазипотенциал описывает противодействие силам устойчивости, которые влияют на динамику точки в
окрестности аттрактора. Функция квазипотенциала и удовлетворяет уравнению Гамильтона - Якоби
ди\ , 1 / ди , „ т ди\
1 ^ Мда-(х)дгх) = 0 (5)
с условием на аттракторе и = 0, а вне аттрактора и > 0. Можно показать, что асимптотика стационарной плотности р(х, е) при малых е имеет вид
р(х,е) и Ке-и(х)/£2.
Решение уравнения (5) также является сложной задачей даже для двумерного случая, поэтому рассматривается аппроксимация квазипотенциала в малой окрестности аттрактора. Положим Д(х) = х — у(х) - отклонение точки х от аттрактора, где у(х) - ближайшая к х точка на цикле. Первым приближением квазипотенциала [15] является функция
ф(х) = (Д(х), Ф(у(х)) Д(х)), и(х) = ф(х) + 0(||Д(х)||3),
1 д2 и
задаваемая функцией Ф(у) = (у), определенной на аттракторе.
2 дх
Введем матричную функцию Ф(у) = Ф+(у) (+ - псевдообращение). В результате асимптотика стационарной плотности распределения р(х, е) при малых шумах в малой окрестности цикла может быть записана в форме нормального распределения
( ) ( ) ^ ( (Д(х), Ф+(у(х))Д(х)) \ р(х,е) и р*(х, е) = Кехр I--^-)
с ковариационной матрицей е2Ф(у). Матрица Ф(у) характеризует реакцию системы (4) вблизи цикла на случайные возмущения.
Для отыскания значений Ф(у) на цикле используется естественная параметризация, задаваемая Т-периодическим решением . Связывая точки цикла с точками интервала [0, Т), получим соотношение Ф(^(г)) = Ш(г). Матрица Ш(г) является единственным решением системы [15]
Ш = г(г)ш + шгт(г) + р(г) в(г)р(г), ш(0) = ш(Т), ш(г)г(г) = о, (6)
состоящей из уравнения Ляпунова, условия периодичности и условия вырожденности. Здесь
г(г) = д1 (кг)), Б(г) = в(т°т№)), г(г) = 1Ш)), Р(г) = РКФ
где Рг = I — ггт/гт г - матрица проектирования на подпространство, ортогональное вектору г = 0.
Матрица Ш (г) играет роль стохастической функции чувствительности цикла. Ее анализ и методы отыскания рассмотрены в работе [16].
3. Двумерный случай. Система уравнений Хиггинса
Рассмотрим Хиггинса [4]
двумерную систему дифференциальных уравнений
¿х/дг = / (х,у) = 1 - ху, ¿у/¿г = д(х, у) = Ру (х -
1 + Я
(7)
Я + у,
где х - число свободных, у - число присоединенных нуклеотидов к биполимеру, находящемуся в процессе редупликации; Р и Я - параметры. В данной работе рассматривается классический случай Я = 1, при этом параметр Р определяет жесткость системы (наличие ярко выраженных быстрых и медленных участков цикла).
На рис. 1, 2 изображена эволюция предельного цикла на фазовой плоскости при различных значениях параметра системы.
Рис. 1. Эволюция предельного цикла системы Хиггинса при квазигармоническом режиме колебаний (^ = 1) для значений Р: 2.01 (наименьший размер); 2.1; 2.2; 2.3; 2.5. Точка покоя (1,1)
Рис. 2. Эволюция предельного цикла системы Хиггинса при квазигармоническом режиме колебаний (^ = 1) для значений Р: 2.2 (наименьший размер); 3; 4; 5. Точка покоя (1,1)
Из рисунков видно, что с увеличением Р (при фиксированном Я = 1) форма предельного цикла достаточно быстро и качественно меняется: от формы, близкой к эллипсу, к форме прямоугольного треугольника с катетами, проходящими близко от осей координат. С увеличением параметра жесткости системы «катеты» треугольного предельного цикла все более прижимаются к осям координат, и тем больше он растет вдоль оси ОУ. Скорость прохождения «гипотенузы» предельного цикла в несколько раз больше скорости прохождения траекторией других участков цикла.
Мультипликаторы. Для двумерной системы (7) один из мультипликаторов тождественно равен 1, а другой имеет вид
р = ехр
гТ
% <*)>+1«-»
¿т
(8)
где х = - Т-периодическое решение.
Мультипликатор р определяет степень устойчивости цикла к возмущению начальных данных. Значение р указывает (асимптотически), во сколько раз траектория приближается к циклу за один оборот.
0
Условием устойчивости цикла является неравенство р < 1 (критерий Пуанкаре).
График р(Р) - зависимости мультипликатора р от параметра жесткости системы Р - представлен на рис. 3. Из рисунка видно, что при увеличении параметра жесткости системы значение мультипликатора уменьшается, следовательно, чем система жестче, тем цикл устойчивее в целом.
Экспериментальная дисперсия. Для экспериментального исследования влияния внешних возмущений рассмотрим стохастическую систему с аддитивными шумами
Рис. 3. Зависимость мультипликатора от параметра жесткости системы
dx/dt = 1 — xy + ew \,
dy/dt = Py (x — 1+ Q ) + ew2, V Q + У
(9)
где wi (£) (г = 1, 2) - независимые стандартные винеровские процессы, г - интенсивность возмущений.
Добавление шумов приводит к стохастической деформации детерминированной орбиты. Под воздействием случайных возмущений траектории системы сходят с детерминированного цикла и формируют вокруг него некоторый пучок. Для получения оценки экспериментальной дисперсии этого пучка ТТ проведены прямое численное моделирование случайных траекторий и подсчет их отклонений по вектору нормали от детерминированной орбиты в каждой точке цикла.
На рис. 4, 5 показана экспериментальная дисперсия Т(^)) в точках цикла при Р = 2.2. На рис. 4 видна высокая степень неравномерности распределения экспериментальной дисперсии на цикле. Отрезок времени [0, Т] на рис. 4 отражен в долях единицы. Значения дисперсии удобно связать с точками цикла в фазовой плос-
Рис. 4. Экспериментальная дисперсия за период; Рис. 5. Экспериментальная дисперсия в фазовой P = 2.2, Q = 1, т = t/T плоскости; P = 2.2, Q = 1
кости; при этом для наглядности сопоставим каждому значению дисперсии в точке круг пропорциональный значению дисперсии. Рис. 5 позволяет определить наиболее чувствительные участки орбиты к случайным внешним возмущениям: наиболее чувствительными участками цикла являются верхний «угол» треугольного предельного цикла и второй по величине нижний «угол». Рис. 5 показывает, что разброс случайных траекторий имеет нерегулярный характер.
Для сравнительного анализа удобно рассматривать функцию максимума экспериментальной дисперсии в зависимости от параметра жесткости Р
36.0
24.0
12.0
H (P )= тахрШ)}.
2.0 3.0 4.0 Р
Рис. 6. Зависимость максимума экспериментальной дисперсии от жесткости системы
На рис. 6 представлен график ма-симума экспериментальной дисперсии по отношению к жесткости. По графику видно, что с увеличением Р значения фикции Н(Р) растут, что говорит об уменьшении стохастической устойчивости отдельных участков цикла. С ростом Р все больше локализуется участок цикла, на котором достигается максимум экспериментальной дисперсии.
При сопоставлении экспериментальной дисперсии (см. рис. 6) с мультипликатором (см. рис. 3) важно отметить, что с ростом жесткости системы глобальная устойчивость цикла увеличивается (значение р(Р) уменьшается), в то время как стохастическая устойчивость отдельных участков уменьшается (значение Н(Р) увеличивается).
Собственные числа системы первого приближения. Проведем анализ локальной устойчивости циклов системы с помощью собственных чисел матрицы Р(Ь) системы первого приближения (2). Динамика поведения значений собственных чисел матрицы Р(Ь) за период для уравнений Хиггинса показана на рис. 7. Рассмотрим изменение максимума значений собственных чисел на цикле в зависимости от жесткости системы
Л(Р) = тах тахШёкШ),КеХ2т}
ге[0,т)
Из табл. 1 видно, что с ростом параметра жесткости Р растет значение Л(Р),
Рис. 7. Собственные числа Xi, X2 матрицы F(t) за период; P = 2.2, Q = 1, т = t/T
Таблица 1
Максимум собственных чисел матрицы Р(Ь), Q = 1
Параметр P Значение Л(Р)
2.2 0.249 х 101
3.0 0.713 х 101
4.0 1.397 х 101
5.0 2.161 х 101
7.0 3.959 х 101
что характеризует общую тенденцию увеличения чувствительности отдельных участков цикла.
Сравнивая рис. 7 с рис. 4, отметим, что 'кг(Ь), ^(Ь) не позволяют выявить чувствительные участки предельного цикла.
Таким образом, несмотря на то, что собственные числа ^ (Ь) системы (2) и отражают рост локальной неустойчивости, они не характеризуют особенности на цикле и не локализуют участки локальной устойчивости или неустойчивости.
Собственные числа матрицы устойчивости. Рассмотрим для двумерного случая метод исследования локальной устойчивости [8], связанный с базисом, состоящим из касательного и нормального векторов к циклу. Матрица устойчивости (3) в этом случае будет иметь вид
(^-а Ъьп
"Ъ'пЬ ^пп
/2/х + д2ду + /д(/у + дх) (д2 - /2)(/у + дх) + 2/д/ + ду)
V/2 + 92\ 0 f 2gy + g2fx - fg(fy + gx)
где fx, fy, gx, gy - соответствующие частные производные, индекс n - обозначение компоненты нормали, индекс t - компоненты касательной.
Описанный в работе [8] метод (см. п. 1 настоящей статьи) разделяет движение по отдельным направляющим векторам, что позволяет анализировать интересующее нас отклонение траекторий по нормали без влияния касательной составляющей. В этом случае, динамика нормальных отклонений целиком определяется собственным числом \nn. Локальная устойчивость характеризуется знаком Xnn:
• на тех участках, где \nn ^ 0, влияние возмущений мало;
• на тех участках, где \nn ~ 0, траектории движутся вдоль цикла, так как нет превалирующей силы отталкивания или притягивания (обычно при мягком режиме колебаний, например, квазигармонические колебания);
• на тех участках, где \nn ^ 0, воздействие внешних шумов велико и вероятны большие отклонения возмущенных траекторий от детерминированной орбиты.
В работе [8] предполагалось, что амплитуда локальных флуктуации достигает своего максимума именно там, где \nn имеет положительный пик. В этом смысле \nn должна характеризовать меру локальной неустойчивости.
На рис. 8, 9 представлены графики Xnn за период и на фазовой плоскости в точках цикла. При сравнении рисунков экспериментальной дисперсии (см. рис. 5) и \nn (см. рис. 9) в фазовой плоскости действительно обнаруживается некоторая корреляция распределения значений на цикле. Сравним графики экспериментальной дисперсии (см. рис. 4) и A,nn (см. рис. 8) за период. Из графиков видно, что максимум экспериментальной дисперсии достигается на том участке цикла, где 'knn имеет наибольший перепад значений при изменении знака. Таким образом, \nn не описывает в точности распределение экспериментальной дисперсии на цикле. Наиболее наглядно это несоответствие видно из сопоставления рис. 5 и рис. 9.
Рассмотрим максимум собственного числа матрицы стабильности
Лпп(Р) = max {hnn(t) ^.
i Рис. 9. Собственное число матрицы устойчиво-
сти \nn в фазовой плоскости. Серый тон - отри-Рис. 8. Собственное число матрицы устойчиво- цательные значения, черный - положительные; сти Inn за период; P = 2.2, Q = 1, т = t/T P = 2.2, Q = 1
При росте жесткости системы Лпп (Р) растет, что согласуется с результатами прямого моделирования пучка случайных траекторий, приведенных в табл. 2.
Итак, полученные численные результаты для уравнений Хиггинса позволяют сделать вывод о том, что Хпп неявно, но локализует наиболее чувствительные участки цикла. Важно отметить, что динамика распределения значений собственных чисел до замены переменных (г), (£)) и после (компоненты матрицы устойчивости В) качественно различна, при этом компоненты матрицы дают более подходящую оценку меры локальной устойчивости.
Функция стохастической чувствительности. Рассмотрим решение уравнения Ляпунова (6) в двумерном случае. Вследствие вырожденности матрица Ш (г) имеет вид
Ш (г) = [(г)Р (г),
здесь Р(г) = р(г)рт(г), гдер(г)±/(%(г)) и\\р(г)\\ = 1.
Функция [(г) - Т-периодическая скалярная функция чувствительности - удовлетворяет уравнению
[ = а(г)[ + Ь(г) (10)
с Т-периодическими коэффициентами, где
а(г) = рт(Ь)[Гт(г) + ^(г)]р(г) и Ь(г) = рт(г) в(г)р(г).
Функция стохастической чувствительности аппроксимирует (при малых возмущениях и в небольшой окрестности) дисперсию пучка случайных траекторий, формирующегося вокруг детерминированной орбиты [17]. Это простое аналитическое описание позволяет определить наиболее чувствительные к случайным возмущениям участки детерминированной орбиты.
Таблица 2 Значение максимума \пп относительно жесткости системы. Q = 1
Параметр P Значение Лпп (P)
2.2 0.173 х 101
3.0 1.041 х 101
4.0 2.943 х 101
5.0 5.606 х 101
7.0 1.315 х 101
Рис. 10. Функция чувствительности ц за период. Метки - соответствующие значения экспериментальной дисперсии; P = 2.2, Q = 1, Рис. 11. Функция чувствительности ц в фазовой т = t/T плоскости; P = 2.2, Q = 1
На рис. 10, 11 показаны графики функции чувствительности на цикле. Сравнивая рис. 11с рис. 5 можно отметить совпадение наиболее чувствительных участков цикла: верхний «угол» треугольного предельного цикла и, менее выраженный, нижний «угол». Из рис. 10 видно, что поведение функции чувствительности (на рисунке - сплошная линия) на цикле хорошо согласуется с поведением экспериментальной дисперсии (на рисунке - метки).
Уравнения Хиггинса также исследовал Ю.М. Романовский [4]. В его работе отмечены наиболее чувствительные участки цикла к случайным флуктуациям: вторая половина нижнего «катета» и нижний «угол» треугольного предельного цикла.
Рассмотрим максимум функции чувствительности для заданного параметра жесткости Р
М(Р) = тах ЫШ. ге[0,т)1К "
Исследуя систему в квазигармоническом режиме колебаний при увеличении Р вблизи бифуркационного значения параметров Р = 2, отметим некоторое уменьшение значения М(Р) (рис. 12). Это согласуется с поведением мультипликатора. Из графика видно, что при дальнейшем росте жесткости системы Р значение М(Р) начинает возрастать. Стохастическая устойчивость цикла уменьшается [18-19]. Этот важный вывод об уменьшении стохастической устойчивости при росте Р также подтвержден прямым моделированием случайных траекторий.
Рассмотрим связь функции чувствительности и компоненты матрицы устойчивости \пп(£). Эти функции связывает следующее соотношение:
Рис. 12. Зависимость максимума функции чувствительности от жесткости системы
Xn
(t) f+g2 (t) l(t) = —ö— a(t)'
2
а так как а(1) лишь параметр, задающий ц(£) (см. формулу (10)), то *кпп (^ лишь косвенным образом характеризует стохастическую устойчивость.
Таким образом, функция чувствительности позволяет наиболее точно определить и локализовать чувствительные к случайным возмущениям участки детерминированной орбиты, а также отражает динамику роста стохастической неустойчивости отдельных участков цикла с ростом жесткости системы. Отметим также, что получение экспериментальной дисперсии для жесткой системы является трудоемкой вычислительной операцией, что делает использование функции чувствительности более удобным математическим аппаратом для исследования стохастической устойчивости.
Работа поддержана грантом РФФИ № 04-01-96098 Урал.
Приложение
Вычисление функции чувствительности
1. Осуществляется построение детерминированного предельного цикла с заданной точностью для исследуемых параметров системы (Р,Я). Для достижения требуемой точности цикла использовалось сечение Пуанкаре. Очередная точка цикла на каждом шаге вычислялась с помощью метода Рунге - Кутты.
2. Вычисляются значения функции чувствительности [(г) в каждой точке цикла из набора, полученного на предыдущем этапе алгоритма.
Приведем аналитический способ решения уравнения (10). Общее решение однородного уравнения д(г) = а(г)д(г) имеет вид Сд(г), где д(г) = ехр(/* а(т)д,т), С - константа. Метод вариации произвольной постоянной приводит к уравнению
* Ь(т)
с(г) = Ь(г)д(г) с решением с(г) = /
Уо
ю д(т)'
Общее решение имеет вид [(г) = д(г)(с + с(г)). (11)
Используя условие периодичности, найдем константу С
[(Т) = д(Т)(С + с(Т)) = д(0)(С + с(0)) = [(0),
откуда = д(Т )с(Т)
= 1 - д(т) .
Функции д(г) и с(г) можно искать численно, решая задачу Коши
( д(г) = а(г)д(г), I с(г) = Ь(г)/д(г)
(12)
с начальным условием д(0) = 1, с(0) = 0. Используя значения д(г), с(г) и С, находим по формуле (11) искомую функцию чувствительности.
Для системы уравнений Хиггинса со случайными возмущениями (9) функции а(г) и Ь(г) в уравнениях (10), (12) имеют вид
2
а(г) =
/2 + д2
(/х - Р/у - ду) + /2 рх -
21_ РЯ(Я +1)
(Я + у)2
2р2
Ь(г) = д + /2
где f = f (x,y), g = g(x,y) - функции из (7), P и Q - параметры системы.
3. Вычисляется M(P) как максимум значений функции чувствительности, посчитанных в каждой точке цикла для исследуемого параметра P(Q = 1).
Библиографический список
1. HigginsJ. // J.Ind.Eng.Chem. 1967. Vol. 59, № 5. P. 18.
2. HigginsJ. // Proc.N.A.S.(USA). 1964. Vol. 51. P. 989.
3. Sel'kov E.E. // Eur.J.Biochem. 1968. № 4. P. 79.
4. Романовский Ю.М., Степанова Н.В., Чернавский Д.С. Математическое моделирование в биофизике. М.: Наука, 1975. С. 61, 258.
5. Гарел Д., Гарел О. Колебательные химические реакции. М.: Мир, 1986.
6. Tomita К., Daido H. Possibility of chaotic behaviour and multi-basins in forced glycolytic oscillations // Physics Letters A. 1980. Vol. 79, № 2, 3. P. 133.
7. Kurrer C., Schulten K. Effect of noise and perturbations on limit cycle systems // Physica D. 1991. Vol. 50. P. 311.
8. Ali F., Menzinger M. On the local stability of limit cycle // Chaos. 1999. Vol. 9, № 2. P. 348.
9. Рытов С.М.Флуктуации в автоколебательных системах томсоновского типа. I // ЖЭТФ. 1955. Т. 29, вып. 3(9). С. 304.
10. Рытов С.М. Флуктуации в автоколебательных системах томсоновского типа. II //Там же. С. 315.
11. Dyckman M., Chu X., Ross J. Stationary probability distribution near stable limit-cycles far from Hopf bifurcation points // Physical Review E. 1993. Vol. 48, № 3. P. 1646.
12. Кузнецов А.П., Капустина Ю.В. Свойства скейлинга при переходе к хаосу в модельных отображениях с шумом // Известия вузов. Прикладная нелинейная динамика. 2000. Т. 8, № 6. C. 78.
13. Копейкин А.С, Вадивасова Т.Е., Анищенко B.C. Особенности процесса установления вероятностной меры на хаотических аттракторах в системах Лоренца и Ресслера с учетом флуктуаций // Известия вузов. Прикладная нелинейная динамика. 2000. Т. 8, № 6. C. 65.
14. Вентцель А.Д., Фрейдлин М.И. Флуктуации в динамических системах под действием малых случайных возмущений. М.: Наука, 1979. С. 424.
15. Мильштейн Г.Н., Ряшко Л.Б. Первое приближение квазипотенциала в задачах об устойчивости систем со случайными невырожденными возмущениями // Прикл. математика и механика. 1995. Т. 59, вып. 1. С. 51.
16. Bashkirtseva I.A., Ryashko L.B. Sensitivity analysis of the stochastically forced Lorenz model cycles under period-doubling bifurcations // Dynamic Systems and Applications. 2002. Vol. 11, № 2. P. 293.
17. Bashkirtseva I.A., Ryashko L.B. Sensitivity analysis of the stochastically and periodically forced Brusselator // Physica A. 2000. 728. P. 126.
18. Башкирцева И.А., Смирнов А.В., Ряшко Л.Б. Стохастическая устойчивость гли-колитического осциллятора // Материалы III Уральской научно-практической
конференции «Математическое моделирование в медицине и биологии». 2001. С. 18.
19. Башкирцева И.А., Смирнов А.В., Ряшко Л.Б. Сравнительный анализ устойчивости гликолитического осциллятора к детерминированным и случайным возмущениям // Материалы 33-й Региональной молодежной школы-конференции «Проблемы теоретической и прикладной математики». 2002. С. 107.
Уральский государственный Поступила в редакцию 1.11.2003
университет После доработки 3.11.2004
DETERMINISTIC AND STOCHASTIC STABILITY ANALYSIS FOR GLYCOLITIC OSCILLATOR
L.B. Ryashko, A.V. Smirnov
The methods of sensitivity analysis of cycles under deterministic and stochastic disturbances for Higgins model describing glycolytic self-oscillations are considered. Two approaches connected with local exponents and stochastic sensitivity function are compared. The most sensitive parts of cycles are discovered. It was found that some parts of cycle lose stochastic stability along with stability increasing of cycle as whole.
Ряшко Лев Борисович - родился в 1953 году, окончил Уральский госуниверситет (1975), кандидат физико-математических наук, зав. кафедрой вычислительной математики Уральского госуниверситета, Соросовский доцент (1996,1998). Область научных интересов - устойчивость, стабилизация и управление стохастическими системами. Автор ряда статей по данной тематике. E-mail: [email protected]
Смирнов Алексей Валерьевич - родился в 1979 году, окончил Уральский госуниверситет (2001), аспирант кафедры математической физики. Область научных интересов - устойчивость стохастических нелинейных систем.
