CC BY
42
ПРИМЕР ЖЕСТКОЙ ТУРБУЛЕНТНОСТИ В СИСТЕМЕ РУСЕЛ И ДЖОКЕРОВ
М.-Г.М. Зульпукаров, Г.Г. Малинецкий, А.В. Подлазов
Рассматривается явление жесткой турбулентности - хаотического режима, отличающегося редкими катастрофическими выбросами на фоне слабых нерегулярных пространственных колебаний. Обсуждается один из примеров упрощенного качественного описания жесткой турбулентности - переключающаяся перемежаемость. Приводится пример решения обратной задачи: на основе временного ряда, порождаемого простой системой, работающей в режиме переключающейся перемежаемости (отображение Ершова), строится детерминировано-вероятностная система (система русел и джокеров), порождающая ряд с аналогичными характеристиками.
Введение
В настоящее время интенсивно ведутся исследования по применению методов нелинейной динамики в области управления риском. Возникающие при этом задачи зачастую трудно формализуемы, отличаются сложностью рассматриваемых систем и отсутствием полных и точных данных об их структуре и поведении. Эти особенности определяют изменение подхода к моделированию. От модели требуется не столько адекватность в традиционном смысле (в силу указанных причин, как правило, недостижимая), сколько воспроизведение качественных особенностей поведения наблюдаемой системы. При этом модель должна быть как можно более простой и наглядной.
Рассмотрим один из классов задач управления риском - разработку методов моделирования редких катастрофических событий. Примером таких событий являются аварии на океанских буровых платформах - сложных, масштабных и дорогостоящих сооружениях. Предполагается, что причиной подобных аварий могут быть нелинейные явления на поверхности атмосфера - океан.
Одно из центральных мест в области моделирования процессов на границе водной и воздушной сред занимает двумерное обобщение уравнения Курамото-Цузуки (Гинзбурга-Ландау) [1,2]. При определенных значениях параметров в системе наблюдается так называемая жесткая турбулентность - хаотический режим с редкими сильными выбросами (пиками), хорошо подходящий для формализованного описания некоторых катастрофических событий.
Для исследования подобных выбросов используется модификация двумерного уравнения Гинзбурга-Ландау - одномерное зависящее от времени уравнение Гинзбурга-Ландау с нелинейностью 5 порядка. Уравнение содержит малый параметр е, характеризующий интенсивность диссипации. В пределе при е ^ 0 данное уравнение переходит в нелинейное уравнение Шредингера. Известные на сегодняшний день уравнения, описывающие рост и спад выброса, получены путем исследования последнего [3].
Остается открытым вопрос описания фоновой (межпиковой) динамики и перехода от фоновой динамики к пиковой. Для этого было предложено использование более простых систем, например, систем, демонстрирующих так называемую переключающуюся перемежаемость - хаотический режим с наличием устойчивого инвариантного многообразия, со временем теряющего устойчивость. Для получения эффекта жесткой турбулентности подобную систему достаточно дополнить механизмом возвращения устойчивости инвариантному многообразию. Примером такой системы является отображение Ершова [4].
Далее обсуждается вопрос моделирования жесткой турбулентности в условиях неполноты информации о наблюдаемой системе. В качестве последней выбрана система Ершова; причем, по условиям задачи, для наблюдения доступна только одна из трех переменных.
Для воссоздания наблюдаемого поведения предлагается использовать метод русел и джокеров, основанный на конструировании модели из «разнородных» компонент - детерминированных и стохастических, непрерывных и дискретных, и т.д. [3-7]. Использование метода оправданно, когда фазовое пространство динамической системы неоднородно, и выделяются области сравнительно простого поведения (русла). В пределах русла система может быть описана простой динамической моделью. Там, где это не представляется возможным, для описания сложного поведения используется джокер - некоторый простейший алгоритм, как правило, вероятностный.
В системе Ершова рост и спад выброса представляют собой пример простого поведения, и им сопоставляются русла. В области устойчивого многообразия система демонстрирует сложное (хаотическое) поведение, для моделирования которого используется джокер.
Главное требование, предъявляемое к разрабатываемой модели, - соответствие основных характеристик временных рядов, порождаемых ею и исходной моделью. В качестве таких характеристик были выбраны распределения величины выбросов и длительности межпиковых интервалов.
1. Уравнение Курамото-Цузуки (Гинзбурга-Ландау) и жесткая турбулентность
Уравнение Курамото-Цузуки (иначе называемое уравнением Гинзбурга-Ландау) представляет собой одну из основных математических моделей, занимающих особое место в нелинейной динамике [1-3]. Данное уравнение демонстрирует большое разнообразие решений, включая стационарные, периодические и хаотические, и используется, в частности, для описания явлений, имеющих место в нелиней-
ных средах, - перехода к турбулентным режимам в гидродинамике, ионно-звуковых волн в плазме, и т.д. Нас интересует его двумерное обобщение, используемое при моделировании процессов, происходящих на границе сред: ветровых волн на воде, морфогенеза, протекающих на поверхности колебательных химических реакций [8].
Двумерное обобщение уравнения Курамото-Цузуки выглядит следующим образом:
Wt = (1 + ico) W + (1 + ici) AW - (1 + ic2) \W|2 W. (1)
Здесь W - некоторая комплекснозначная функция действительных переменных; co, ci, С2 - параметры (действительные постоянные). При определенном соотношении значений параметров в системе (1) можно наблюдать жесткую турбулентность (в англоязычных источниках используется термин hard/strong turbulence) - хаотический режим с редкими и исключительно сильными выбросами (пиками жесткой турбулентности). Для исследования данного режима используется специальная модификация уравнения Курамото-Цузуки - одномерная, но с более сильной нелинейностью
Wt = eW + (i + e) Wxx + (i - e) \W\4 W. (2)
Данное уравнение известно под названием зависящего от времени уравнения Гинзбурга-Ландау с нелинейностью 5 порядка (quintic time-dependent Ginzburg-Landau equation, QTDGL). При e = 0 оно переходит в нелинейное уравнение Шре-дингера
Wt = iWxx + i \W\4 W. (3)
Важное свойство последнего - наличие ряда законов сохранения, существенных для
описания жесткой турбулентности. Это закон сохранения:
«массы»
M = j \W\2 dx,
«импульса»
P = 1.J WW - W*Wx) dx
и «энергии»
E = J (\Wx\2 - 3 \W\6) dx.
На сегодняшний день также известны следующие свойства решений уравнения (2). Во-первых, при малых e решение уравнения (3) представляет собой хорошую аппроксимацию профиля W (x, t) в области роста больших пиков. Во-вторых, распределение максимумов пиков h описывается (в промежуточном интервале масштабов) степенной асимптотикой p (h) ~ h-a, где а - константа, а ~ 7 + 8. В области малых h имеет место иное распределение. Также максимумы пиков ограничены сверху, что, предположительно, объясняется воздействием множителя e в уравнении (2).
На практике, как правило, из всех характеристик распределения наибольший интерес представляют моменты первого и второго порядка. Здесь существенна интерпретация W. Если имеют смысл величины \W\ и \W\2, то аномально большие
пики не оказывают принципиального влияния на среднее и дисперсию, то есть распределение р (К) можно считать эффективно ограниченным. Если же рассматривается величина \Ш|4, то имеет место картина, характерная для некоторых разновидностей катастрофических событий, когда редкие аномальные выбросы оказывают существенное влияние на среднее и обращают в бесконечность дисперсию.
В численных экспериментах жесткую турбулентность можно наблюдать, решая уравнение (2) в области длины Ь с периодическими граничными условиями, с параметром е в интервале 0.01 + 0.0001. Поведение системы зависит от размера области: при Ь < 5 пики жесткой турбулентности отсутствуют, в областях промежуточной длины (10 < Ь < 80) наблюдаются возникновения одиночного пика, в более длинных областях возможно возникновение нескольких пиков одновременно. Внешне жизненный цикл одиночного пика жесткой турбулентности выглядит следующим образом (здесь и далее предполагается рассмотрение в области промежуточной длины).
• Вначале в системе наблюдаются медленные нерегулярные пространственно-временные колебания. В каждый момент времени график Ш (х) выглядит гладким, причем \Ш\ < 1.
• Затем в какой-то момент времени начинается рост пика. Скорость роста значительно превышает характерную величину Ш вне области пика. Ширина пика при этом сокращается.
• Далее, после того, как пик достигает максимальной величины (обратно зависящей от е), начинается его распад. Вначале фазы распада величина пика уменьшается, причем уменьшение происходит быстрее предшествующего роста. Ширина пика при этом увеличивается, а его форма сглаживается.
• Распад пика продолжается его превращением в быстро осциллирующий волновой пакет с колоколообразной огибающей. Высота пика при этом имеет величину порядка единиц.
• После этого, волновой пакет распространяется в ширину до заполнения всей области рассмотрения. Амплитуда волн продолжает уменьшаться, а высшие гармоники затухают.
• В конце концов, в области рассмотрения вновь наблюдаются медленные колебания. Характерное время сглаживания профиля имеет порядок е-1. Новый пик может образоваться не прежде, чем происходит сглаживание. Характерный вид профилей Ш в различных фазах пика жесткой турбулентности показан, в частности, в [9].
В ходе эволюции пика масса, импульс и энергия ведут себя следующим образом. В начале выброса происходит всплеск энергии - практически мгновенный рост на несколько порядков. Затем в ходе распада пика и далее энергия медленно и монотонно убывает. Новые пики не возникают, пока энергия не достигнет величины порядка 1.
Масса после возникновения пика начинает убывать, пока не перестанет выполняться условие Е > М. После этого начинается рост. В целом масса меняется значительно слабее, чем энергия: минимальное и максимальное значения различаются менее чем на порядок.
Импульс практически не отражает развитие пика: в расчетах, где его начальное значение близко к нулю, его дальнейшее изменение незначительно.
В настоящее время не известны простые способы точного предсказания момента и координаты начала следующего пика, а также его высоты. В этой связи представляется естественной попытка дополнить динамическое описание системы статистическим, возможно, с привлечением более простых моделей жесткой турбулентности, чем уравнение (2).
2. Переключающаяся перемежаемость и отображение Ершова
Один из вариантов упрощенного (качественного) описания жесткой турбулентности основан на использовании систем, демонстрирующих так называемую переключающуюся перемежаемость (оп-оff ¡ШетшШепсу). Приведенный термин используется для обозначения разновидности динамического режима, в котором хаотическое движение чередуется с резкими сильными выбросами. В данном случае такое поведение системы обусловлено наличием устойчивого инвариантного многообразия, теряющего устойчивость в результате бифуркации. Кроме того, система дополнена механизмом возвращения устойчивости инвариантному многообразию. Отметим следующее: предполагается, что хаотическое движение длится достаточно долго, чтобы точное определение основных характеристик выброса (момент начала, длительность, высота, направление) было невозможно.
Примером системы, в которой можно наблюдать переключающуюся перемежаемость, является отображение Ершова [3,4]. Данная система строилась с учетом следующих требований.
1. Наличие одной медленной переменной, играющей роль параметра в отображении для быстрых переменных. Данная переменная, по аналогии с уравнением (2), называется «энергией» и обозначается Е, так как предполагается, что именно энергия среди всех интегралов движения уравнения (2) играет наиболее существенную роль.
2. Наличие быстрых переменных, описывающих соответствующее межпиковой фазе движение на хаотическом аттракторе.
3. Изменение свойств аттрактора в зависимости от значения медленной переменной. При переходе медленной переменной некоторого критического значения должен происходить кризис аттрактора - появление неустойчивой области, так называемой «дырки», через которую изображающая точка покидает аттрактор, что соответствует росту пика. Здесь и далее имеется в виду изображающая точка в фазовом пространстве быстрых переменных.
4. Прекращение роста пика и возвращение изображающей точки на аттрактор (аттрактор вновь становится глобально притягивающим). Существует два очевидных способа определения момента возвращения аттрактору свойства притяжения. Первый, параметрический, предполагает зависимость от значения энергии (рост прекращается при превышении энергией некоторого порогового значения). Второй, динамический, предполагает зависимость от значений быстрых переменных (рост прекращается при попадании изображающей точки в заданную область фазового пространства). Динамический способ считается в большей степени соответствующим сути поведения системы, заданной уравнением (2). Здесь следует отметить, что для простоты было бы естественно ограничиться единственной быстрой переменной. Однако в сочетании с динамическим способом ограничения роста пика это
привело бы к колебаниям вокруг порогового значения (при превышении порогового значения рост прекращается, начинается спад, после обратного перехода через порог вновь начинается рост и так далее). Поэтому необходимы наличие как минимум двух быстрых переменных и несколько более сложный механизм остановки роста пика. Например, следующий: изображающая точка, приближается к аттрактору по первой координате и продолжает удаляться по второй, пока не произойдет возврат по первой координате в область аттрактора. После этого начинается возврат по второй координате.
5. Максимально возможная простота функции отображения быстрых переменных. Целесообразно использование кусочно-линейной функции.
Система строится следующим образом. Есть две быстрые переменные, х и у, и соответствующие им два связанных одномерных отображения с параметром Е. Когда Е выше порогового значения Екр, имеется хаотический аттрактор с некоторой конечной областью притяжения. При уменьшении Е до порогового значения и далее происходит кризис аттрактора - границы аттрактора и его области притяжения соприкасаются, после чего появляется «дырка». При попадании в «дырку» изображающей точки начинается уход последней на бесконечность. Затем при превышении переменной у некоторого порогового значения укр начинается возвращение на аттрактор по координате х. Когда изображающая точка достигает аттрактора (и более не покидает его) по координате х, начинается возвращение по координате у.
Далее энергия Е рассматривается не как параметр, но как медленная переменная. Пока изображающая точка движется в пределах аттрактора, Е медленно убывает, сходясь к некоторой неподвижной точке Е* < Екр. При этом, когда условие Е > Екр нарушается, происходит кризис аттрактора, изображающая точка начинает уходить на бесконечность, и вместе с этим начинает расти Е. Как только условие Е > Екр вновь начинает выполняться, «дырка» закрывается, и аттрактор становится локально притягивающим (росту значений быстрых переменных это не препятствует). После остановки роста пика по координате х рост энергии также прекращается, и она вновь начинает медленно убывать.
Соотношение быстрого и медленного масштаба времени определяется параметром е, задающим скорость изменения энергии. Стоит отметить, что при е ^ 0 рост пика начинается при значении энергии, очень мало отличающемся от критического (Е ^ Екр — 0). Поскольку «дырка» возникает при Е = Екр и увеличивается в размерах по мере уменьшения Е, постольку в данном случае размеры «дырки» оказываются малы. То есть разброс начальных значений при росте пика мал и для энергии, и для быстрых переменных. Это не соответствует реально наблюдаемому поведению динамической системы, заданной уравнением (2).
Данная проблема решается путем некоторого усложнения механизма образования пика. При переходе энергии в область Е < Екр открывается «дырка» и начинается рост только по одной координате х. «Дырка» по у открывается только при достижении х некоторого порогового значения хкр. Параметры отображений подбираются таким образом, чтобы переменная х достигала порогового значения за количество итераций, достаточное для образования разброса по у порядка размера аттрактора. Изображающая точка попадает в область открывшейся «дырки» по у спустя некоторое время Т\. Разброс по Т тем больше, чем больше разброс по у к моменту открытия «дырки», и чем меньше размер «дырки».
Время, которое проходит до момента достижения порогового значения по у, обозначается Т2. Таким образом, после достижения жкр рост по х продолжается в течение времени Т1 + Т2. Если размер «дырки» по у мал, то разброс по Т2 также мал, и разброс максимальных значений по х определяется разбросом по Т1.
Окончательно, с учетом изложенных соображений, записывается трехмерное отображение следующего вида:
хп+1 — / (хпт кх,пт ап) + Ул ,| И
1 + \Уп\
Уп+1 — / (Уп, ку,п, ао) - 70 л гХп—т,
1+ е\Хп \
Еп+1 — Еп — е (Еп + хп),
Е
а (Е) — аж-
кх —
ку —
кх > 1, У > укр;
кх < 1, У < укр;
кУ+ > 1, х < хкр
ку < 1, х > хкр
(4)
Е + Екр( ^ — 1 а
кр
Здесь величины у и уо - параметры связи одномерных отображений быстрых переменных, а / - кусочно-линейная функция, определяющая характер этих отображений,
2 1 + а
/ (х, к, а) — <
а — 1 3
х +
х +
1 — а' 2а
2а + 2 2а + 2' кх + к — ^,
х > а, — 1 < х < а, х < -1.
(5)
График этой функции показан на рис. 1, а. Можно видеть, что одномерное отображение, заданное этой функцией, всегда имеет неустойчивую неподвижную точку. При к > 1 появляется еще одна неустойчивая неподвижная точка, обозначаемая х*. Эта точка представляет собой границу области притяжения аттрактора данного отображения.
Рис. 1. Графики функций, входящих в отображение Ершова: а - функция /, входящая в одномерные отображения быстрых переменных, х* - неподвижная точка; б - зависимость параметра а фунции / от энергии
Функция а (Е) определяет влияние энергии на параметр а функции / и, как следствие, на геометрию аттрактора и его области притяжения, ограничивая это влияние в области больших значений Е. Величины акр, а- параметры этой функции; ее график показан на рис. 1, б.
Динамика системы (4) по построению подразделяется на следующие фазы. В межпиковой фазе происходит хаотическое движение на аттракторе; значения обеих быстрых переменных имеют порядок 1. Энергия медленно убывает, соответственно значение параметра а в отображении для х увеличивается и, следовательно, наклон правой ветви / возрастает.
Фаза выброса по координате х начинается, когда энергия, а значит, и наклон правой ветви / достигает порогового значения. В результате этого значение хп может оказаться в области х < х*, что означает выход изображающей точки из области притяжения аттрактора и экспоненциальный рост: хп ~ (к+) . Энергия быстро растет вместе с ростом |х|, Еп ~ е (к+) . В это же время динамика по координате у приближенно описывается отображением
Уп+1 = / (уп, к+, а0) ,
пока выполняется условие 1хп1 ^ 1/е. Аттрактор этого отображения - локально притягивающий, уход с него на бесконечность невозможен, несмотря на то, что левая ветвь / стала растягивающей, когда 1хп1 превысило |хкр|.
Фаза выброса по у наступает, когда координата х изображающей точки достигает величин порядка 1/е. По мере роста х, появляется необходимость учета добавочного члена в отображении, описывающем динамику по координате у,
уп+1 = / (уп,к+,ао) + гп, гп = Уо 1 , •
1 + 1ехп1
Параметры подобраны таким образом, что отображение
уп+1 = / (уп,к+ ,ао) + 70>
соответствующее хп = 1/е, уже не имеет аттрактора. «Дырка» в аттракторе возникает при достижении добавочным членом Гп некоторого критического значения
Г = Г =_1/2_
Г = " (k + —1)1+° — 1 ^
1 — ао
Ее размеры невелики, так что изображающая точка попадает в область «дырки» не сразу, но спустя некоторое количество итераций. Все это время добавочный член и, как следствие, геометрия «дырки» продолжают меняться, так что момент начала экспоненциального роста уп ~ и его начальные условия имеют некоторый
разброс.
Фаза спада по х наступает, как только |у| превышает критическое значение |укр|. Левая ветвь / в отображении по х становится сжимающей, вследствие чего 1х„\ начинает убывать: хп ~ (к-)п. Параметры подбираются таким образом, что убывает быстрее, чем рос до того. Рост энергии сменяется медленным уменьшением: Еп ~ (1 — е)п. Рост по у при этом продолжается.
Фаза распада пика начинается при переходе уменьшающегося | через пороговое значение |хкр|. Левая ветвь / в отображении по у становится сжимающей, рост
Igbi
lg E 14.0
12.0
10.0
8.0
6.0
4.0
2.0
0.0
-2.0
а
14
!: \ — ig|x| .....Ш ....... lg E
i • 1
Ii.W.
f; v /i: Л 1 •• '1
f' / I i * 'i f:! i 1 / .i ii i '' ' l
1 -// V t 1
"ii X. ' 1:' '"Щ fli Л .<J\ Äi» ÄJrflT*. 4 ,-UU
8450
8500
8550
800 1000 t. mc
Рис. 2. Жесткая турбулентность в отображении Ершова и в эксперименте: а - графики E, x, у в полулогарифмическом масштабе для абсолютных значений. Цифрами обозначены фазы: 1 - межпиковая, 2 - выброс по x, 3 - выброс по у, 4 - спад по x, 5 - распад пика. б - экспериментальный график зависимости напряженности электрического поля от времени для ленгмюровского солитона огибающей, связанного с локальным радиовсплеском III типа [10]
\yn\ сменяется еще более быстрым уменьшением: yn ~ {к—n. Пока значение \yn\ превышает пороговое \укр\, уменьшение \x\ продолжается. Параметры должны быть подобраны так, чтобы изображающая точка успела войти по координате x в область аттрактора к моменту перехода yn через порог, когда правая ветвь f в отображении по x станет растягивающей. Уменьшаясь по модулю, yn входит в область аттрактора. Система вновь вступает в межпиковую фазу.
Общий вид происходящего в системе показан на рис. 2, а. На рис. 2, б для сравнения приведен пик жесткой турбулентности, зарегистрированный в описанном в [10] эксперименте по измерению напряженности электрического поля при локальном солнечном радиовсплеске III типа. Показанное соответствует следующим приведенным в [3] значениям параметров:
е = 0.02, ao = 0.5, акр = -0.2668, ах = -0.8, Ер = 0.3, у = 0.122,
4
Yo = 0.5, к+ = 4, к- =0.1, к+ = 2, к- =0.1, Хкр = -з, Укр = -100.
(6)
Данный набор значений параметров используется и во всех дальнейших расчетах. Для системы (4) аналитически (в предположении е ^ 0) показано, что распре-
деление максимумов1 пиков по координатам x и у - степенное
ax (-x)-bx , x < Bx, 0, x > Bx,
Px
(x) =
(7)
1
1 Здесь и далее слова «максимум пика» используются для обозначения максимума по абсолютной величине.
а распределение длительностей межпиковых интервалов - экспоненциальное
Г апе~ь"п, п > Вп, Рм (п) = (8)
[0, п < Вп
(здесь ах, Ьх, Вх, ап, Ьп, Вп - некоторые константы). В дальнейшем, будем рассматривать эти распределения как основные характеристики режима жесткой турбулентности.
3. Реконструкция отображения Ершова: случай одной переменной
Итак, пусть имеется система (4), работающая в режиме жесткой турбулентности. Пусть система рассматривается как «черный ящик», и из трех порождаемых ею временных рядов, Еп, хп и уп, для наблюдения доступен только один - хп. Наша задача заключается в построении простой (модельной) системы, порождающей временной ряд с аналогичными характеристиками.
Заметим, что теоретически, в условиях поставленной задачи, можно воспользоваться теоремой Такенса для восстановления фазового портрета по одномерной наблюдаемой [3]. Однако из-за наличия двух различных масштабов времени (наблюдаются резкие и короткие пики и длинные межпиковые интервалы, в пределах которых энергия Е меняется медленно по сравнению с переменными х и у) при таком подходе потребовалась бы выборка неприемлемого объема. Таким образом, на практике необходимо использование какого-либо упрощения.
Поскольку информация наблюдателя об исходной системе неполна, будем учитывать как динамические, так и статистические характеристики временного ряда. Следовательно, при построении модельной системы можно совместно использовать (хотя бы, для простоты) динамические и вероятностные компоненты, то есть задействовать метод русел и джокеров. Иными словами, предполагается, что данный метод может упростить построение модельной системы, и в то же время дать возможность использовать максимум доступной информации о наблюдаемой системе. В этом состоит отличие от распространенной практики моделирования, предписывающей игнорировать как можно большую часть доступной информации, с целью упрощения модели.
Первичный анализ информации о наблюдаемой системе. Обсудим, какие характеристики временного ряда хп могли бы быть рассмотрены при построении модели. Главная особенность временного ряда - наличие редких сильных выбросов, чередующихся с длительными межпиковыми интервалами. Следовательно, среди интересующих нас характеристик временного ряда можно выделить характеристики выбросов (распределение максимумов пиков, их продолжительность, динамика роста и спада) и характеристики межпиковых интервалов (продолжительность, межпиковая динамика). Также нас может заинтересовать процесс перехода от межпиковой (фоновой) динамики к росту пика и обратно - от спада к фоновой динамике.
Прежде всего нам понадобится правило, позволяющее отличить значения хп, относящиеся к фазе выброса, от значений хп, относящихся к межпиковым интервалам. На рис. 2 можно видеть, что первые превышают вторые на много порядков
величины. В то же время, длительность выброса крайне мала по сравнению с длительностью межпикового интервала. Отсюда следует, что медиана хм временного ряда [\хп\] даст оценку характерной величины переменной состояния в межпиковой фазе. Для выделения пиковых значений хп может использоваться пороговое значение - некоторая величина хп > 0, превышающая медиану на один или несколько порядков. Надпороговые значения хп, то есть, такие, что \хп\ > хп, считаются относящимися к фазе выброса.
Согласно (7), эмпирическая функция плотности распределения максимумов пиков в логарифмических координатах (1п(—х), 1п р) должна (в идеальном случае) представлять собой прямую. Оценив свободный член со и угловой коэффициент С1 соответствующей линейной функции, например, методом наименьших квадратов, получим оценки параметров ах и Ьх
ах = еС0,
Ьх = —С1.
(9)
Здесь и далее надстрочный знак «» используется для обозначения оценки. Чтобы оценить Вх, воспользуемся условием нормировки
Г Вх
ах ( „Л-Ьх + 1
Рхтах (х) ййх = ( —х)
ах (—Вх)-Ьх + 1 = 1
—оо
Ьх — 1
(предполагается выполнение условия сходимости интеграла Ьх > 1). Решив уравнение относительно Вх и подставив оценки ах и Ьх, получим
Вх = — Г^ . (10)
Ьх — 1
Точно так же, согласно (8), построив эмпирическую функцию плотности распределения длительностей межпиковых интервалов в полулогарифмических координатах (п, 1п р) и воспользовавшись формулами, полностью аналогичными (9), получим оценки ап и Ьп. Решив уравнение нормировки
ап -Ьпп
рм (п) йп = е
J —Ьп
= а^. е-ЬпВп = 1
Вп Ьп
относительно Вп и подставив ап и Ьп, получим оценку
1 ап
Вп = г 1п (11)
Ьп Ьп
Рост пика, согласно (4) и (5), в области больших \х\ описывается отображением
хп+1 кх хп + кх 0. (12)
1
Здесь мы пренебрегли членом, содержащим yn. В силу его сравнительной малости при заданном наборе значений параметров (6), а также ненаблюдаемости переменной y по условиям задачи, воздействие данного члена можно считать слабой флуктуацией. Отображение (12) можно переписать в виде
(Хп - x*) = (xo - x*) (к+)n , (13)
где x* - его неподвижная точка. При отсутствии флуктуаций эта точка ограничивает область роста справа. Таким образом, точки на графике в полулогарифмических координатах (n, ln (—xn)), относящиеся к росту пика, в области больших \x\ (порядка максимума пика) должны аппроксимироваться некоторой прямой. Имеем оценку к+ для коэффициента к+
к+ = eCl, (14)
где Ci - угловой коэффициент данной прямой.
Спад пика в области больших \x\ описывается уравнением
xn ~ (кх ) , xn+1 = кх xn (15)
(начальное значение для данного отображения - точка максимума пика).
Оценку x* неподвижной точки x* несложно вычислить, зная оценки свободного члена Со и углового коэффициента ci линейной функции (12)
x* = Со • (16)
1 - С1
Оценки Со и ci можно получить, аппроксимировав прямой точки, относящиеся к росту пика, на графике (xn,xn+i). Данные оценки будут весьма неточными по ряду причин. Во-первых, в силу экспоненциального роста xn, для построения аппроксимирующей прямой можно использовать очень небольшое количество соседних точек: 3 или 4 точки для набора значений параметров (6). Некоторого повышения точности можно добиться, накапливая точки, относящиеся к разным пикам, а также заменив в формуле (16) ci более точной оценкой к+. Во-вторых, оценка по точкам \x\ ~ \x*\ неточна, так как в этой области возрастает влияние слабых флуктуаций в (12). В-третьих, оценка по точкам \x\ ^ \x*\ неточна, поскольку погрешность оценки свободного члена линейной функции нарастает по мере удаления от начала координат.
С другой стороны, нас интересует не столько получение точной оценки x*, сколько воспроизведение динамики роста пика. Поэтому, оставим за собой право определять x*, исходя из соображений простоты построения модели.
Процедура получения оценки к х для коэффициента к х полностью аналогична описанной выше. Роль неподвижной точки сжимающего отображения не представляется существенной.
Фоновая динамика внешне выглядит как случайное отображение некоторой ограниченной области на себя, при выходе из этой области начинается рост пика. Будем считать эту область отрезком [Dx>min, Dx)max]. Поскольку пик растет в отрицательном направлении, величину Dx max можно оценить как максимальное значение xn за время наблюдения
Dx,max = max {xn} . (17)
{n}
Величину можно оценить как минимальное значение хп в промежутке между
какими-либо двумя соседними пиками
О X
= ШШ {хп N <п< N2 }
(18)
где N1 - момент окончания спада первого пика, а N2 - момент начала роста второго.
Для определения N1 и N2 можно поступить следующим образом. Рассматривая поведение системы в области больших \х\ (|х| ^ хп), нетрудно выделить некоторый момент N0, относящийся к спаду первого пика. Считая, что в процессе спада значение хп непрерывно возрастает, N1 можно оценить как первый, начиная с N0, момент нарушения этой закономерности
хп < хп+1,П = ^,N0 + 1,... , хМ1 > хИ1 + 1.
Рис. 3. Фоновая динамика, оценки границ области фоновой динамики. График хп(п) в межпиковой фазе. Цифрами обозначены точки: 1 - точка с ординатой х = Д
согласно (17), 2 - точка с ж,тт согласно (18), 3 - точка с
/х,тах
ординатой х = Д абсциссой п = N0, 4 - точка с абсциссой п = М1 согласно (19), 5 - точка с абсциссой п = N3, 6 -точка с абсциссой п = N2 согласно (20). Пунктирной линией обозначен порог \х\ = хп
(19)
Аналогично, выделяется некоторый момент N3, относящийся к росту второго пика, и определяется N2
хп-1 > хп,П = ^,N3 - 1,... ,
хм2-1 ^ хм2 ■
) ) (20)
Фоновая динамика и процедура получения оценок Е>х,тп и Е>х,тах показана на рис. 3. Для повышения точности оценки Охт\п можно повторить описанную процедуру для нескольких межпиковых интервалов.
Предварительные соображения по схеме русел и джокеров. Итак, анализируя временной ряд хп, можно выделить:
• фазу фоновой динамики,
• событие перехода от фоновой динамики к росту пика,
• фазу роста,
• событие остановки роста (достижения максимума пика),
• фазу спада,
• событие перехода от спада к фоновой динамике.
Динамика в фазах роста и спада описывается простыми уравнениями (13) и (15). Исходя из этого, имеет смысл назначить для моделирования роста и спада два русла, С1 и С2, соответственно.
Фоновая динамика хаотична, а момент перехода от фоновой динамики к росту пика - случайная величина (переменная Е, инициирующая переход, не наблюдается). Следовательно, имеет смысл использовать джокер для их описания. Обозначим данный джокер Как будет видно далее, этот же джокер можно использовать и для описания события перехода от спада к фоновой динамике.
Достижение максимума пика - также случайное событие (смена динамики не следует из уравнения роста; инициирующая ее переменная у не наблюдается). Для его описания используем джокер, обозначенный
Попутно заметим, что алгоритмы моделирования случайных событий должны быть рассчитаны на использование стандартного датчика случайных чисел, показания которого считаются случайной величиной, равномерно распределенной на отрезке [0,1]. В дальнейших выкладках мы будем обозначать ее как Z, а соответствующую функцию плотности распределения как
( 1, г е [0,1], Р^ (*) = < (21)
\ 0, г е [0,1].
Русла С1 и С2. Динамика на русле С1 задается уравнением роста пика (13). Из этого уравнения следует, что рост пика не может начинаться с точки х = х*. Следовательно, имеет смысл ограничить область русла неравенством
х ^^ -Ах << .
В некоторой окрестности х* влияние слабых флуктуаций является определяющим, и динамика хаотична. Область роста пика примыкает к области фоновой динамики слева, и с точки зрения наблюдателя, данная окрестность х* относится к области фоновой динамики. Фактически, точка х = Ах ограничивает окрестность слева, поэтому можно считать, что
Ах = ^х,шт- (22)
Будем полагать, что начальное значение хо, согласно намеченной схеме русел и джокеров, представляет собой результат работы джокера ^ - случайную величину. Функцию плотности распределения х0 обозначим рХ0 (х).
Динамика на русле С2 задается уравнением спада (15). Поскольку рост и спад пика происходят в одной и той же области фазового пространства, область русла С2 будет совпадать с областью русла С1. Таким образом, требуется уточнять, какое из русел в данный момент считается действующим.
Джокер Разработка алгоритма джокера подразумевает решение двух задач. Во-первых, необходимо задать правило выбора момента срабатывания джокера. Во-вторых, требуется определить саму процедуру срабатывания.
Вторая задача в данном конкретном случае представляется более простой. Результат срабатывания джокера ^ должен заключаться в том, что экспоненциальный рост значения наблюдаемой переменной сменяется экспоненциальным спадом. Таким образом, процедура срабатывания представляет собой изменение уравнения отображения (в терминах намеченной схемы русел и джокеров - переключение между руслами С1 и С2). Здесь просматривается аналогия с предложенным в [7] джокером типа «шов», также меняющим уравнение движения без перемещения изображающей точки.
Разница между джокером типа «шов» и джокером ^ должна состоять в определении момента срабатывания. Первый срабатывает при пересечении фазовой траекторией линии, представлающей собой область джокера. В случае второго, это про-
стое правило неприменимо. Координата изображающей точки в момент срабатывания джокера представляет собой максимум пика - случайную величину. Продолжительность роста пика - также случайная величина, что не дает возможности «настроить» джокер на срабатывание в фиксированные моменты времени.
Таким образом, момент срабатывания джокера 32 должен определятся некоторым вероятностным предикатом. Его аргументами могут быть как время п, так и координата х. Остановимся на определении момента срабатывания, исходя из координаты, так как этот способ в большей степени соответствует «традиционному» понятию области джокера.
В результате, приходим к следующему алгоритму. В процессе роста пика на каждой итерации принимается решение: либо срабатывает джокер и рост пика останавливается (пик достигает максимального значения), либо рост продолжается в соответствии с уравнением русла Сь Вероятность срабатывания джокера зависит от значения х на данной итерации, р = рз2 (х). Задача сводится к поиску функции рз2 (х), обеспечивающей распределение максимумов, соответствующее распределению (7). По этой причине будем ссылаться на (7) как на целевое распределение.
Для удобства дальнейших выкладок, перейдем к логарифмической системе координат
X = 1о§д,+ (х* — х), х < х*. (23)
Искомую зависимость вероятности срабатывания джокера от координаты далее будем обозначать рз2 (х). Удобство (23) состоит в том, что отображение, описывающее рост пика, примет очень простой вид (рис. 4)
хп+1 = хп + 1. (24)
Рис. 4. Рост пика, срабатывание джокера и спад пика в логарифмической системе координат. Штриховая линия обозначает отображение, описывающее рост пика (движение в русле С1). Штрих-пунктирная линия обозначает отображение, описывающее спад пика (движение в русле С2). Символом ▲ отмечены значения переменной состояния в процессе роста. Символом ■ отмечена точка срабатывания джокера. Символом ▼ отмечены значения переменной состояния в процессе спада
Найдем вид целевого распределения в новых координатах. Для этого понадобится формула распределения функции случайной величины. Пусть имеется некоторая случайная величина V с функцией плотности распределения ру (у). Пусть задана некоторая детерминированная функция д, непрерывная, строго монотонная и имеющая обратную функцию д-1 на всей области рассмотрения. Функция от случайной величины Ш = д (V) будет случайной величиной с функцией плотности распределения
р^ (ад) = ру (д-1
(25)
Сопоставим хтах V, хтах Ш, а функцию (23) д. Отметим, что в области х < Вх < х*, во-первых, функция д отвечает всем предъявляемым требованиям, а во-вторых, величиной х* в выражении (23) можно пренебречь в силу неравенства
|х*| < |Вх|. Отсюда
д-1 (х) = — (к+)х ,
( х
— д-1 (х) = — (к+)х\п к+,
Рхтах {д-1 (х)) = ах (к+)
-ЬхХ
Далее, можно видеть, что произошел переход от степенного распределения к экспоненциальному
Pimax (X) = ax = ax ln к+,
axe
0,
x > Bx, x < Bx,
(26)
bx = (bx - 1) ln k+.
В дальнейшем понадобится соотношение коэффициентов, которое нетрудно получить из уравнения нормировки,
ax - bxe
bx Bx
(27)
Отметим следующее важное свойство отображения (24). Пусть дан некото-
рый интервал Qx,Qx + sj, Qx — const, s < 1. Очевидно, что он не пересекается
со своим образом Qx + 1,Q x + s + 1 ) и прообра3ом Qx 1,Qx + s - 1). Обозначим pqqx s (X) условную функцию плотности распределения значений X (условием является попадание в данный интервал в процессе роста пика). Из (25) следует, что условная функция плотности распределения на ¿-образах (прообразах) Qx + ¿,Qx + s + ^ данного интервала имеет вид
pQx+i,s (X)= PQXs (X - ¿) G Z
Интервал Qx + u,Qx + u + s j, u G R, в общем случае, не является образом
(прообразом) Q x, Qx + s ) . Однако при условии s - 1 первый интервал можно
разбить на два подинтервала - образы (прообразы) соответствующих подинтервалов второго интервала с числом итераций, различающимся на единицу. Следовательно, обозначив pqx (X) = pqx 1 (X), можно записать соотношение
P^x (X) = PQx (X - Sx + Qx)
(28)
В дальнейших выкладках будем предполагать, что в процессе роста пика изображающая точка гарантированно попадает в интервал Вх ,ВХ + и, как следствие, функция плотности распределения значений х на данном интервале совпадает с РВ (X). Для этого необходимо, чтобы начальные значения находились полностью в области х < Вх, то есть
Вх
J Pxo (X) dX = 1
(хотя бы, с некоторой приемлемой точностью). Тогда функция плотности распреде-
ления максимумов пиков в интервале X G
Bx,Bx + 1) будет равна
Pxmax (X) = PBx (X) № (X) > X G
Bx,Bx + 1
(29)
x
Согласно (28), для образов указанного интервала можно записать последовательность уравнений следующего вида:
Pimax (X + 1) = РBx (X)(1 - J (X)) J (X + ^ >
Pimax (X + 2)= Р Bx (X) (1 - J (X)) (1 - J (X + ^ pJ2 (X + 2)
(30)
Pimax (X + П) = PBx (X)
n— 1
П (1 - PJ2 (X + i))
i=0
Pj2 (X + n) ,
Рис. 5. Типичный вид распределения максимумов пиков для джокера Л2. На графике показан общий вид распределения максимумов пиков, обеспечиваемый джокером Л2 в соответствии с формулами (29) и (30). Для наглядности, в качестве рвх (X) и рз2 (х) взяты произвольные функции
где х е Вх,Вх + 1) ,п е N. Здесь учтено то обстоятельство, что изображающая точка попадает в следующий интервал, только если джокер не срабатывает на данном интервале. Общий вид функции рхтах (х) показан на рис. 5. Отсюда имеем уравнение
Pimax (Х + П + 1)= Pimax (Х + n) (1 - J (Х + n))
Pj2 (X + n + 1) Pj2 (X + n) '
X e Bi,Bi + 1
n e N.
Учитывая, что из (26) следует
Pimax (X + n +1)= e —bx , X > Bi, n e N,
PiBmax (X + n)
получаем для p j2 (X) рекуррентное соотношение
pj2 (X + n + 1) =
J (X + n) e—~bx
(1 - J (X + n))
X
или, что эквивалентно,
J (X +1) =
J (X)
(1 - J (X))
Bi,Bi + 1
X > BT.
nN
(31)
Будем искать решение p j2 (X) этого уравнения в простейшем виде pj2 (X) = pi = const, X > Bi. Подставив в (31), получим выражение для искомой константы
Pi = 1 - e—Bx.
Проверим, что найденное нами решение не нарушает условие нормировки
Bx+1
J PBx (X) dX = 1.
(32)
(33)
x
Из (26) и (29) следует выражение для функции плотности распределения
Подставив его в (33), получим уравнение
X Е
Вх, Вх + 1 .
вх+1
РВ (X) йх = ах —е
х Рх — Ьх
ах
вх+1
Вх Р х Ьх
1е-
Вх
е-ЬхВх = 1
Преобразовав с учетом (27), получим необходимое условие для (33)
Рх = 1 - е-Ьх.
Можно видеть, что (32) не противоречит данному условию.
Теперь, рассмотрим процесс роста пика в целом. В силу определенной свободы выбора границы Ах русла С1, а также распределения начальных значений Хо, положим, что все возможные начальные значения находятся в интервале Ах, Ах + .
Отсюда, с учетом (28), (27) и (34), следует выражение для распределения начальных значений
Рго(Х) = Рлх (Х) = ахе~
х е
Ах, Ах + 1
ЪхеЬхЛх
(35)
ах —
Записав для интервала
Рх
Ах, Ах + 1
последовательность уравнений, аналогичную
(29)-(30), с ограничением X + п < Вх и сравнив с (26), можно доопределить р^2 (X) на интервале Ах,Вх ) нулем. Окончательно, формально распространив область джо-
кера 32 на все русло С1 и выполнив обратное преобразование координат, имеем
Рл2 (х) =
Рх = Рх = 1 - е Ь 0, х е (Вх, Ах].
х < Вх,
(36)
Джокер Прежде всего на джокер ^ возлагается ответственность за моделирование фоновой динамики и за переход от фоновой динамики к росту пика. Сначала рассмотрим данные задачи по отдельности.
Пусть моделирование фоновой динамики заключается в случайном отображении отрезка Дх>тт, Дх,тах] на себя. Для простоты положим, что хп на каждой итерации присваивается значение некоторой случайной величины, функция плотности распределения которой полностью находится в данном отрезке (то есть, значения хп независимы в совокупности). С той же целью используем случайную величину Z с распределением (21). Итак, в межпиковой фазе
хп = (Д
х,тах Дх,тт) х Z + Дх
п
0, 1,
(37)
Следующая задача джокера 31 - обеспечение длительности межпиковой фазы в соответствии с экспоненциальным распределением (8). Пусть на каждой
х
х
п-й итерации возможно прекращение межпиковой фазы джокером с вероятностью, определяемой некоторой функцией рз1 (п). Чтобы найти эту функцию, рассмотрим пуассоновское приближение биномиального распределения
Рт (п) =
(пр)г
т!
-е-пр, п » 1, р < 1,
где Рт (п) - вероятность ровно т успехов в п испытаниях Бернулли, а р - вероятность успеха в отдельном испытании. Известно, что промежуток (количество испытаний) между успешными испытаниями распределен экспоненциально
Рм (п) = ре
—пр
р < 1.
Следовательно, взяв вероятность прекращения джокером межпиковой фазы на данной итерации р = Ьп, получим искомый вид распределения межпиковых интервалов2. Доопределив функцию рз1 (п) нулем в начале межпиковой фазы, окончательно получим
( рп = Ьп, п > Вп,
р31 (п) = (38)
[ 0, п = 0,...,Вп - 1.
Теперь, рассмотрим переход от фоновой динамики к росту пика. В этот момент, джокер 31 вместо выполнения очередной итерации, согласно (37), должен перевести изображающую точку в исток русла Сь Таким образом, имеет смысл возложить на ^ обязанность формирования в истоке русла распределения начальных значений вида (35), необходимого для получения заданного распределения максимумов пиков.
Для этого требуется решить частную задачу поиска зависимости Хо = д z (^), позволяющей получить распределение (35) на основе (21). Согласно (25), имеем уравнение
Рхо(х) = Рz (д—1 (х))
Тхд—1 (х)
Тхд—1 (х)
0,
д—1 (х) е [о, 1] д—1 (х) е [о, 1]
Предположив, что выражение под знаком модуля положительно на всем интервале Ах, Ах + , перейдем к следующему уравнению:
^ д—1 (х) = ахе-ЬхХ,
д— 1 (X) е [о, 1],
х
А х
,А х + 1) .
Проинтегрировав, получим
д— 1 (X) = - -Ле—ЬхХ + С. Ьх
Замечание. На первый взгляд, в аналогичном случае для джокера Л2 целевое распределение зависит не только от рл2 (х), но и от распределения начальных значений х0. Поэтому, данная схема, предположительно, не даст существенного уменьшения объема выкладок.
Можно видеть, что на интервале
Ах,Ах + ^ полученная функция строго
возрастает, что соответствует сделанному выше предположению относительно положительности выражения под знаком модуля (рис. 6).
Чтобы найти константу интегрирования С, запишем граничные условия:
д—1 Ш = 0,
д—1 (^Ах + 1)=1.
Отметим, что одно из граничных условий избыточно и, следовательно, может использоваться для проверки непротиворечивости решения. Имеем
х „~ьхлх
С = Ь е Ьх
Выполнив обращение найденной функции, получим
Х0 = дЕ (Я) = -11п -2 + е"^ Ьх \ «х
о- ъхлх
Рис. 6. Получение распределения начальных значений на основе равномерного распределения. 1 -функция плотности распределения (21) случайной величины Z; 2 - искомая зависимость (39) начального значения х0 от случайной величины Z; 3 - заданная функция плотности распределения (35) начального значения х0
или после обратного преобразования координат и подстановки определения ах из (35) получим
Х0 = Х* - (&+)Лх (1 - рхЯ) Ьх
(39)
Теперь можно определить значение Х* таким образом, чтобы воспроизвести динамику в начале фазы роста пика. Как было отмечено выше, получение точной оценки х* затруднено, и задача заключается, скорее, в воспроизведении внешнего вида фронта пика. Поэтому имеет смысл при определении оценки пользоваться наиболее грубыми и простыми методами.
Рассмотрим начало роста некоторого пика. Момент выхода хп из области
[Дс,т1п, Дт,тах] обозначим ^старт
ХМстарт — 1 — Dx,min,
(40)
Момент превышения абсолютной величиной хп порогового значения хп обозначим N
\Хп\ < ХШ Жстарт < П < Ып,
(41)
\Хп\ > Хп, П = Ып. Число итераций между этими двумя моментами обозначим Ыфронт
Ыфронт — Ып Ыстарт-
(42)
+ х
1п к
Среднее по всем наблюдаемым пикам значение Лфронт обозначим (Лфронт). Воспользуемся обозначением (Лфронт,м) для модельной системы и (Лфронт,н) для наблюдаемой системы. Будем считать, что внешний вид фронта пика воспроизводится приемлемо, если
( Афронт,м) ~ ( Афронт,н) •
Значение (^фронт) модельной системы можно грубо оценить, воспользовавшись аналогичным (13) равенством
(хп - х,) = (Ы - х,) (к+)(Мфронт'м) , (43)
где (хо) обозначает математическое ожидание начального значения. Чтобы найти последнее, выполним обратное (23) преобразования координат для распределения начальных значений (35). Согласно (25), имеем
рзд (х) = /( ,ч1_Ь—Ц-^-г~ (х* - х)-Ьх ' х е (к+Ах>Ах] •
{к+)1-Ьх - 1) (х, - Ах)1-Ьх
Отсюда,
Ах
(хо) = J Рхо (х) йх = х, - Со (х, - Ах) •
(44)
где
1 - Ьх (к+)2-Ьх - 1
С0 = 2-ЬХ (к+)1-Ьх - 1 • (45)
Приравняв (Лфронт,н) и (Лфронт,м), подставив (22) и (44) в (43) и решив уравнение относительно х,, получим оценку
Сп (к+)(^фрон1'н) П —т
- _ С0 \кх )_Пх,Ш1П хп
х* = /м \ • ( )
со (к+) (^фронт'н) - 1
Наконец еще одна потенциальная обязанность джокера ^ - прекращение спада пика согласно уравнению русла С2. Фаза спада считается законченной, когда изображающая точка достигает области фоновой динамики. По этой причине, а также с учетом написанного выше имеет смысл объявить отрезок [Пх Ш1П, Пх шах] областью джокера Таким образом, двигаясь в соответствии с уравнением русла С2, изображающая точка покинет русло и перейдет в область джокера, после чего произойдет активизация последнего.
Построение системы русел и джокеров. Окончательно цикл работы системы русел и джокеров выглядит следующим образом. В межпиковом интервале в области джокера ^ выполняются итерации случайного отображения (37) в течение времени, определяемого вероятностью начала выброса (38). Выброс начинается с задания в истоке русла С1 джокером ^ начального значения согласно (39). После этого выполняются итерации соответствующего руслу С1 отображения (13). При этом вероятность прекращения роста определяется джокером ^ согласно (36).
кх Ах
Рис. 7. Система русел и джокеров. Прямой штриховкой обозначены совпадающие области русел С1 и С2 и джокера Л2. Косой штриховкой обозначена область джокера Л1. Символом ▲ отмечены значения переменной состояния в процессе роста пика. Символом ■ отмечена точка срабатывания джокера Л2. Символом ▼ отмечены значения переменной состояния в процессе спада пика. Цифрами обозначены: 1 - срабатывание джокера Л1: случайное перемещение изображающей точки в соответствии с уравнением (37); 2 - срабатывание джокера Л1: перевод изображающей точки в исток русла С1; 3 - движение в русле С1 в соответствии с его уравнением (13); 4 - срабатывание джокера Л2: перевод изображающей точки в русло С2; 5 - движение в русле С2 в соответствии с его уравнением (15); 6 - изображающая точка вновь попадает в область джокера Л1, движение в соответствии с уравнением русла С2 прекращается
В момент прекращения роста русло С1 перестает считаться действующим, активизируется русло С2, и начинают выполняться итерации соответствующего ему отображения (15). Это продолжается до тех пор, пока изображающая точка не достигнет области джокера после чего начинается новый межпиковый интервал (рис. 7).
В соответствии с постановкой задачи, исходные данные для построения системы представляют собой временной ряд хп (предполагается, что данных достаточно для получения всех необходимых статистических оценок).
Ниже приводится предлагаемый порядок построения системы русел и джокеров.
1. Вычисляется медиана хм временного ряда {|хп|} и пороговое значение хп.
2. Определяются границы областей русел и джокеров.
2.1. По формуле (17) вычисляется оценка Ох,тах (используется весь имеющийся временной ряд, либо данные по одному или нескольким межпиковым интервалам).
2.2. По надпороговым значениям хп выделяются два следующих непосредственно друг за другом пика.
2.3. Выделяются соответствующие надпороговым значениям моменты N0 и N3 в фазах спада первого пика и роста второго, соответственно.
2.4. По формуле (19), с использованием значений хп, начиная с п = определяется момент N1 начала межпикового интервала.
2.5. По формуле (20), с использованием значений хп, начиная с п = N3 в сторону уменьшения п, определяется момент N2 конца межпикового интервала.
2.6. По формуле (18) вычисляется оценка Ох>т\п. Это же значение, согласно (22), является правой границей области русел Ах.
2.7. Возможен повтор данной процедуры для остальных межпиковых интервалов с целью уточнения оценок.
3. Вычисляются оценки ап, Ьп и Вп.
3.1. Строится гистограмма длительностей межпиковых интервалов в полулогарифмических координатах (п, 1пр).
3.2. Для построенной гистограммы прокладывается аппроксимирующая прямая.
3.3. Вычисляются оценки ап и Ьп путем подстановки коэффициентов аппроксимирующей прямой в формулу, аналогичную (9).
3.4. Вычисляется оценка Вп путем подстановки полученных оценок ап иЬп в (11).
4. Вычисляются оценки ах, Ъх и Вх.
4.1. Строится гистограмма максимумов пиков хтах в логарифмических координатах (1п (-х), 1п р).
4.2. Для построенной гистограммы прокладывается аппроксимирующая прямая.
4.3. Вычисляются оценки ах и Ъх путем подстановки коэффициентов аппроксимирующей прямой в (9).
4.4. Вычисляется оценка Вх путем подстановки полученных оценок ах и Ъх в (10).
5. Вычисляются оценки к+ и к-.
5.1. Для одного из пиков строится график в полулогарифмических координатах (п, 1п(-хп)).
5.2. Прокладывается аппроксимирующая прямая по точкам в области больших 1х1 (порядка Вх ), относящимся к фазе роста пика.
5.3. Вычисляется оценка к+ путем подстановки углового коэффициента аппроксимирующей прямой в (14).
5.4. Прокладывается аппроксимирующая прямая по точкам в области больших |х| , относящимся к фазе спада пика.
5.5. Вычисляется оценка к- путем подстановки углового коэффициента аппроксимирующей прямой в формулу (14).
5.6. Возможен повтор данной процедуры для других пиков, с целью уточнения оценок.
6. Вычисляется оценка х*.
6.1. Выполняется подстановка оценки Вхт\п в правило (40). Согласно данному правилу, на основе значений хп, начиная с момента п = N2 конца одного из межпиковых интервалов, определяется значение ^харт.
6.2. Определяется значение Nп по правилу (41) и вычисляется значение ^ронт по формуле (42).
6.3. Аналогичным образом определяются значения ^ронт для остальных пиков. Вычисляется среднее значение (^р^^).
6.4. Вычисляется значение Со путем подстановки оценок Ъх и к+ в формулу (45).
6.5. Вычисляется оценка х* путем подстановки значений с0, (^р^^), хп и оценок Дц)тт, к+ в формулу (46).
7. Вычисляются коэффициенты Ах, Вх, ах, Ъх и вероятность рх.
7.1. Ах и Вх - путем подстановки х = Ах и х = Вх в (23).
7.2. ах и Ъх - путем подстановки ах, кх и к+ в (26).
7.3. рх - путем подстановки оценки Ъх в (32).
8. Окончательно определяются уравнения схемы русел и джокеров.
8.1. Уравнение (13) русла Ci - путем подстановки X* и к,+.
8.2. Уравнение (15) русла C2 - путем подстановки к-.
8.3. Функция отображения (37) джокера Ji - путем подстановки Dx,min и
Dx, max.
8.4. Функция вероятности прекращения межпиковой фазы (38) джокера Ji -путем подстановки bn и Bn.
8.5. Функция формирования джокером Ji начального значения для роста пика (39) - путем подстановки k+, Ax, px и bx.
8.6. Функция вероятности срабатывания джокера J2 (36) - путем подстановки Ax и Bx.
Результаты моделирования, сравнение. Для проверки предложенной модели использовался временной ряд xn, порождаемый системой (4) с набором параметров (6) на протяжении 5107 итераций. За время наблюдения произошел 8221 выброс.
На основе собранной статистики были получены следующие оценки: Dx,max = 1-06, Dx,min = -1.17, ^ = 4.33 • 10-4, bn = 2.5 • 10-4, Bn = 2199, ax = 9.36 • 10-2, bx = 1-04, Bx = -2.47 • 109, k+ = 4.0, к- = 0.1, X* = -1.16. Вероятностные параметры джокеров Ji и J2 составили pn = 2.5 • 10-4 и px = 5.35 • 10-2, соответственно.
Общий вид динамики модельной системы Ji-Ci-J2-C2 показан на рис. 8 (ср. рис. 2).
Гистограммы абсолютных величин максимумов пиков и длительностей межпиковых интервалов в системах (4) и Ji-Ci-J2-C2 показаны на рис. 9. Можно видеть, что модельная система соответствует наблюдаемой (исходной) системе по выбранным для сравнения статистическим характеристикам.
Рис. 8. Динамика системы русел и джокеров. а - фоновая динамика и одиночный пик в системе Ji-C1-J2-C2, полулогарифмический масштаб. б - фронт пика, подробно: сплошная линия соответствует системе J1-C1-J2-C2, пунктирная линия - наблюдаемой системе; штрихпунктирной линией показан уровень х = Dmin, штриховой - уровень x = хп
20.0 60.0 100.0 140.0 180.0 1п|х| 0.0 1.0 2.0 3.0 и10г4
а б
Рис. 9. Сравнение статистических характеристик исходной и модельной систем: а - гистограммы максимумов пиков, логарифмический масштаб; б - гистограммы длительностей межпиковых интервалов, полулогарифмический масштаб. Сплошная линия соответствует системе Л1-С1-Л2-С2, пунктирная линия - наблюдаемой системе
Заключение
Рассмотрен пример построения модели системы с жесткой турбулентностью в условиях неполноты доступной информации. Исходные данные представляют собой один из трех временных рядов, порождаемых трехмерным отображением Ершова (4). Характерная особенность рассмотренной системы заключается в наличии двух различных масштабов времени, что делает использование теоремы Такенса невозможным на практике. С учетом этого при решении задачи использовался метод русел и джокеров, заключающийся в совмещении динамических и вероятностных методов моделирования.
Наблюдаемый временной ряд отличается наличием длительных промежутков хаотической динамики, чередующихся с сильными выбросами (экспоненциальный рост, сменяющийся экспоненциальным спадом). В качестве главного показателя соответствия системы русел и джокеров исходной системе выбрано сходство распределений максимумов пиков и длительностей межпиковых интервалов порождаемых временных рядов.
Предложена схема, включающая два русла С1 и С2 и два джокера 31 и 32. Уравнения русел описывают экспоненциальный рост и спад. Области русел совпадают, поэтому специально оговаривается, какое из русел в данный момент считается действующим.
Джокер 3 отвечает за моделирование фоновой динамики и за переключение между различными типами динамики. Фоновая динамика моделируется путем случайного отображения области джокера на себя (используется равномерно распределенная случайная величина). Выброс начинается в результате перевода джокером изображающей точки в исток русла С1 (после этого русло С1 объявляется действующим). Спад заканчивается после достижения изображающей точкой, следующей в русле С2, области джокера 31. Область джокера смежна с областью русел.
Джокер 32 представляет собой вероятностное правило определения момента перехода от роста к спаду, то есть в результате его срабатывания русло С1 перестает, а русло С2 начинает считаться действующим. Для формирования требуемого вида распределения максимумов пиков оказалось достаточно задаться постоянной веро-
ятностью срабатывания джокера на каждой итерации. Область джокера совпадает с областью русел.
Следует подчеркнуть две особенности предложенной схемы. Во-первых, она предусматривает совпадение областей русел и джокеров, и в этом смысле представляет собой дальнейшее развитие схемы, предложенной в [7]. Во-вторых, в отличие от предлагавшихся ранее схем, наибольшую часть времени изображающая точка проводит не в области русел, а в области джокера.
Среди возможных направлений доработки предложенной системы можно отметить изменение джокера Ji для более точного моделирования фоновой динамики. Например, воспроизведение распределения значений наблюдаемой величины, ее функции автокорреляции и т.д. Также возможно решение более сложной задачи, когда для наблюдения доступны две переменные x и y из трех. В этом случае потребуется не только обеспечение заданного вида распределения для каждой переменной, но и соблюдение соотношения между параметрами распределений, учет распределения запаздывания выбросов по переменной y относительно выбросов по переменной x, учет кросс-корреляции переменных в межпиковой фазе и т.д.
Библиографический список
1. Kuramoto Y. Chemical oscillations, waves, and turbulence. Berlin: Springer-Verlag, 1984. 156 p.
2. Ахромеева Т.С., Курдюмов С.П., Малинецкий Г.Г., Самарский А.А. Нестационарные структуры и диффузионный хаос. М.: Наука, 1992. 544 с.
3. Малинецкий Г.Г., Потапов А.Б. Современные проблемы нелинейной динамики. М.: УРСС, 2002. 360 с.
4. Малинецкий Г.Г., Потапов А.Б., Подлазов А.В. Нелинейная динамика: подходы, результаты, надежды. М.: УРСС, 2006. 280 с.
5. Белайчук Л.В., Малинецкий Г.Г. Проделки джокеров на одномерных отображениях. Препринт / Институт прикладной математики им. М.В. Келдыша РАН. М., 1997. № 24. 29 с.
6. Малинецкий Г.Г., Потапов А.Б. Русла и джокеры: о новых методах прогноза поведения сложных систем. Препринт / Институт прикладной математики им. М.В. Келдыша РАН. М., 1998. № 32.
7. Зульпукаров М.-Г.М., Малинецкий Г.Г., Подлазов А.В. Метод русел и джокеров на примере исследования системы Розенцвейга-Макартура. Препринт / Институт прикладной математики им. М.В. Келдыша РАН. М., 2006. № 21. 32 с.
8. Гарел Д., Гарел О. Колебательные химические реакции. М.: Мир, 1986. 148 с.
9. Малинецкий Г.Г. Управление риском. Риск, устойчивое развитие, синергетика. М.: Наука, 2000.
10. Thejappa G., MacDowall R.J.Ulysses observations of nonlinear wave-wave interactions in the source regions of type III solar radio bursts // Journal of Astrophysics and Astronomy. 2000, September, December. Vol. 21. № 3, 4. P. 447.
Институт прикладной математики Поступила в редакцию 10.12.2007
им. М.В. Келдыша РАН
AN EXAMPLE OF HARD TURBULENCE IN THE CHANNELS AND JOKERS SYSTEM
M.-G.M. Zulpukarov, G.G. Malinetskii, A.V. Podlazov
Hard turbulence, a chaotic mode distinguished by infrequent catastrophic outbreaks on the weak irregular space oscillation background, is considered. On-off intermittency as one of the possible ways of simplified qualitative definition of hard turbulence is discussed. The paper introduces an example solution of the inverse problem, constructing a deterministic-probabilistic system (a channels and jokers system) generating time series with characteristics similar to the ones of the time series generated by a simple system working in on-off intermittency mode (the Ershov mapping).
Зульпукаров Магомед-Герей Меджидович - родился в 1978 году в Москве. Окончил Московский инженерно-физический институт (2002) и аспирантуру МИФИ (2005). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (2007). Работает младшим научным сотрудником в Институте прикладной математики им. М.В. Келдыша РАН. Автор и соавтор 9 научных работ.
Малинецкий Георгий Геннадьевич - родился в Уфе (1956), окончил физический факультет МГУ (1979), защитил кандидатскую диссертацию на тему «Нестационарные диссипативные структуры в нелинейных средах» (1982) и докторскую диссертацию на тему «Диффузионный хаос и новые типы упорядоченности в нелинейных средах» (1990) в Институте прикладной математики РАН. В настоящее время работает заместителем директора ИПМ, а также профессором кафедры прикладной математики МФТИ. Автор большого количества статей в области исследования хаоса и нелинейных явлений, а также учебника «Структуры, хаос, вычислительный эксперимент. Введение в нелинейную динамику». E-mail: gmalin@spp.keldysh.ru
Подлазов Сергей Викторович - родился в Москве (1973), окончил Московский физико-технический институт (1996), к.ф.-м.н. В настоящее время работает в Институте прикладной математики им. М.В. Келдыша РАН. Область научных интересов - нелинейная динамика, теория самоорганизованной критичности, математическая история, теоретическая демография, теория управления риском. E-mail: Tiger@Keldysh.ru
