CC BY
20
Ампилова Наталья Борисовна, Терентьев Сергей Валерьевич
УДК 517.9
РЕАЛИЗАЦИЯ ИНТЕРВАЛЬНОГО АЛГОРИТМА ЛОКАЛИЗАЦИИ ИНВАРИАНТНЫХ МНОЖЕСТВ ДИНАМИЧЕСКИХ СИСТЕМ
Аннотация
Работа посвящена разработке и реализации основанного на интервальной арифметике алгоритма локализации инвариантных множеств динамических систем. Используется метод аппроксимации системы с помощью символического образа, представляющего собой ориентированный граф, построенный по системе. Ячейки разбиения рассматриваются как интервальные вектора в пространстве соответствующей размерности. Оценка параметров символического образа позволяет определить точность построения. Приведены результаты численных экспериментов и сравнение с алгоритмами локализации, основанными на обычной арифметике.
Ключевые слова: динамические системы, инвариантные множества, символический образ, компьютерное моделирование, интервальная арифметика.
ВВЕДЕНИЕ
Теория динамических систем является фундаментальной математической дисциплиной, тесно связанной с большинством областей математики. Динамические системы все чаще используются во многих областях как модели реальных процессов. Они делятся на два основных класса: непрерывные и дискретные. Непрерывные динамические системы задаются автономными системами дифференциальных уравнений, а дискретные - рекуррентными соотношениями (разностными уравнениями).
Гладкие динамические системы с непрерывным временем являются основными моделями в классической механике со времен Ньютона. В XX веке они вновь привлекли к себе внимание благодаря работе А. Пуанкаре, который при решении проблемы трех тел обнаружил, что в ней возникают сложные режимы, появление которых невозможно объяснить, опираясь на принципы классической механики, а аналитическое решение получить не удается. Такие режимы характеризуются
© Н.Б. Ампилова, C.B. Терентьев, 2011
существованием апериодических траекторий, при этом решения с близкими начальными данными оказываются существенно различными (чувствительная зависимость решений от начальных данных). А. Пуанкаре разработал геометрический подход, который положил начало современной динамике, изучающей долгосрочное асимптотическое поведение системы при помощи методов, не основанных на знании ее решений в явном виде [5].
В 1963 году Е. Лоренц рассмотрел гидродинамическую модель для предсказания погоды. Рассмотренная им система дифференциальных уравнений 3 порядка ведет себя непредсказуемо (хаотично) в том смысле, что малейшие ошибки в измерении состояния атмосферы быстро растут, что делает долгосрочный прогноз погоды неправильным. Системы с хаотическим поведением траекторий были обнаружены во многих прикладных задачах, в динамике популяций и в хорошо известных уравнениях Дуффинга.
Развитие компьютерной техники привело к активному использованию компьютерного моделирования при изучении динамики систем со сложным поведени-
ем траекторий, в частности, ее инвариантных множеств. Один из основных классов компьютерных методов исследования динамических систем представляют так называемые методы, основанные на множествах (set-oriented methods), которые используют конечное разбиение фазового пространства на клетки (ячейки) и отслеживают динамику системы по попаданию точек траекторий исследуемой системы в эти клетки. Для выбранной области фазового пространства K строятся последовательные подразбиения, из которых удаляются те ячейки, образы которых заведомо не принадлежат K. При стремлении диаметра ячеек к нулю мы можем получать все более точный фазовый портрет. На этих методах основано большинство алгоритмов локализации разных видов инвариантных множеств и, в частности, аттракторов динамических систем ([1, 10, 12]).
В 1983 г. Г. Осипенко [2] предложил метод символического образа, который являлся естественным развитием описанной выше идеи и позволил применить к исследованию динамических систем методы теории графов. Символический образ есть конечная аппроксимация динамической системы, представляющая собой ориентированный граф. Он строится по заданному покрытию фазового пространства ячейками Ci, вершины графа соответствуют ячейкам, дуги - связям между ними, а именно: вершины i и j соединяются дугой, если образ ячейки Ci при действии динамической системы пересекается с ячейкой Cj. В работах [2, 3] приводятся доказательства сходимости метода при решении различных задач, например, при построении инвариантных, в частности, цепно-рекуррентных множеств динамических систем.
Наряду с развитием методов исследования динамических систем, основанных на обычной арифметике, все чаще появляются работы, использующие арифметику интервальную, основанную на проведении вычислений с получением точных границ для искомого значения, то есть получения в качестве ответа интервала, содержащего это значение.
Компьютерные методы решения различных задач с помощью интервальной арифметики приведены в работах [7, 6, 8, 17]. Следует отметить работу У. Такера [18], в которой было доказано, что для классических значений параметров в системе Лоренца имеется аттрактор. У. Та-кер разработал алгоритм, основанный на использовании интервальной арифметики с направленным округлением, позволяющий находить точные решения большого класса обыкновенных дифференциальных уравнений. Подход Такера при исследовании системы Лоренца был основан на комбинации аналитических и компьютерных методов: строгих вычислений и теории нормальных форм. Авторы работы [16] успешно применяют методы интервальной арифметики для компьютерного доказательства существования хаоса в системе Лоренца. Статья [13] посвящена подробному исследованию отображения Икеда, оценке верхней границы инвариантного множества, локализации неблуждающего множества и периодических орбит небольшого периода. Последнее выполняется с использованием интервального оператора Кравчука. Следует отметить, что в работах, посвященнх данной тематике, методы интервальной арифметики используются в основном как доказательные вычисления.
В [9] была решена задача создания программного комплекса для компьютерного исследования динамических систем методами интервальной арифметики и символического образа. Ячейка покрытия естественным образом может быть представлена интервальным вектором, а интервальное расширение функций, описывающих систему, можно рассматривать как «отображение символического образа». При этом сам граф символического образа не строится, связь ячейка-образ(-то есть дуга на графе) хранится в специальной структуре данных. При последовательном подразбиении покрытия вычисляются параметры символического образа, позволяющие оценить точность приближения.В данной работе мы описываем различные реализации алгоритмов
локализации инвариантных множеств динамических систем. Основное внимание уделено оптимизациям в представлении данных и процессе вычислений. Приведены результаты численных экспериментов, убедительно показывающие преимущества данного алгоритма по сравнению с известными алгоритмами, основанными на обычной арифметике.
1. ОСНОВНЫЕ ОПРЕДЕЛЕНИЯ
Символический образ динамической системы строится по заданному покрытию фазового пространства и системе и представляет собой ориентированный граф, вершины которого соответствуют ячейкам покрытия.
Для динамической системы /, покрытия {М}, г =1, ..., п на символическом образе С строится дуга (г, у) если/(М.) п/(М.) ф 0.
Для непрерывной системы символический образ строится для отображения сдвига.
Символический образ характеризуется следующими параметрами:
- максимальный диаметр ячейки -
- нижняя грань г - наименьшее из расстояний между образами ячеек и теми ячейками, с которыми образы не пересекаются,
- максимальный из диаметров образов ячеек - д.
Определение 1. Интервалом называется множество вида:
х = [х, х] = е Я | х < ~ < х}.
Радиус интервала х определяется следующим образом:
rad (х) = х х
2 ■
Пространство всех интервалов над Я обозначается 1Я.
Определение 2. Пусть х1, ..., хт е 1Я.
Тогда х = (х1,...,хт) называется интервальным вектором. Пространство интер-вальныгх векторов размерности т обозначается 1Ят.
Радиус интервального вектора определяется как
тай(х) = (тай(х1), ..., тай(хт)),
хр хт е 1Я-
Нетрудно заметить, что интервальный вектор имеет довольно естественное геометрическое представление: прямое произведение составляющих его интервалов.
Определение 3. Пусть а, Ь с 1Я. Тогда функция
р(а,Ь) = тах{| а - Ь |,| а - Ь |} задает расстояние между интервалами а и Ь.
Так введенное расстояние является метрикой на 1Я. Для а, Ь е 1Ят определим р(а, Ь) = тахге ^ р (а,, Ь), где (а1, Ь.) - компоненты соответствующих интервальных векторов.
Пусть /: Ят ® Ят, а е 1Ят. Определим область значений / на а :
V/, а)= {/(х), х е а }.
Определение 4. [17] Интервальная (то есть выгчисленная по правилам интервальной арифметики) функция В называется интервальным расширением/, если для любого интервального вектора х вытолняет-ся включение V(/, х) с В( х).
Поскольку для вещественнозначной функции / существует много интервальных расширений ([15]), удобно выбрать то, которое получается, если вычислять / по правилам интервальной арифметики (естественное расширение). Обозначим такое расширение /(а) . По известному свойству монотонности по включению V(/, а) с/(а). Разобьем а на более мелкие части а1, ..., ап. Тогда из указанного свойства следует
V/, а) с ип=0 /(а1) с/(а). Поскольку при стремлении диаметра подразбиения к нулю р(У (/, а), / (а)) оценивается через диаметр а (как множества), а для дифференцируемой функции / через квадрат диаметра [18], то при последовательном подразбиении интервальных векторов мы получаем оценку для V (/, а) с любой заданной точностью.
Таким образом, интерпретируя ячейки покрытия как интервальные вектора в пространстве соответствующей размерности, можно рассматривать интервальное расширение исходной функции f как «отображение символического образа» и строить образ ячейки, используя правила и функции интервальной арифметики.
2. ОПИСАНИЕ АЛГОРИТМА
Пусть {M} есть покрытие фазового пространства M динамической системы f, где i = 1, ..., n. Обозначим за Mt интервальный векто(р, соответствующий ячейке Mi, тогда Mt = aj х...хaj, где aj - j-ая компонента интервального вектора. Таким образом, ячейка покрытия является многомерным параллелепипедом. Введем также обозначения
с, = {M\М} n f (M .) Ф0], С,. = {jMj n f(M,)Ф0] и R = Uj£CMj.
Определим диаметр ячейки diam (Mi) = maxke[i,m]{2 * rad(M,)}k и обозначим через d наибольший из таких диаметров. Пусть r = minikdist (f (M,), Mk), где f (M,) n Mk = 0 и dist обозначает ха-усдорфово расстояние между множествами. (Отметим, что в данном случае удобнее использовать именно это расстояние, а не интервальное, так как последнее дает большее значение для непересекающихся интервалов.) Таким образом, для оценки параметров символического образа нужно вычислять d и r.
Пусть Dj = {Mj,...,Mn} - начальное покрытие компакта M ячейками M, с IRm; i = 1, ..., n. Для каждой ячейки Mi вычислим множества Ri . Построим множество-представление графа G1, соответствующее символическому образу G1 отображения f относител ьно покрытия D1, а именно GDj = Uie[i, „]R. Это множество содержит все ячейки покрытия, которые участвуют в построении графа G1. Вычисляем значения d1 и r1. На следую-
щем шаге разбиение применяется к элементам множества О01, в результате повторения описанной процедуры получается множество-представление Ов для графа 02. Вычисляем d2, г2. Множества локализуются с точностью е > d2.
На к-ом шаге мы рассматриваем разбиение вида
^ = {¥л,...,м1т\мп сс0к_1>,
строим множество-представление О
Dk
И
вычисляем соответствующие dk и rk. Согласно [3], при стремлении диаметра разбиения к нулю мы локализуем инвариантные множества с точностью e > dk.
3. ПРЕДСТАВЛЕНИЕ СИСТЕМЫ И ДИНАМИЧЕСКАЯ ГЕНЕРАЦИЯ КОДА
Вычисление арифметических выражений, описывающих систему, может быть оптимизировано следующим образом: программный модуль, создаваемый динамически на основе специального описания системы с помощью специально разработанной библиотеки. Такой подход позволяет избежать многократного обхода дерева разбора при вычислении арифметического выражения [9]. При реализации алгоритма были созданы 2 библиотеки динамической генерации кода:
- для Microsoft .Net платформы [9];
- с поддержкой POSIX стандартов.
При проектировании данных учитывались следующие требования:
1) быстродействие,
2) инкапсуляция данных,
3) поддержка вещественной и интервальной арифметик,
4) простота и гибкость в использовании.
Вычисление арифметических выражений, а также доступ к описанию системы должны осуществляться максимально эффективно. Описание системы должно инкапсулироваться в одном объекте на протяжении всего жизненного цикла программы. Такой подход позволяет достичь не только предельной локализации изменений (при их необходимости), но и прогно-зируемости последствий изменений. Опи-
сания системы для двух разных арифметик выглядят одинаково, арифметические выражения вычисляются однообразно. Поэтому для достижения этой тождественности в реализации методы создания класса-системы должны опираться на использование принципов обобщенного программирования, то есть объект-системы можно применять к различным арифметикам, не меняя внутреннюю архитектуру этой структуры данных. При использовании модуля динамической генерации кода необходимо ознакомиться только с форматом входных данных, способом инициализации и с набором исключительных ситуаций, возникающих при некорректной работе с библиотекой, не вникая в подробности внутреннего устройства библиотеки.
Платформа Microsoft .Net. В качестве промежуточного языка используется C# [22]. Благодаря технологии .Net Code Document Object Model [20], компилятор и компоновщик представлены в виде некоторого объекта в адресном пространстве процесса. С помощью технологии .Net reflection [19] сборка (Assembly [22]) производится сразу в адресном пространстве запущенного инструмента без создания физического DLL-файла (библиотеки). Взаимодействие с объектом-системой осуществляется через специальный параметризованный интерфейс:
internai interface ISystem<T> {
IInterval [] Bounds { get; } int Dimension { get; } IFunction<T> Function { get; } string InputCoveragePath { get; } string OutputCoveragePath { get; }
}
Многомерная функция представлена параметризованным интерфейсом, который имеет следующий вид:
public interface IFunction<T> {
T[] In { set; } T[] Out { get; } Evaluate(); Node Root { get; }
}
Аргументы функции задаются через свойство In, а возвращаемое значение хранится в Out.
Свойство Root является корневым узлом дерева разбора заданного арифметического выражения. Это дерево используется при генерации кода класса-функции, а также при решении ряда специальных задач. Например, при вычислении константы Липшица (ее формула рекурсивна, поэтому нужна рекурсивная структура данных).
Поддержка POSIX стандартов. Эта
библиотека реализована исключительно на C++ , с поддержкой POSIX-платформ. В качестве компилятора была использована утилита g++.exe, а компоновщика -dllwrap.exe. Кроме того, для обработки сценариев компиляции применяется утилита make.exe. Указанные программы являются частью набора инструментов MinGW [26].
В этой реализации система также является параметризованным классом. Параметр - операнд в арифметике( тип double или Interval). Класс-система имеет следующий интерфейс доступа:
template<typename T> class System
{
public:
SystemBounds: :const_iterator boundsBegin() const;
SystemBounds::const_iterator boundsEnd()const;
T* Evaluate( T* );
std::string GetInputCoveragePath()const;
std::string GetOutputCoveragePath()const;
};
4. ПРЕДСТАВЛЕНИЕ ЯЧЕЙКИ И РАЗЛИЧНЫЕ РЕАЛИЗАЦИИ АЛГОРИТМА
Основной задачей при локализации инвариантных множеств является реализация алгоритма построения образа ячейки. Существует несколько таких алгоритмов, которые представляют ячейку некоторым произвольным набором точек (например вершинами и центром), а в качестве образа рассматривают объединение ячеек, в которые эти образы попадают. При этом алгоритмы, которые являются наиболее предпочтительными с теорети-
ческой точки зрения, оказываются наиболее трудоемкими, поскольку значительно увеличивается количество вычислений функций, описывающих систему.
Таким образом, использование интервальной арифметики является хорошей альтернативой арифметики обычной. Мы вычисляем образ ячейки как образ множества, то есть вычисляем значение интервальной функции. Для заданной системы время выполнения операции вычисления образа ячейки есть некоторая константа.
Оптимизация при выборе структуры данных позволяет значительно сократить время вычислений по сравнению с алгоритмами локализации, основанными на обычной арифметике. Анализ алгоритмов локализации показывает, что наибольшее время занимает поиск ячейки в покрытии. Мы рассматриваем ячейки одного размера. В описываемых далее реализациях использовались представления ячейки целым числом, многомерным целочисленным вектором и элементом R-дерева. Проводились сравнения быстродействия. При решении задачи о локализации инвариантных множеств было реализовано 4 версии алгоритма:
1) базовая (I0) [11],
2) модификация с распределенными вычислениями и использованием модуля динамической генерации кода (I) [9],
3) модификация I с поддержкой R-де-ревьев (II) [9],
4) версия III с поддержкой оптимального хранения ячеек.
Все приложения реализованы на языке C++. При реализации алгоритма была использована библиотека BOOST интервальной арифметики [24]. Текущая реализация этой библиотеки работает на платформах POSIX, Win32 и Macintosh Carbon.
Визуализация инвариантных множеств осуществляется с помощью технологии GNUPLOT [25].
4.1. ПРИЛОЖЕНИЕ I
В I0 была реализована поддержка дискретных систем размерности 2. Вычисление арифметических выражений выполня-
лось традиционно - путем обхода дерева разбора. В реализацию I была добавлена поддержка непрерывных систем размерности 2, а также дискретных и непрерывных систем размерности 3. При организации процесса обработки непрерывных систем в I был реализован интервальный метод Рунге-Кутта 4-го порядка, разработанный Ю. Шокиным [7], и использованы приведенные оценки ширины интервала-решения, полученного данным методом.
Использование модуля динамической генерации кода позволило существенно оптимизировать процесс локализации инвариантных множеств.
В I, так же как и в I0 для хранения ячеек покрытия используется структура данных множество^ЕТ), а операции поиска, вставки и удаления производятся за время пропорциональное log N, где N - количество элементов в множестве. Моделью данных, соответствующей SET, является двунаправленный список, все элементы которого различны.
Так как операциями, существенно влияющими на скорость работы I0, являются вычисление образа ячейки (Img) и пересечение образа ячейки с покрытием (Intersect), то на их выполнение тратится основная часть общего времени работы программы.
При реализации распределенных алгоритмов решались следующие задачи:
1) распараллеливание операций Img и Intersect;
2) синхронизация доступа к ячейкам текущего покрытия (каждая ячейка должна обрабатываться однократно, например, если операция Img для одной ячейки будет выполняться в нескольких потоках одновременно, тогда быстродействие обработки ДС может существенно снизиться);
3) синхронизация доступа к новому покрытию (для корректной модификации нового покрытия поток-обработчик должен единолично владеть этим ресурсом).
4.2. ПРИЛОЖЕНИЕ II
Так как операция Intersect сводится к перебору ячеек покрытия, то предлагае-
мая оптимизация ориентирована на выявление способов организации представления покрытия таким образом, чтобы поиск ячейки становился наиболее эффективным. Мы используем подход, предложенный в работе [21]. Он основан на понятии Я-деревьев. В этом случае фазовое пространство логически представляется как последовательность иерархически вложенных многомерных кубов, а процесс поиска ячейки упрощается за счет выявления и отсечения тех областей покрытия, в которых ячейки быть не может.
В отличие от задач, связанных с использованием баз данных, нам достаточно хранить только координаты ячейки, поэтому потребность в хранении индексов многомерных объектов отпадает. Так как во всех реализациях этого алгоритма рассматриваются ячейки одного размера и используется целочисленная система координат, единица длины в которой есть размер ячейки, нам кажется естественным представление индексной записи в виде массива целых чисел - многомерного индекса ячейки, а листового узла - в виде упорядоченного списка таких кортежей, элементы которых различны. Не листовой узел в нашем случае представлен в виде списка пар: указатель на дочернюю вершину и целочисленный интервальный вектор. Если дочерняя вершина - листовой узел, тогда ограничивающий многомерный куб равен минимальному интер-
Рис. 1. Инвариантное множество отображения Ван-Дер-Поля
вальному вектору, содержащему все номера ячеек этого листа. Иначе ограничивающий интервальный вектор равен минимальному интервальному вектору, содержащему все интервальные вектора из списка пар дочернего узла.
4.3. ПРИЛОЖЕНИЕ III С ПОДДЕРЖКОЙ ОПТИМАЛЬНОГО ХРАНЕНИЯ ЯЧЕЕК
Каждой ячейке сопоставляется уникальное натуральное число - индекс. Ячейки можно упорядочить по возрастанию согласно их порядковому номеру. Совокупность индексов ячеек представлена в виде одномерного массива машинных слов, длина которого вычисляется по
следующей формуле: [count + 8], где count - количество ячеек в покрытии.
При таком представлении ячейка покрытия может быть представлена 1 битом : значение 1 - ячейку нужно обработать. Операции пересечения и добавления осуществляются за время пропорциональное 0(1). В обоих случаях производится операция обращения к элементу массива по индексу (то есть адрес элемента с индексом index равен адресу массива + index * размер элемента массива (в байтах)). Для того, чтобы определить, принадлежит ли ячейка покрытию, необходимо вычислить 2 величины: индекс слова, в котором ячейка могла бы располагаться, а также бит в слове, который соответствует этой ячейке. Обозначим за index порядковый номер некоторой ячейки. Тогда индекс слова вычисляется по формуле: word_index = [index + 8]. Позиция в слове: pos_index = index mod 8, проверка вхождения ячейки в покрытие осуществляется следующим образом: (1 << pos_index) & coverage_array[word_index]. Операция сдвига ( << ) позволяет эффективно возводить 2 в степень pos_index. У полученного числа разряды, не стоящие на позицииpos_index, имеют нулевое значение. Результат операции & равен нулю, если у coverage_array[word_index] на позиции pos_index нулевой бит.
В этой реализации был использован упрощенный вариант распределенного
алгоритма локализации инвариантных множеств. Мютех используется только при выборе ячейки для обработки. Стек обработанных ячеек и поток-наполнитель отсутствуют. При таком подходе удалось избежать снижения производительности, связанного с наличием взаимных блокировок при наполнении нового покрытия. При потребности хранения всех ячеек покрытия (даже тех, которые не участвуют в процессе локализации) адресное пространство процесса существенно сократилось за счет оптимального представления ячейки в памяти (1 бит). При выборе ячейки был использован алгоритм с отсечением, основанный на специфике представления хранилища индексов ячеек: если машинное слово равно нулю, значит, данный блок не хранит ячеек и, значит, перебирать его разряды не нужно. Осуществляется переход к следующему слову.
4.4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ
Все вычисления использовали следующую программную и аппаратную конфигурацию:
- CPU: core 2 DUO: E6420 2.13 Ghz,
- Memory 2 Gb,
- OS: Microsoft Windows XP SP 2.
4.4.1. Система Ван-дер-Поля
Рассмотрим следующую систему: x = у
у = 1.5(1 - x2)у - x.
Инвариантное множество отображения Ван-Дер-Поля содержит точку покоя (0, 0) и периодическую орбиту. Вычисления проводились в области
M =[-3.5, 3.5] х [-3.5, 3.5]. Функция системы для построения представления символического образа задается при помощи 100 итераций метода Рунге-Кутта с шагом 0.001 (оператор сдвига на 0.1) (рис. 1).
Рис. 2. Инвариантное множество системы О. Юнге
Начальная область М = [-0.38, 0.98] х [0.05,1.45] х [-0.38,0.98]. Начальное разбиение - [100,100]. Количество подразбиений -10. Результат показан на рис. 2.
4.5. СРАВНИТЕЛЬНЫЕ ХАРАКТЕРИСТИКИ ПРОИЗВОДИТЕЛЬНОСТИ ПРИЛОЖЕНИЙ I, II, III
Быши рассмотрены следующие системы.
Отображение с задержкой
ХЛ
= У,
у1 = 2.27 ау (1 - х).
Отображение Икеда х1 = 2 - 0.9(х cosт(x, у) - у sinт(x, у)), у1 = 0.9(у cosт(x, у) + х sinт(x, у)),
т(х, у) = 0.4 - --6-г.
1 + X + у
Система Дуффинга X = у ,
у = х - х3 - 0.15у .
Для системы Дуффинга вычисление образа ячейки производилось по методу Рунге-Кутта с шагом г, равным 0.03. Количество итераций - 100.
В приведенных далее табл. 1 и диаграмме (рис. 3) используются следующие обозначения:
Табл. 1. Сравнительные характеристики работы различных реализаций алгоритма
4.4.2. Система О. Юнге
x1 = y - mx,
У1 = Ау(1 - ^ z1 = x - gz, где l = 2.35, m = 0.5, g =0.1.
Система I II III
A 61.875 sec 48.045 sec 18.047 sec
B 237.192 sec 198.437 sec 71.609 sec
C 135.501 sec 107.782 sec 43.641 sec
D 540.328 sec 429.593 sec 166.767 sec
E 386.615 sec 350.890 sec 191.822 sec
F 1417.016 sec 1263.472 sec 593.032 sec
А - отображение с задержкой (в разбиении 4 миллиона ячеек),
В - отображение с задержкой (в разбиении 16 миллионов ячеек),
С - отображение Икеда (в разбиении 4 миллиона ячеек),
Э - отображение Икеда (в разбиении 16 миллионов ячеек),
Е - система Дуффинга (в разбиении 4 миллиона ячеек),
Б - система Дуффинга (в разбиении 16 миллионов ячеек).
Рис. 3. Диаграмма сравнительных характеристик работы различных реализаций алгоритма
Литература
5. ЗАКЛЮЧЕНИЕ
В работе описаны реализации алгоритмов локализации инвариантных множеств динамических систем, основанные на методах интервальной арифметики и символического образа. Основное внимание уделено способам оптимизации, а именно представлению данных и оптимизации вычислений. Реализованы версии алгоритма с различными оптимизациями и проведено сравнение их эффективности. Результаты численных экспериментов показали, что оптимальным является представление ячейки одним битом, означающим ее присутствие в текущем покрытии. При оптимизации вычислений используется динамическая генерация кода. Это позволяет отказаться от использования дерева разбора арифметического выражения при вычислениях. Описанные алгоритмы могут быть применены как в учебном процессе при моделировании динамических систем, так и в исследовательской работе.
1. Dellnitz M., Junge O. Set oriented numerical methods for dynamical systems. Handbook of dynamical systems. Vol. 2. Ed. B. Fiedler, 2002.
2. Осипенко Г.С. О символическом образе динамической системы. Сб. Граничные задачи. Пермь, 1983. С. 101-105.
3. Осипенко Г.С., Ампилова Н.Б. Введение в символический анализ динамических систем. СПб.: Изд. СПбГУ, 2005.
4. Ахо А., Сети Р., Ульман Дж. Компиляторы: принципы, технологии, инструменты. М.: Издательский дом «Вильямс», 2003.
5. Каток А.Б. и Хассельблат Б. Введение в современную теорию динамических систем // М.: Факториал, 1999.
6. МеньшиковГ.Г. Интервальный анализ и методы вычислений. СПб.: НИИХ СПбГУ, 19982001. Вып. 1. Введение в интервальную организацию вычислений // www. apmath. spbu. ru/ru/ education/courses/elective/menshikov.html.
7. Шокин Ю.И. Интервальный анализ. Новосибирск: Сибирское отделение изд-ва «Наука», 1981. С. 64-68 // http://www.nsc.ru/interval/index.php?j=Library/InteBooks/index.
8. Шарый С.П. Конечномерный интервальный анализ, «XYZ», 2009.
9. Терентьев C.B. Разработка и реализация основанных на интервальной арифметике алгоритмов компьютерного исследования динамических систем / Дисс. на соиск. к.ф.м.н, СПб.: СПбГУ, 2010.
10. Dellnitz M. and Hohmann A. A Subdivision Algorithm for the Computation of Unstable Manifolds and Global Attractors // Numerische Mathematik, 1997. Vol. 75, № 3. P. 293-317.
11. Ампилова Н.Б. Терентьев C.B. Применение методов интервальной арифметики к задаче построения символического образа // Электронный журнал «Дифференциальные уравнения и процессы управления», 2006. № 4. http://www.neva.ru/journal/j/index.html.
12. Dellnitz M., Junge O. An adaptive subdivision technique for the approximation of attractors and invariant measures // Comput. Visual. Sci., 1998. № 1. P. 63-68.
13. Galias Z. Rigorous investigation of the Ikeda map by means of interval arithmetic // Nonlinearity 15(2002). P. 1759-1779.
14. Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system, Opt. Commun. 30 (1979). P. 257-261.
15. Jaulin L., Kieffer M., Didrit O., Walter E. Applied interval analysis. Berlin: Springer, 2001.
16. Mischaikow K. and Mrozek M. Chaos in the Lorenz equations: A computer assisted proof. Part II: Details // Mathematics of Computation, 1998, July. Vol. 67, № 223. P. 1023-1046.
17. Neumaier A. Interval Methods for systems of equations, Cambridge University Press, 1990.
18. Tucker W. A rigorous ODE solver and Smale's 14th problem. Found. Comput. Math., 2002. № 2. P. 53-117.
19. .NET Framework Developer's Guide. Reflection. Microsoft Developers Network // http:// msdn.microsoft.com/en-us/library/cxz4wk15.aspx.
20. Kathleen Dollard. Code Generation in Microsoft .NET. Apress, 2004. P. 760.
21. Guttman A. R-Trees: A Dynamic Index Structure for Spatial Searching. Proc. 1984 ACM SIGMOD International Conference on Management of Data. P. 47-57 // http://www.sai.msu.su/ ~megera/postgres/gist/papers/gutman-rtree.pdf.
22. MSDN ( http://msdn.microsoft.com/library/).
23. Standard Template Library (http://msdn2.microsoft.com/en-us/library/c191tb28(VS.80).aspx).
24. BOOST ( http://www.boost.org).
25. Standard Template Library (http://gnuplot.sourceforge.net/).
26. MinGW (http://www.mingw.org/).
27. Antlr (http://www.antlr.org/).
Abstract
The paper is devoted to the design and implementation of the algorithm for localizing invariant sets of dynamical systems. The algorithm is based on interval arithmetics methods. The system is approximated by using a symbolic image which is an oriented graph constructed in accordance with the system. The cells of a covering are supposed to be interval vectors. The degree of approximation to the invariant set is estimated through the symbolic image parameters. The results of numerical experiments and the comparison with the algorithms based on real arithmetics are given.
Keywords: dynamical systems, invariant sets, symbolic image, computer simulation, interval arithmetics.
Ампилова Наталья Борисовна, кандидат физико-математических наук, доцент кафедры информатики математико-механического факультета СПбГУ, nataly@is1483.spb.edu,
Терентьев Сергей Валерьевич, кандидат физико-математических наук, инженер-программист Peter-Service, waterq@yandex.ru.
© Наши авторы, 2011. Our authors, 2011.
