CC BY
108
Изв. вузов «ПНД», т. 13, № 5-6, 2005 УДК 530.18
ПРИМЕР РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ТЕОРИИ БИФУРКАЦИЙ В ДИНАМИЧЕСКОЙ СИСТЕМЕ С ШУМОМ
М.-Г.М. Зульпукаров, Г.Г. Малинецкий, А.В. Подлазов
Рассматриваются бифуркации в нелинейных системах, испытывающих воздействие слабого шума. Описаны случаи локальных бифуркаций: «седло-узел», транскритическая бифуркация, суперкритическая «вилка», субкритическая «вилка».
На основании известного явления роста и насыщения уровня шума по мере приближения к точке бифуркации поставлена обратная задача - по наблюдаемому изменению шума (характер роста, уровень насыщения, плотность распределения) определить положение точки предстоящей бифуркации и ее тип. Предложен алгоритм решения обратной задачи.
Введение
В настоящее время построение и анализ математических моделей нелинейных явлений находят все более широкое применение в исследовании поведения сложных систем, в том числе и в областях, традиционно не относящихся к точным наукам -медицине, социологии, экономике, истории и т.д. [1]. Изучение таких нелинейных явлений, как бифуркации, в этих случаях имеет большое прикладное значение, так как они соответствуют кризисным событиям в моделируемых системах.
Замечательным свойством нелинейных систем является то, что и в простых, и в сложных системах возникают бифуркации одного и того же типа. Это дает возможность свести задачу исследования бифуркаций исходной системы к изучению более простой (модельной) системы.
Отметим также следующее важное свойство нелинейных систем. Известно, что во время прохождения точки бифуркации снижается устойчивость системы по отношению к вариации начальных данных [1,2]. Это делает систему чувствительной к малым воздействиям. Подчеркнем, что малые воздействия могут носить случайный характер. Из сказанного следует, что выбор конкретного пути развития в точке бифуркации может осуществиться под влиянием малого управляющего воздействия либо может быть вообще случайным.
Отсюда вытекает наличие двух практических задач.
Во-первых, задача ретроспективного анализа поведения системы с целью определения критических моментов (точек бифуркации) в прошлом, а также исследование альтернативных вариантов развития. Результаты решения этой задачи, помимо чисто теоретического интереса, могут представлять собой опыт, который стоит учитывать в будущем.
Во-вторых, задача определения точек бифуркации в будущем и сосредоточения управляющих усилий в определенном направлении [3]. Исходными данными для задачи, как правило, являются наблюдаемое состояние и поведение системы. Пример решения такой задачи - создание устройств, отключающих авиационные и ракетные двигатели при возникновении в них аномальных шумов, чтобы избежать взрыва и катастрофы.
При рассмотрении поведения сложной системы часто принимается во внимание, что на нее действует случайный шум - малое нерегулярное внешнее воздействие неопределенной природы. Также говорят, что в системе присутствуют слабые флуктуации. Например, если математическая модель описывает систему на каком-то конкретном уровне организации, то для учета влияния других уровней в модель вводятся стохастические составляющие.
Слабые флуктуации не оказывают заметного воздействия на поведение дисси-пативной системы, отличающейся высокой устойчивостью. Если же система находится вблизи точки бифуркации, то следствием ее чувствительности к малым изменениям начальных данных является усиление флуктуаций. Отметим, что усиление флуктуаций может играть конструктивную роль - приводить к образованию упорядоченных структур на макроуровне (порождение так называемых диссипативных структур, или упорядочение через флуктуации [4, 5]).
Ранее было выяснено, что приближение системы к точке бифуркации сопровождается появлением характерных особенностей шума - шумовых предвестников нелинейных неустойчивостей. Например, в [6] описаны предбифуркационные изменения спектра шума для систем дифференциальных уравнений, имеющих периодические решения. Описание предбифуркационного усиления шума в нелинейных отображениях дано в работах [2, 7, 8], в [9] приведен пример экспериментального исследования предбифуркационных шумов; в работе [10] описаны и другие примеры шумовых предвестников. Нас интересует обратная задача - как по нарастанию шума в системе определить, проходит ли она точку бифуркации.
Суммируя сказанное выше, кратко сформулируем цель данной работы как изучение поведения нелинейных систем с шумом вблизи точки бифуркации и поиск возможностей обнаружения точек бифуркации на основе анализа временного ряда, отражающего поведение системы.
1. Состояние проблемы
На сегодняшний день особенно хорошо изучены предбифуркационные шумы в системах с дискретным временем (системах, моделируемых нелинейными отображениями). В работах [2, 7, 8, 11, 12] отмечен ряд следующих явлений.
Явление роста и насыщения уровня шума в окрестности точки бифуркации исследовано в работах [2, 7, 8] на примере нелинейных отображений вида
Хп+1 = F(Хп, ц) + /и , (1)
где ц - параметр, а /п - возмущение (последовательность одинаково распределенных, независимых в совокупности, центрированных случайных величин с дисперсией о?). Отображение имеет устойчивую неподвижную точку Х = Х(ц) (в силу устойчивости математическое ожидание (х) = Х). Рассматривается шум - центрированная случайная величина п = Хп - Х. Бифуркационное условие определяется равенством |у| = 1, где у = ¥'х(х) (например, у = —1 является условием бифуркации удвоения периода), поэтому вводится параметр а = 1 — |у|, имеющий смысл расстояния до точки бифуркации.
Показано, что в линейном приближении дисперсия шума при подходе к точке бифуркации растет как
2 2
(o2) d-f t2 \ _ Of _ Of (2) )iin=iin) _ 1—y2 _ 2a • (2)
Значения параметра, при которых система может быть линеаризована, определяются соотношением
л/3 |e| Of ,„ч
a > amin _-^— > (3)
где e = F'J.x(Х)/2. Вблизи точки бифуркации, когда это условие перестает выполняться, дисперсия шума прекращает расти и выходит на уровень насыщения
(°2) nonlin
Коэффициент усиления возмущения К = о? / о2 при этом достигает максимального
V3 |e| l/°2
значения: Ктах ~ 1/о/.
Важное явление роста и насыщения интервала корреляции шума описано в работах [2, 11]. Показано, что в области применимости линейного приближения интервал корреляции к растет как 1 /а, а вблизи точки бифуркации выходит на уровень насыщения
к ~ 1 2
ктах
amin V3of |e|
В нелинейных системах с непрерывным временем наблюдаются аналогичные явления. Например, в работах [7, 8] рассмотрен частный случай такой системы -нелинейный осциллятор, заданный уравнением
ди
4 + = /®, где в - коэффициент затухания, и - потенциал
и (I) = 1ВI2 +1Л^4,
а шум / (Ь) представляет собой 8-подобный случайный процесс. Интервал корреляции / (Ь) обозначается Tf. Считается, что интервал Tf мал по сравнению с периодом колебаний осциллятора. В данной системе происходит суперкритическая бифуркация типа «вилка» (спонтанное нарушение симметрии) при достижении параметром В критического значения Вс = 0.
Линейное приближение в этом случае даст зависимость
Ни
а/Ч
2|ЗВ
в области
В > Вт1и = а^ Ч
в
а вблизи точки бифуркации дисперсия наблюдаемого сигнала выходит на уровень насыщения
2А = _. ГЧ
\ Ч) поп1т : =
/п0п11п V 12Ав
После прохождения точки бифуркации дисперсия наблюдаемого сигнала некоторое время держится вблизи уровня насыщения, а затем также убывает по линейному закону [2].
Итак, для нелинейных систем с дискретным временем (и ряда частных случаев систем с непрерывным временем) показано, что по мере приближения к точке бифуркации дисперсия наблюдаемого сигнала нарастает строго определенным образом, выходя вблизи точки бифуркации на уровень насыщения. Несколько усложняет наблюдение данного явления то обстоятельство, что интервал корреляции наблюдаемого сигнала, определяющий время, необходимое для сбора статистики, также возрастает в окрестности точки бифуркации.
Естественным продолжением описанных выше исследований было бы решение следующих задач: во-первых, получение аналогичных результатов для как можно более широкого класса нелинейных систем с непрерывным временем; во-вторых, изучение закона распределения наблюдаемого сигнала в зависимости от близости системы к точке бифуркации; и в-третьих, обратная задача - определение положения точки бифуркации и выяснения ее типа по изменению характеристик наблюдаемого сигнала.
2. Постановка задачи
В качестве объекта исследования выберем некоторую диссипативную динамическую систему, поведение которой описывается автономным обыкновенным дифференциальным уравнением с параметром
х = у(х, X), (4)
где х = х(Ь) - переменная состояния, Ь имеет смысл времени, X - параметр, V - некоторая нелинейная функция (фазовая скорость). Также будем пользоваться понятием потенциала
X
и(х, X) = v(z, Х)йг. (5)
о
Для данного типа математических моделей имеется ряд строгих результатов, представляющих собой основу для изучения качественных трансформаций широкого класса динамических систем. В сочетании с относительной простотой, это делает динамическую систему вида (4) особенно удобной для исследования.
Решения уравнения (4) в зависимости от V, X и начального условия х = хо могут вести себя следующим образом [13]: либо стремиться к постоянному значению - аттрактору (фазовое пространство в данном случае одномерно, и единственный допустимый вид аттрактора - устойчивая особая точка х(Ь) ^ х* при Ь ^ ж, v(x*) = 0), либо неограниченно возрастать (убывать) за бесконечное время, либо неограниченно возрастать (убывать) за конечное время.
Соответственно выбору однопараметрической одномерной системы для исследования, мы ограничиваемся рассмотрением наиболее простых бифуркаций - локальных бифуркаций положений равновесия. Рассматриваемые нами бифуркации будут иметь коразмерность 1. Чем ниже коразмерность бифуркации, тем более типичной для разнообразных динамических систем она является. Бифуркационные задачи для систем высокой размерности часто сводятся к рассмотрению систем размерности 1 или 2 [14]. Таким образом, сделанные для системы (4) выводы будут достаточно общими.
Перечислим наиболее распространенные бифуркации коразмерности 1, которые будем рассматривать далее. Это бифуркация типа «седло-узел» с нормальной формой у = ^ — у2 (у - независимая переменная, ^ - параметр), бифуркация типа «обмен устойчивостью» (транскритическая бифуркация) с нормальной формой у = ^у — у2, суперкритическая (мягкая) бифуркация типа «вилка» (у = ^у — у3), субкритическая (жесткая) бифуркация типа «вилка» (у = —цу + у3).
Рассмотрим, как можно в математической модели нелинейной системы учесть влияние шумов. В случае непрерывного времени в качестве математической модели самого шума используют какой-либо случайный процесс, обычно винеровский. Обозначим его Ч(Ь) и назовем вносимым шумом. Математической моделью системы будет стохастическое дифференциальное уравнение. Существует классификация по способу ввода Ч(Ь) в уравнение. Наиболее часто рассматриваются системы с аддитивным шумом - х = и(х, X) + Ч(Ь). Также различают мультипликативный шум -х = v(x, X) х Ч(Ь), внутренние шумы - системы вида х = v(F(х, Ч(Ь)), X), параметрические шумы - системы вида х = и(х, F(X, Ч(Ь))) и т. д.
Мы будем рассматривать броуновский шум, действующий по следующей схеме. Пусть в определенные моменты времени ^ переменная состояния х суммируется с реализацией Ч некоторой случайной величины Н. Назовем это возмущением. Для простоты будем иметь дело только с периодическим шумом (то есть возмущения происходят периодически, с интервалом Тш). Полученная в результате возмущения сумма используется в качестве начального условия в задаче Коши для дифференциального уравнения (4). Система интегрируется в течение времени Тш, затем происходит очередное возмущение, и т. д. В расширенном фазовом пространстве (х, Ь) это выглядит как случайное мгновенное перемещение между интегральными кривыми. В результате получится некоторый случайный процесс х(Ь) марковского типа. Назовем х(Ь) наблюдаемым сигналом. Броуновский шум представляется достаточно простым и удобным для анализа и, в то же время, он не столь хорошо исследован, как, например, аддитивный шум.
Функцию плотности распределения случайной величины Н обозначим /ш и наложим на нее требование четности (снос отсутствует). Стандартное отклонение Н будем обозначать ош.
В дальнейших расчетах будем использовать равномерно распределенный периодический броуновский шум
. №) = / 1/Л, ПК Л/2 , I 0, |£| > Л/2,
Л (6)
ш 2^3'
¿¿+1 — и = Тш, г = 0,1,2,...
Исследуем воздействие шума на систему с параметром X, изменяющимся с течением времени. Пусть изменение параметра происходит в пределах [Хт;п, Хтах] дискретно с постоянным шагом по X, обозначаемым АХ (АХ ^ Хтах — Хтщ), и постоянным интервалом по времени ТХ. Поскольку значительно проще изучать наблюдаемый сигнал как установившийся процесс, нам понадобится, чтобы его интервал корреляции был мал по сравнению с временем изменения параметра. Для этого, как минимум, необходимо выполнение условия
Тш < Тх . (7)
Для простоты в теоретических расчетах полагаем, что параметр изменяется квазистатически: Тх ^ то. Также нам понадобится учитывать еще одно характерное время - время выхода на аттрактор Та. Мы не будем его фиксировать, так как при изменении X время Та меняется в широких пределах, и важно рассмотрение его во всем диапазоне.
Окончательно поставим задачу следующим образом.
Во-первых, в условиях действия периодического равномерно распределенного броуновского шума для случая одной из наиболее типичных бифуркаций требуется получить закон распределения наблюдаемого сигнала и выяснить характер его изменения по мере приближения параметра к точке бифуркации.
Во-вторых, исходя из этого, требуется исследовать возможность решения обратной задачи: на основании изменения характеристик наблюдаемого сигнала в процессе медленного изменения параметра системы оценить расстояние до точки бифуркации и сделать предположения относительно ее типа.
3. Экспериментальные данные
С целью наблюдения явления предбифуркационного нарастания дисперсии наблюдаемого сигнала для перечисленных ранее основных типов бифуркаций был поставлен ряд экспериментов. В этих экспериментах выполнялось численное моделирование поведения динамических систем, заданных дифференциальными уравнениями непосредственно в соответствующих нормальных формах, в присутствии броуновского шума. Для каждой системы была выполнена серия экспериментов и
получена эмпирическая зависимость дисперсии наблюдаемого сигнала от бифуркационного параметра и параметров вносимого шума. В экспериментах рассматривались значения параметров, лежащие в области применимости линейного приближения, так как серийные эксперименты в окрестности точки бифуркации связаны с большими вычислительными затратами, что вызвано упоминавшимся выше предби-фуркационным увеличением интервала корреляции наблюдаемого сигнала.
Эксперименты строились следующим образом. Бифуркационный параметр ц изменялся в пределах от —10 до —1 с шагом 0.02 (точка бифуркации во всех экспериментах ц = 0). Для каждого значения ц выполнялось моделирование воздействия шума на систему. Шаг модели представлял собой интегрирование в течение времени АЬ = 0.001. В начальный момент модельного времени система помещалась в устойчивое положение равновесия у*, вычисляемое для данного ц. В начале шага, если оно совпадало с началом периода возмущения, производилось возмущение. В конце каждого шага положение системы у регистрировалось как реализация случайной величины. По окончании моделирования для зарегистрированных значений у вычислялось стандартное отклонение ау, интерпретируемое как характеристика интенсивности шума. Зависимость ау от ц, полученная в эксперименте, откладывалась на графике. Для каждого значения ц выполнялось 50000 шагов; число шагов подобрано опытным путем таким образом, чтобы быть достаточным для сбора статистики во всех экспериментах для всех значений ц (при уменьшении числа шагов до 1000 получаемые данные начинают заметно меняться от повторения к повторению).
Характеристики вносимого шума в пределах серии экспериментов варьировались по следующей схеме. Как было указано выше, параметрами равномерно распределенного броуновского шума являются амплитуда возмущения А и период возмущения Тш. Величина А принимала значения 10"2, 10"3, 10"4, 10"5 и 10"6, а Тш - 10"3, 10"2 и 10"1. Таким образом, каждая серия включала 15 экспериментов. В процессе моделирования для Тш = 10"3 возмущение производилось на каждом шаге, для Тш = 10"2 - на каждом десятом шаге, а для Тш = 10"1 - на каждом сотом шаге.
Нарис. 1 показаны графики ау (ц) (для разных А и Тш) в двойном логарифмическом масштабе. Видно, что экспериментальные точки аппроксимируются прямыми, то есть зависимость (вне ние наблщцаем°г° отгаага °у области насыщения), как и в случае дискретного времени, является степенной. Для нахождения параметров аппроксимирующей прямой использовался метод наименьших квадратов. Угловой коэффициент аппроксимирующей прямой определяет степень ц в уравнении для конкретного графика. Видно, что его значение неизменно в
10":
■Зе *
10
1-3
10
г4
10":
*******
* **
*****
* X
Н- + + "
+ + 1
С О
..............^ - - -
* * ** ..............С?)
X X х ¥
+. +
+ +
+ ■З"
1
ОООФО
оо-о-а
-ю-1
-10°
Рис. 1. Данные серии экспериментов для случая суперкритической бифуркации типа «вилка» для следующих значений: (1) - А = 10-4, Тш = 10-3; (2) - А = 10-4, Тш = 10-2; (3) - А = 10-3, Тш = 10-3; (4) - А = 10-3, Тш = 10-2. По оси абсцисс отложены значения бифуркационного параметра ц, по оси ординат - стандартное отклоне-
пределах серии. Сравнивая свободные члены в уравнениях прямых для разных графиков серии, можно получить множитель, зависящий от характеристик вносимого шума.
В итоге для транскритической бифуркации и суперкритической «вилки» получена эмпирическая зависимость вида
Л1 У VТш у/Ш'
а для бифуркации «седло-узел» -
Л1
оу уш • (9)
Подчеркнем следующее: уже можно видеть, что определение типа бифуркации в отдельных случаях возможно на достаточно большом расстоянии от точки бифуркации (вне области насыщения). Особенно удобно, что интервал корреляции наблюдаемого сигнала в этом случае сравнительно мал, и данные для статистики можно набрать достаточно быстро.
Таким образом, для нелинейных систем с непрерывным временем, претерпевающих наиболее типичные бифуркации, экспериментально проверено и подтверждено, что в окрестности точки бифуркации имеет место степенное нарастание шума. Выведенные на основании экспериментальных данных эмпирические формулы в отдельных случаях позволяют различать типы бифуркаций еще до выхода дисперсии наблюдаемого сигнала на уровень насыщения. В дальнейшем предстоит обосновать полученные зависимости теоретически, получив функции плотности распределения наблюдаемого сигнала для различных бифуркаций.
4. Теоретическое обоснование
Рассмотрим простейшую диссипативную динамическую систему вида (4), заданную уравнением
Х = х(Х — х2). (10)
Данная система удобна тем, что нет необходимости приводить ее к нормальной форме: можно видеть, что в точке X = 0 происходит суперкритическая бифуркация типа «вилка».
Потенциал (5) в данном случае имеет вид
Х
4
Хх2
и(х, Х) = т — — . (11)
Потенциальный минимум в точке х = 0, имеющий место при X < 0, после прохождения параметром X бифуркационного значения X = 0 заменяется на локальный максимум с одновременным появлением пары симметрично расположенных «потенциальных ям».
Считаем, что на систему действует равномерно распределенный броуновский шум вида (6). Пусть управляющий параметр X < 0, и X приближается к точке бифуркации с соблюдением условия (7). Попробуем оценить, как при этом меняется дисперсия наблюдаемого сигнала х(Л). Воспользуемся следующими очевидными свойствами функции у(х, X) для случая (10):
у(х, X) и Хх, |Х| >> х2 , (12)
у(х, X) и—х3, |Х| << х2 . (13)
Здесь х обозначает характерную величину х. При статистическом исследовании в качестве х можно взять, например, стандартное отклонение о х. Упрощения (12) и (13) позволяют разделить анализ поведения системы (10) по трем множествам значений параметра X: в одном из них зависимость поведения системы от параметра является наиболее простой, в другом отсутствует вообще, а в третьем упрощенный анализ невозможен.
Чтобы получить выражение для функции плотности распределения наблюдаемого сигнала, удобно вместо отдельно взятой системы вида (4) ввести в рассмотрение ансамбль таких систем. Если одной системе соответствует точка в одномерном фазовом пространстве, то образом ансамбля систем будет совокупность точечных частиц. Пусть системы подвергаются возмущению описанным выше образом, и эти возмущения не согласованы (то есть у каждой системы - свой независимый генератор шума). Движение каждой частицы будет включать две составляющие: детерминированную (движение под действием разности потенциалов) и стохастическую (броуновское движение). В случае потенциала (11) можно считать, что с течением времени установится некоторое стационарное распределение линейной плотности частиц. Предположим наличие свойства эргодичности: распределение вероятностей для положения отдельной частицы совпадает с распределением плотности частиц ансамбля.
Воспользуемся следующими определениями и обозначениями:
(М Си
р = —--линейная плотность частиц,
ах
и - скорость потока частиц (скорость частиц в точке в отсутствие шумов),
Си
д = — = ри - поток через точку,
аЛ
К(х, в, Л) - относительная скорость переноса.
Относительная скорость переноса К(х, в, Л) (плотность вероятности переноса частицы за единицу времени) вводится в качестве характеристики стохастической составляющей движения частиц. Пусть в момент времени Л координата частицы равна х. Рассмотрим вероятность случайного переноса частицы на отрезок [в, в + Св] в течение интервала времени СЛ. Определение К заключается в том, что эта вероятность равна
р = к (х, в, л)авал.
Поскольку возмущения разных частиц не согласованы, то для выбранной наугад частицы неизвестен фазовый сдвиг возмущений относительно выбранного момента начала отсчета Л = 0. Поэтому впредь будем считать стохастическую составляющую движения частиц автономной: К(х, в, Л) = К(х, в).
Также, с учетом сделанных предположений о четности функции плотности распределения возмущения, запишем
£ = |в — x|, К (х,в) = К (£).
Главной интегральной характеристикой вносимого шума является постоянная, называемая коэффициентом диффузии [15]:
2
к = I К (ШЧ •
В частности, для вносимого шума вида (6) относительная скорость переноса будет очевидным образом равна
К (Ю^ /ш (£),
Тш
а коэффициент диффузии, соответственно
к = тш2 / = 2Тш . (14)
шш
—оо
Известно [15], что для броуновских частиц поток в отсутствие потенциала зависит от плотности следующим образом:
дв(х,Л) = —к^^ • (15)
дх
Составляющая потока частиц, обусловленная уравнением (10), очевидным образом равна
др = ру.
Записав эти выражения совместно и приняв во внимание, что нас интересует установившееся стационарное распределение плотности (а значит, суммарный поток частиц д = дв + др в каждой точке равен нулю), получим
д р
д = —к^- + ур = 0. (16)
дх
Решив это дифференциальное уравнение и приняв во внимание свойство эргодичности, получим выражение для функции плотности вероятности наблюдаемого сигнала
/ (х) = /о ехр(—К и (хЛ , (17)
где /о - нормирующий множитель
1
/ 1 \
/ ехр I — и (х) I йх к
Для определения характерной величины ж воспользуемся стандартным отклонением для плотности вероятности /(ж):
х = °/
\
[ ж2/(х)(1х. (19)
Воспользуемся полученными результатами для анализа системы (10) при условии X < 0. Как и прежде, упростим систему в соответствии с (12) и (13), запишем результаты для каждого случая и определим области, где полученные результаты будут корректны.
Для системы, упрощенной в соответствии с (12), первообразная у(х) примет
вид
и (ж) = . (20)
Далее нам понадобится формула
У хв-1 ехр(—Лха)^х = Л-Р/аГ ^а)
хр_1 ехр(—Лха)йх = Л-^аГ[ а ) (21)
(можно убедиться, что заменой Лха = Ь интеграл в левой части приводится к Г-функции [16]). Согласно (18), с учетом (20) и (21), подставив в = 1, а = 2 и Л = |Х| /2к, найдем нормирующий множитель /о. Аналогичным образом, с учетом (17), подставив в (21) в = 3 (так как в формуле для дисперсии в подынтегральное выражение входит множитель х2), а = 2, получим выражение для дисперсии наблюдаемого сигнала. Окончательно, подставив выражение для нормирующего множителя в выражение для дисперсии, получим с учетом (19) и свойств Г-функции:
„/ = /к х . Р2)
Отсюда следует, что упрощение (12) применимо, пока выполняется условие
|Х| » л/к. (23)
Подставив (14), нетрудно убедиться, что (22) соответствует полученной ранее эмпирической зависимости (8).
Перейдем к рассмотрению системы, упрощенной в соответствии с (13). Выражение для потенциала при этом примет вид
и (ж) = £ . (24)
В этом случае а = 4 и Л = 1/4к. Проделав манипуляции, аналогичные описанным выше, получим
',г(3/4) „
'/ = Г Н174Гк- (25)
Условие (13) запишется в виде
Ш 2Г(3/4) ^
|Х|<< 2 Г(Т/4Г к •
На рис. 2 можно видеть пример графика функции плотности распределения наблюдаемого сигнала для данного случая. Для сравнения там же показано описанное в работе [9] экспериментально полученное распределение шума в модифицированном осцилляторе Ван дер Поля (управляющий параметр находится в точке суперкритической бифуркации типа «вилка»).
Поскольку данные на графике относятся к различным системам, различным значениям параметров вносимого шума (и, что имеет смысл подчеркнуть особо, различным видам вносимого шума), экспериментальный график пришлось подвергнуть перемасштабированию. Однако можно видеть, что теоретическая и экспериментальная зависимости совпадают (с точностью до масштаба).
Осталось рассмотреть ситуацию, когда X находится в промежуточной области и ни одно из упрощений не применимо (|Х| по порядку величины сопоставима с \/к). Полагаем, что множитель в подынтегральном выражении в (18) отвечает условию Т/к ^ Т. С учетом этого, для вычисления (18) и (19) уместно применить лемму Ватсона [16], которая позволяет для функции
№ 3000
Ш' 2000 Ш 1
Ш 1000 рв 1 1 1
/ -о/! 1 г 1
-310"4-2-Ю"4-10"4 О 10"4 2-10'4 3-10'4 х Рис. 2. Плотность распределения наблюдаемого сигнала в точке суперкритической бифуркации типа «вилка». Масштаб по координатным осям дан для теоретического графика. Теоретическая зависимость (е-е-е) получена путем подстановки (24) в (17), Of рассчитано согласно (25). Коэффициент диффузии к определен по формуле (14) для вносимого шума вида (6) с параметрами А = 10-8, Тш = 10-2
а
Ф(Л) = у хв-1^ (х)ехр(-Лха)^х
(26)
при Л ^ то воспользоваться асимптотическим разложением
ф(Л) = Т £ Л-(к±£) ^
а ^ \ а ) к!
к=0
(27)
В качестве примера найдем таким образом /о в (18). Подставим в (26) в = Т, а = 4. Положим
Л = Т/4к,
^ (х)= е-(2Л'х') х 2 • (28)
Для простоты вычисления старших производных ^ (х) введем вспомогательную переменную w и функцию п('ш):
w = у/2ЛЩх,и№) = е~'ш2, ^(к)(х) = (2Л |X|)k/2u(k)(w) •
Обратим внимание на то, что ^ - четная функция, и потому для нечетных к выполняется равенство п(к) (0) = 0. Имея в виду этот факт, получим согласно (27)-(29)
- = 2Ф(Л) = 2 £ Л-^г (^) ^дм (2л 1X1)* = I: а» |х\к. (30)
к—0 к—О
Этот степенной ряд сходится для любых X: нетрудно доказать методом индукции, что множитель п является знакопеременным и растет не быстрее, чем 22кк!.
Чтобы получить функцию плотности распределения наблюдаемого сигнала, подставим (30) в (17). Проделав простейшие преобразования, получим
ехр ( —1 и (х)
/ (х) =-V К( т ^ , (31)
^щЦ Ж
где Шл - функция, определяемая следующим образом:
ш (,) ш± £ г (2к±1) п^. (32)
^ ' \/2»—0 V 4 } (2к)! к '
Исходя из полученного закона распределения, найдем стандартное отклонение наблюдаемого сигнала О/ в (19). Для этого сначала запишем аналогичное (30) выражение для интеграла, определяющего дисперсию для плотности вероятности /(х) (разница будет заключаться лишь в в = 3). Затем выполним нормировку, разделив полученное выражение на (30) почленно. Окончательно получим
О/ = ^ X М) . (33)
где
V Г( 2к п(2к)(0) ,
Итак, поведение наблюдаемого сигнала в окрестности \ X | = у/К определяется посредством функций Ш и Шл, заданных с помощью бесконечных степенных рядов. Указанные функции являются эталонными и могут быть построены численно единственный раз, заранее. Это серьезно облегчает решение обратной задачи.
Попутно отметим, что по мере увеличения \ X \ затраты на вычисление Ш и Шл для соответствующих в быстро нарастают, и использование общих формул (31) и (33) для расчетов в области применимости упрощений (12) и (13) нецелесообразно.
На рис. 3 совместно показаны оценки дисперсии наблюдаемого сигнала согласно (22) и (25), график (33), а также измеренные экспериментально значения дисперсии. Можно видеть, что полученные зависимости (в пределах своей области применимости) находятся в согласии с экспериментальными данными. Также ясно видно, что на практике общая формула (33) не заменяет (22), так как ей можно пользоваться только в области выхода дисперсии наблюдаемого сигнала на уровень насыщения: при вычислении согласно (34) Ш(в) в конечном приближении, существует
Рис. 3. Дисперсия наблюдаемого сигнала: оценки и экспериментальные данные. Точками показаны экспериментальные значения дисперсии, полученные при численном моделировании поведения системы, заданной уравнением (10), под действием вносимого шума вида (6). Параметры шума: A = 10~3, Тш = 10~3 (к = 4.17 х 10~Б). Штрихпунктир-ная линия соответствует оценке уровня насыщения согласно (25), штриховая линия соответствует (22). Сплошной линией показан график (33) (числитель и знаменатель W(s) вычислены в 278 приближении)
такое критическое значение s, свыше которого вкладом остатков рядов в числителе и знаменателе пренебрегать уже недопустимо (на приведенном графике критическое значение выглядит как точка локального минимума). В данном случае применение (33) имеет смысл при |X| < 10\/К, так как при вычислении W(s) используются полиномы 278 степени (ограничение вызвано переполнением типа double при вычислении соответствующих факториалов).
5. Пример решения обратной задачи
Задача определения точки бифуркации в пространстве параметров и выяснения ее типа по наблюдаемой траектории динамической системы с шумом является частным случаем задачи реконструкции динамической системы по временному ряду. Так называемые обратные задачи, заключающиеся в определении некоторых количественных характеристик явления по результатам измерений их косвенных проявлений, обычно относятся к некорректным, в силу, например, разрывности обратного преобразования [17].
Приближенное решение некорректно поставленной задачи, устойчивое к малым изменениям исходных данных, может быть найдено, если условие задачи включает некоторую дополнительную (априорную) информацию об исследуемой системе. Как правило, эта информация либо позволяет сузить класс возможных решений, либо определяет общий характер решения (например, характер его гладкости). В частности, для задачи нахождения точек бифуркации и/или определения типа бифуркации, имеют значение сведения относительно параметров системы, характеристик вносимого шума и т.д. Большой класс методов приближенного решения обратных задач (так называемые методы регуляризации) был предложен в научной школе академика А.Н. Тихонова [17].
Рассмотрим вариант задачи на нахождение точки бифуркации и проверку гипотезы о ее типе на примере суперкритической бифуркации типа «вилка». (Прочие из рассмотренных простых бифуркаций связаны со сближением устойчивого и неустойчивого положений равновесия, следовательно, в их случае вблизи точки
бифуркации при достаточном времени наблюдения должен отмечаться уход сигнала х(Ь) в область аномальных значений (шумовое разрушение); анализ таких ситуаций представляет собой отдельный класс задач.)
Допустим, некоторая нелинейная система генерирует сигнал х(Ь), в котором присутствует случайная составляющая (х(Ь) суть наблюдаемый сигнал). Пусть также имеется следующая априорная информация о системе. Во-первых, система не является колебательной (это упрощает дальнейший анализ). Во-вторых, есть возможность измерения некоторой величины X, которую мы будем считать параметром системы. При этом характерное время изменения X значительно превышает интервал корреляции наблюдаемого сигнала Тх:
Тх » Тх •
Пусть также отмечено, что наблюдаемый сигнал при фиксированном X имеет стационарное распределение, причем при изменении параметра X в определенном направлении происходит рост дисперсии наблюдаемого сигнала оХ.
Требуется проверить предположение, что при дальнейшем изменении параметра, по достижении им определенного значения X,, в исследуемой системе происходит бифуркация, и определить XI,. Также требуется проверить, является ли предстоящая бифуркация суперкритической бифуркацией типа «вилка».
Предположим, что параметр находится достаточно далеко от точки бифуркации, так что система может быть линеаризована. Тогда, согласно соотношению (22), если график зависимости дисперсии наблюдаемого сигнала от параметра, оХ (X), построен в двойном логарифмическом масштабе, то точки должны ложиться на прямую с угловым коэффициентом +Т.
Если это не наблюдается, то следует найти некоторое преобразование параметра X, чтобы удовлетворить этому требованию:
X' = Щ) (35)
(фактически данное преобразование означает переградуировку шкалы прибора, измеряющего X). Например, если экспериментальные точки лежат на прямой с отличным от +Т угловым коэффициентом, то достаточно возведения X в соответствующую степень и, возможно, изменения знака. Подчеркнем, что в общем случае такое преобразование неочевидно, и его поиск представляет собой отдельную задачу.
Известно, что линейное приближение дает бесконечный рост уровня шума в точке бифуркации. Воспользуемся этим, чтобы приближенно найти точку бифуркации X,. Из (22) следует соотношение
2к
оХ -
- - Ч)'
Итак, мы получили уравнение с двумя параметрами - бифуркационным значением X, и коэффициентом диффузии к (напомним, что по условиям задачи он неизвестен). Чтобы найти значения параметров, приведем его к виду уравнения прямой:
т., X:
-2 = — X' + ^ • (36)
от
х к к
Воспользовавшись, например, методом наименьших квадратов для экспериментальных точек (Xi, (о-2)*), найдем угловой коэффициент и свободный член аппроксимирующей прямой
о-2 = aX + b, (37)
откуда, согласно (36), получим оценки для к и X'b:
1
к = — ,
a ' (38)
\'ь = bK.
Выполним следующее преобразование координат для удобства использования (33):
s = , У = 4. (39)
кк
Область применимости линейного приближения (22) в этих координатах будет определяться соотношением
S - 4)»1. (40)
Если это соотношение не выполняется, то исходное предположение об удаленности параметра от точки бифуркации неверно, и преобразование (35) и основанные на нем оценки (38) не имеют смысла; данный случай требует отдельного рассмотрения.
Следует отметить, что на практике по мере удаления параметра от точки бифуркации точность оценки коэффициента диффузии по формуле (38) возрастает (так как поведение нелинейной системы все менее отличается от поведения линеаризованной системы), в то время как точность оценки координаты точки бифуркации снижается: оценка координат центра гиперболы по удаленным точкам, близким к асимптоте (в масштабе погрешности), затруднена. Таким образом, для получения более точной оценки координаты точки бифуркации необходимо исследование наблюдаемого сигнала в области выхода его дисперсии на уровень насыщения. Относящиеся к этой области экспериментальные точки (si, yi) следует аппроксимировать функцией
y = W{s' - 4) , (41)
где сдвиг s'b - единственный параметр аппроксимации - определяет положение точки бифуркации. Если удалось выполнить аппроксимацию с заданной точностью, то можно считать, что положение точки бифуркации в новых координатах успешно определено.
Предположение о типе бифуркации можно проверять после оценки положения точки бифуркации. Для этого следует сравнить построенную по результатам наблюдения (за время, в течение которого изменение X можно считать незначительным) функцию распределения наблюдаемого сигнала с функцией (31). Поскольку в области насыщения закон распределения является уникальным для данного типа бифуркации, при совпадении сравниваемых функций предположение о типе можно считать подтвержденным.
Итак, на основании сказанного, алгоритм решения данной задачи можно сформулировать следующим образом.
1. Построить график для некоторого набора экспериментальных точек (1« А^,
1« [(°Х)*]).
2. Предположив выполнение (23), подобрать преобразование (35) таким образом, чтобы точки (1« А^, 1« [(оХ)г]) ложились на прямую с угловым коэффициентом +1.
3. Построить график для экспериментальных точек (А^, (о-2)г).
4. Методом наименьших квадратов определить коэффициенты в уравнении аппроксимирующей прямой (37).
5. По формулам (38) получить оценки для коэффициента диффузии к и координаты предполагаемой точки бифуркации "к'ь.
6. Выполнить преобразование координат согласно (39).
7. Проверить выполнение условия (40). Если условие не выполняется, то следует изменить выбор экспериментальных точек для п. 1. Например, можно попробовать выбрать точки с минимальной дисперсией наблюдаемого сигнала (предположительно, они наиболее удалены от точки бифуркации).
8. Рассчитать теоретически плотность распределения наблюдаемого сигнала для линеаризованной системы, подставив (20), А = А' — А'ь и нормирующий множитель, рассчитанный по формуле (21), в (17).
9. Проверить соответствие плотности распределения наблюдаемого сигнала, полученной экспериментально, теоретическим расчетам. При несоответствии считать алгоритм неприменимым в данном случае. Вероятной причиной несоответствия является слишком быстрое изменение параметра.
10. Если с ростом в' отмечается замедление роста дисперсии наблюдаемого шума по сравнению с оценкой для линеаризованной системы, то попытаться аппроксимировать экспериментальные точки функцией (41). В случае, если достигнута заданная точность аппроксимации, считать координату точки бифуркации определенной.
11. С учетом данных о положении точки бифуркации, рассчитать плотность распределения по формуле (31) и сравнить с экспериментальными данными. В случае совпадения считать предположение о типе бифуркации подтвержденным.
Приведем пример применения данного алгоритма. В качестве объекта для испытания была выбрана система, заданная уравнением (10). В силу выбора системы для исследования, п. 2 алгоритма выполнялся тривиальным образом. Система подвергалась воздействию периодического равномерно распределенного броуновского шума (6) со следующими характеристиками:
Тш = 10-2, А = 10-3, ош = 2.9 х 10-4, /к = 2.0 х 10-3.
Рассматривались серии по 10 точек ("•, о2), отстоящих друг от друга по А на ДА = 2 х 10-4. Дисперсия о2 наблюдаемого сигнала х(£, А^-) вычислялась по результатам измерений XI = х(и, А^), таких что ¿¿+1 — и = 10-3. Приведем результаты по двум случаям, различающимся скоростью А изменения параметра А.
В первом случае А = 2 х 10-9. Для измерений вблизи значений параметра А = —0, 02 п. 5 алгоритма дал оценки
/к = 2.0 х 10-3, АЬ = 5.3 х 10-3.
Можно видеть, что условие (40) выполнено. При этом оценка к является достаточно точной, а ошибка оценки XЬ превышает размер области насыщения. По результатам серии измерений для X на отрезке [—3 х 10_3; — 1 х 10_3] , согласно п. 10 алгоритма была получена оценка
4 = 0.05, Хь = —10"4-
Итак, ошибка составила 5% от размера области насыщения. На рис. 4, a можно видеть, что экспериментальная плотность распределения х хорошо согласуется с теоретическими расчетами по формуле (31), даже с учетом погрешности определения положения точки бифуркации. Можно считать, что гипотеза о типе предстоящей бифуркации подтверждена.
Во втором случае скорость изменения параметра в 1000 раз больше: X = 2 х 10_6. В пп. 1-2 алгоритма угловой коэффициент прямой, аппроксимирующей экспериментальные точки (^ Xj, ^ а2), был найден равным 10. Использование таких экспериментальных данных повлечет неверный выбор преобразования (35), и все дальнейшие оценки будут неадекватны. В качестве иллюстрации на рис. 4, б показана экспериментальная плотность распределения, найденная по измерениям в районе X = —0.02. Можно видеть ее несоответствие теоретическим результатам. В данном случае причиной этого является нехватка данных для расчета статистики.
Тем не менее, можно видеть, что при соблюдении условия медленного (квазистатического) изменения параметра возможно получение достаточно точного прогноза бифуркации и проверка гипотезы о ее типе весьма простым образом на основе минимума априорной информации.
Рис. 4. Плотность распределения наблюдаемого сигнала. Сравнение экспериментальной плотности распределения с теорией. Сплошной линией показан теоретический график функции плотности распределения f (х), точками - экспериментальный. а - случай медленного изменения параметра (X = 2 х 10~9); измерения в области насыщения уровня шума (X = —0.002); точка бифуркации определена с достаточной точностью, экспериментальные данные согласуются с теорией. б - быстрое изменение параметра (X = 2 х 10~6); измерения в области линейного поведения системы (X = —0.02); можно видеть, что положение точки бифуркации не может быть определено
Заключение
В настоящей работе рассмотрены вопросы поведения нелинейных динамических систем, испытывающих воздействие слабого шума, в окрестности точки бифуркации. Дан краткий обзор имеющихся результатов для систем, моделируемых нелинейными отображениями (систем с дискретным временем). Рассмотрены аналогичные явления в апериодических системах, моделируемых нелинейными дифференциальными уравнениями - системах с непрерывным временем. В частности, получена зависимость, описывающая нарастание дисперсии наблюдаемого сигнала по мере приближения к точке бифуркации.
Для изучения были выбраны наиболее типичные бифуркации - локальные бифуркации коразмерности 1. Из них были рассмотрены бифуркации типа «седло-узел», транскритическая бифуркация и две бифуркации типа «вилка» (субкритическая и суперкритическая).
Поставлен ряд экспериментов по изучению изменения наблюдаемого сигнала по мере приближения к точке бифуркации. По результатам экспериментов выведены эмпирические зависимости уровня шума от бифуркационного параметра и параметров вносимого шума для разных типов бифуркаций.
Для отдельно взятой бифуркации (суперкритическая бифуркация типа «вилка») получена точная зависимость дисперсии наблюдаемого сигнала от вносимого шума и бифуркационного параметра. Выведенная теоретически зависимость находится в согласии с полученными ранее экспериментальными данными.
Приведен пример алгоритма решения обратной задачи - обнаружения и определения типа бифуркаций посредством анализа наблюдаемого сигнала.
На основании полученных и имеющихся результатов можно сделать следующие выводы.
Во-первых, бифуркацию типа «седло-узел» можно предвидеть, основываясь на результатах наблюдений задолго до самой бифуркации, то есть, изучая рост дисперсии наблюдаемого сигнала до его выхода на уровень насыщения (рассматривая линеаризованную систему). Это особенно важно с учетом того факта, что интервал корреляции наблюдаемого сигнала при этом относительно мал, и данные для статистики можно набрать достаточно быстро.
Во-вторых, для остальных рассмотренных типов бифуркаций задача определения типа может быть решена, исходя из данных о характеристиках наблюдаемого сигнала на этапе выхода его дисперсии на уровень насыщения. Применение диффузионной модели отчасти снимает проблему увеличения интервала корреляции наблюдаемого сигнала. Действительно, если для исследования одновременно доступно множество систем, длительный сбор статистики для одной системы может быть заменен изучением установившегося распределения (или установившейся динамики, если распределение нестационарно) по ансамблю.
В заключение, считаем приятным долгом поблагодарить О.Я. Бутковского и Е.Д. Суровяткину, которые привлекли наше внимание к этому важному классу задач, за полезное обсуждение.
Работа выполнена при поддержке РФФИ (проекты № 04-01-00510, 05-0100852) и гранта Президента РФ НШ-374.2003.1.
Библиографический список
1. Малинецкий Г.Г., Потапов А.Б. Современные проблемы нелинейной динамики. М.: УРСС, 2002.
2. Кравцов Ю.А., Бильчинская С.Г., Бутковский О.Я., Рычка И.А., Суровятки-на Е.Д. Предбифуркационное усиление шума в нелинейных системах // Журнал экспериментальной и теоретической физики. 2001. Т. 120, вып. 6 (12). С. 1527-1534.
3. Малинецкий Г.Г., Подлазов А.В., Кузнецов И.В. О национальной системе научного мониторинга. Препринт Института прикладной математики им. М.В. Келдыша РАН. 2004, № 47.
4. Николис Г., Пригожин И. Самоорганизация в неравновесных системах. М.: Мир,
5. Капица С.П., Курдюмов С.П., Малинецкий Г.Г. Синергетика и прогнозы будущего. Изд. 3. М.: УРСС, 2003.
6. Wiesenfeld K. Virtual Hopf phenomenon: A new precursor of period-doubling bifurcations // Phys. Rev. A. Vol. 32, № 3, September 1985. P. 1744-1751.
7. Kravtsov Yu. A. and Surovyatkina E.D. Nonlinear saturation of prebifurcation noise amplification // Phys. Let. A. Vol. 319, Issues 3-4, 8 December 2003. P. 348-351.
8. Surovyatkina E. Prebifurcation noise amplification and noise-dependent hysteresis as indicators of bifurcations in nonlinear geophysical systems // Nonlinear Processes in Geophysics (2005) 12. P. 25-29.
9. Juel A., Darbyshire A.G., Mullin T. The effect of noise on pitchfork and Hopf bifurcations // Proc. R. Soc. Lond. A (1997) 453, 2627-2647.
10. Anishchenko VS., Neman A.B. Structure and Properties of Chaos in the Presence of Noise // «Nonlinear Dynamics of Structures» / Edited by R. Z. Sagdeev, U. Frisch, F. Hussain, S.S. Moiseev and N.S. Erokhin. Singapore-New Jersey-London-Hong Kong: World Scientific, 1991. P. 21-48.
11. Surovyatkina E.D. Rise and saturation of the correlation time near bifurcation threshold // Phys. Let. A. Vol. 329, Issue 3, 23 August 2004. P. 169-172.
12. Surovyatkina E., Kurths Ju. Pre-bifurcational noise-dependent phenomena as diagnostic instrument for revealing bifurcations in geophysical systems // Geophysical Research Abstracts. 2005. Vol. 7. P. 00462.
13. Малинецкий Г.Г. Хаос. Структуры. Вычислительный эксперимент. Введение в нелинейную динамику. М.: Наука, 1997.
14. Йосс Ж., Джозеф Д. Элементарная теория устойчивости и бифуркаций. М.: Мир, 1983.
15. Зельдович Я.Б., Мышкис А.Д. Элементы математической физики. Среда из невзаимодействующих частиц. М.: Наука, 1973.
16. Федорюк М. В. Метод перевала. М.: Наука, 1977.
17. Тихонов А.Н., Арсенин В.Я.Методы решения некорректных задач. М.: Наука,
1979.
1979.
Институт прикладной математики РАН
Поступила в редакцию 02.06.2005 После доработки 23.09.2005
BIFURCATION THEORY INVERSE PROBLEM IN A NOISY DYNAMICAL SYSTEM. EXAMPLE SOLUTION
M.-G.M. Zulpukarov, G.G. Malinetskii, A.V. Podlazov
Bifurcations in nonlinear systems with weak noise are considered. The local bifurcations are discussed: the saddle-node bifurcation, the transcritical bifurcation, the supercritical and subcritical pitchfork bifurcations.
Basing on the known prebifurcational noise rise and saturation phenomenon, the inverse problem is introduced: the problem of the bifurcation detection and determining it's type by the observed noise change (noise deviation growth fashion, saturation level, probability density). The inverse problem solution algorithm is suggested.
Зульпукаров Магомед-Герей Меджидович - родился в Москве (1978), окончил Московский инженерно-физический институт (2002). После окончания аспирантуры МИФИ в 2005 году работает младшим научным сотрудником в Институте прикладной математики им. М.В. Келдыша РАН. Е-таП:угеетё@таП.ги
Малинецкий Георгий Геннадьевич - родился в Уфе (1956), окончил физический факультет МГУ (1979), защитил кандидатскую диссертацию на тему «Нестационарные диссипативные структуры в нелинейных средах» (1982) и докторскую диссертацию на тему «Диффузионный хаос и новые типы упорядоченности в нелинейных средах» (1990) в Институте прикладной математики РАН. В настоящее время работает заместителем директора ИПМ, а также профессором кафедры прикладной математики МФТИ. Автор большого количества статей в области исследования хаоса и нелинейных явлений, а также учебника «Структуры, хаос, вычислительный эксперимент. Введение в нелинейную динамику». E-mail: gmalin@spp.keldysh.ru
Подлазов Сергей Викторович - родился в Москве (1973), окончил Московский физико-технический институт (1996), к.ф.-м.н. В настоящее время работает в Институте прикладной математики им. М.В. Келдыша РАН. Область научных интересов - нелинейная динамика, теория самоорганизованной критичности, математическая история, теоретическая демография, теория управления риском. E-mail: Tiger@Keldysh.ru
