CC BY
69
Сер. 10. 2009. Вып. 4
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
ПРИКЛАДНАЯ МАТЕМАТИКА
УДК 517.9
Н. Б. Ампилова, С. В. Терентьев
О ПРИМЕНЕНИИ ИНТЕРВАЛЬНОЙ АРИФМЕТИКИ
ПРИ ЧИСЛЕННОМ ИССЛЕДОВАНИИ ДИНАМИЧЕСКИХ СИСТЕМ
1. Введение. За последние годы динамические системы стали одним из основных инструментов моделирования многих реальных процессов. Практика показала, что методы компьютерного моделирования систем со сложным поведением траекторий являются неотъемлемой частью исследования.
Среди численных методов наиболее известны так называемые методы, основанные на множествах (set oriented methods) [1]. В их основе лежит идея об аппроксимации фазового пространства конечным набором клеток (ячеек). Построение образов ячеек дает приближенное представление о поведении системы. Для получения более точной информации применяется процесс последовательного подразбиения. При стремлении диаметров ячеек к нулю и выполнении разумных ограничений на систему такой процесс сходится и позволяет построить фазовый портрет системы с нужной точностью.
К довольно перспективным относится предложенный в 1983 г. Г. С. Осипенко [2] метод изучения поведения динамических систем с помощью символического образа, представляющего собой ориентированный граф, построенный по системе и заданному покрытию. Траекториям системы сопоставляются пути на графе, при этом одному пути может соответствовать несколько траекторий, вообще говоря. Символический образ — конечная аппроксимация исходной системы, при этом более мелкое покрытие порождает более точную аппроксимацию.
В последние годы в вычислительной математике широкое распространение получил метод интервального анализа. Методы численного решения задач с помощью интервальной арифметики подробно описаны в работах [3-5]. В работе [6] с целью обеспечения доказательных вычислений на ЭВМ был разработан и реализован алгоритм аппроксимации вещественных функций линейными сплайнами, на основании которого был создан пакет функционально-интервальной арифметики, позволивший решить
Ампилова Наталья Борисовна — кандидат физико-математических наук, доцент кафедры информатики математико-механического факультета Санкт-Петербургского государственного университета. Количество опубликованных работ: 48. Научные направления: компьютерное моделирование динамических систем, символическая динамика. E-mail: nataly@is1483.spb.edu.
Терентьев Сергей Валерьевич — аспирант кафедры информатики математико-механического факультета Санкт-Петербургского государственного университета. Научный руководитель: кандидат физико-математических наук, доц. Н. Б. Ампилова. Количество опубликованных работ: 6. Научные направления: компьютерное моделирование динамических систем, символическая динамика. E-mail: waterq@.yandex.ru.
© Н. Б. Ампилова, С. В. Терентьев, 2009
ряд вычислительных задач. В настоящее время известен ряд пакетов интервальных вычислений: XSC (extensions for scientific computations) расширения языков C, PASCAL, FORTRAN; SUN STUDIO - реализация интервальной арифметики для языков C++, FORTRAN; GAOL; BOOST [7].
В данной работе интервальные методы применяются для численного исследования динамических систем. Ячейки покрытия представляются как интервальные вектора, а их образы строятся с помощью интервального расширения функций, описывающих исследуемую динамическую систему. При таком подходе можно легко получать оценки параметров символического образа, гарантирующие определенную точность построения. Описанный метод применяется к задаче локализации инвариантных множеств и задаче о построении е-траектории, проходящей через заданные точки.
2. Символический образ. Символический образ строится по заданному покрытию фазового пространства и системе и представляет собой ориентированный граф, вершины которого соответствуют ячейкам покрытия. Для динамической системы f, покрытия {Mi}, i = 1,...,п, на символическом образе G строится дуга (i,j), если
f (Mi) n Mj = 9.
^мволический образ характеризуется следующими параметрами:
- максимальный диаметр ячейки d;
- нижняя грань r - наименьшее из расстояний между образами ячеек и теми
ячейками, с которыми образы не пересекаются;
- максимальный из диаметров образов ячеек q.
Определение 1 [8]. Бесконечная в обе стороны последовательность точек {xi, —ж < i < +ж} называется е-траекторией системы f (или псевдотраекторией), если для любого i выполняется неравенство dist(f (xi),xi+i) < е. Если последовательность {xi} является периодической, то она называется е-периодической траекторией. Нетрудно заметить, что любая точная траектория динамической системы является е-траекторией для Уе > 0. В работе [8] показано, что пути на графе G соответствуют е-траекториям, где величина е может быть оценена через параметры символического образа. Для уточнения приближений к траекториям системы используется процедура последовательного подразбиения и строится последовательность символических образов. Точность приближения определяется величинами d и г, которые стремятся к нулю при последовательном подразбиении.
3. Методы интервальной арифметики.
Определение 2. Интервалом называется множество вида х = [х,х\ =
{х G R | х ^ х ^ х} . Тогда rad(x) = 2^® называется радиусом интервала х.
Пространство всех интервалов над R обозначается IR.
Определение 3. Пусть xi,..., xm G IR. Тогда X = (xi,..., xm) называется интервальным вектором. Пространство интервальных векторов размерности m обозначается IRm.
Радиус интервального вектора определяется как
rad(x) = (rad(xi),..., rad(xm)), xi,..., xm G IR.
Нетрудно заметить, что интервальный вектор имеет довольно естественное геометрическое представление: прямое произведение составляющих его интервалов.
Определение 4. Пусть а, b С IR. Тогда функция р{а,Ъ) = max {|а — Ъ\,\а — Ц } задает расстояние между интервалами а и b. Для a,b G IRm расстояние есть максимум из расстояний по компонентам.
Так введенное расстояние является метрикой на соответствующих пространствах. Пусть / : Ят ^ Ят, а € 1Ят и V(/, й) = {/(х), х € а} есть область значений / на а.
Для вещественнозначной функции / рассмотрим ее интервальное расширение /(а), которое получается, если вычислять / по правилам интервальной арифметики. По известному свойству монотонности по включению верно следующее: V(/, й) С /(й). Разобьем й на более мелкие части а1 ,...,ап. Тогда из указанного свойства следует V(/, й) С уП=о /(аг) С /(й). Поскольку при стремлении диаметра подразбиения к нулю р(У(/, й),/(й)) оценивается через диаметр а (как множества), а для дифференцируемой функции / - через квадрат диаметра [9], то при последовательном подразбиении интервальных векторов получаем оценку для V(/, й) с любой заданной точностью.
Таким образом, интерпретируя ячейки покрытия как интервальные вектора в пространстве соответствующей размерности, можно рассматривать интервальное расширение исходной функции / как «отображение символического образа» и строить образ ячейки, используя правила и функции интервальной арифметики.
4. Локализация инвариантных множеств. Пусть {Ыг} есть покрытие фазо-
вого пространства Ы динамической системы /, где г = 1,...,п. Обозначим через Ыг интервальный вектор, соответствующий ячейке Ыг, тогда Ыг = а\ х ... х агт, где а* - ]-я компонента интервального вектора. Таким образом, ячейка
покрытия является многомерным параллелепипедом. Введем также обозначения: С = |м,-\Ыу П /(Ыг) = 0} , с* = ^\Ыу П /(Ыг) = 0} и Я* = Ujecí Ыу. Определим диаметр ячейки ё1аш(Ыг) = шахке[1т]{2гаё(Ыг)}к и обозначим через d наибольший из таких диаметров. Пусть г = штгкdist(/(Ыг),Ык), где /(Ыг) П Ык = 0 и dist обозначает хаусдорфово расстояние между множествами. (Отметим, что в данном случае удобнее использовать именно это расстояние, а не интервальное, так как последнее дает большее расстояние для непересекающихся интервалов.) Таким образом, для оценки параметров символического образа нужно вычислять d и г.
Пусть = |Ыь...,Ып} - начальное покрытие компакта Ы ячейками Ыг С
1Ят,г = 1,...,п. Для каждой ячейки Ыг вычислим множества Яг. Построим
множество-представление графа 61, соответствующее символическому образу 61 отображения / относительно покрытия &1, а именно Оог = и*е[1 п] йг. Это множество содержит все ячейки покрытия, которые участвуют в построении графа 61. Вычисляем значения dl и г1. На следующем шаге разбиение применяется к элементам множества Оог, в результате повторения описанной процедуры получается множество-представление 0^2 для графа О^. Вычисляем с],2,г2. Траектории локализуются с точностью £ > d2.
На к-м шаге рассматриваем разбиение вида
& к = |ЫЛ ,..., ЫУт \ ЫУ1 С ,
строим множество-представление 0^к и вычисляем соответствующие и г к. Таким образом, при стремлении диаметра разбиения к нулю локализуем инвариантные множества с точностью £ > dk.
5. Построение псевдотраектории.
5.1. Постановка задачи. Пусть задан гомеоморфизм / : Ы ^ Ы, Ы С Яп и точки х*,х**. Определить, с каким £ можно построить приближенную траекторию, проходящую через точку х* вблизи точки х**. (Предполагается, что £ достаточно мало.)
Построим траекторию точки х*, т. е. последовательность хк такую, что хк+1 = /(хк), хо = х*, к = 1, 2,.... Если существует такое число N, что /(х^) = х**, то мы построили точную траекторию точки х*, проходящую через точку х**. Она является £-траекторией для У£ > 0.
Если же существует такое число N, что р(/(х^),х**) < £, то построена точная траектория точки х*, проходящая вблизи точки х**. Она является £-траекторией с заданным £, в силу последнего неравенства.
Таким образом, задача имеет очевидное решение при малых значениях £: отрезок точной траектории, который либо соединяет заданные точки, либо проходит вблизи второй точки.
Кроме того, могут существовать другие решения задачи. Для этого рассмотрим полученную траекторию хо,х1,...,х^ ,х** и множества и (хк )-£-окрестности точек хк, к = 0,...,N — 1. (Будем рассматривать такие окрестности как интервальные вектора с диаметром £.) Составим следующую последовательность: хо = хо,х1 € и(х1),. ..,х^ € и(х^), х** и оценим 6, при котором она будет некоторой ¿-траекторией. Заметим, что по выбору £ точка хN лежит в £-окрестности точки х**.
Так как \/(хо) — х1 \ < £, т. е. /(хо) € и(х1), то точки хо,х1 принадлежат £-траектории. Следовательно, для того чтобы оценить значение 6, нужно оценить величину \/(хк) — хк+1 \, к = 2,...^ — 1.
Обозначим як = ё1аш/(и(хк)).
Возможны следующие ситуации:
а) / (хк) € / (и (хк)) П и (хк+1), к € [2, N — 1], тогда \/(хк) — хк+1 \ < £ + Як, где хк € / (и (хк )),хк+1 € и (хк+1).
Положим я = шахк як, к € [2^ — 1], тогда траектория {хк} - (£ + я)-траектория системы /;
б) /(хк) € /(и(хк)) П и(хк+1), к € [2^ — 1]. Тогда траектория {хк} есть
£-траектория системы /.
5.2. Использование константы Липшица. Как показано в [10], если / удовлетворяет условию Липшица в х € 1Ят, то справедливо неравенство гаё(/(х)) ^
Xf (х)гаё(х), где Xf - константа Липшица, которая вычисляется через константы Липшица для функций, входящих в /.
Таким образом, Аf может быть вычислена в процессе анализа правых частей системы и диаметр образа ячейки естественным образом оценивается через диаметр ячейки
и А1.
5.3. Описание алгоритма.
1. Строим траекторию точки х*, вычисляя при этом £г = ё1в1(х**, /г(х*)), г = 1,...,к.
2. Выберем £ = ш1пг{£г} и обозначим через т индекс, при котором этот минимум достигается. Если т < к, то количество шагов уменьшается, значение параметра т становится новым значением количества шагов.
3. Введем множество '^гкзе!. В начале работы "мэгкве! = и (х*). Для г € {1,. ..,т — 1} это множество переопределяется по правилу "йюгкве! := "йюгк8е1и(/(и(хг)) р| и(хг+1))). В конце работы workset содержит последовательность интервалов и!г € 1Яп, г € {0,1,...,т}. Любая траектория Е, построенная по правилу Е = {хо, х1,..., хт} , хг € и!г, является £-траекторией.
6. Реализация.
6.1. Смешанные вычисления. Этот термин был введен для обозначения техники оптимизации компьютерных программ. Обозначим через Рд программу, которая принимает на вход любые данные из множества данных I и на выходе выдает
данные из набора данных O, а через Pi - программу, которая принимает на вход любые данные из Ii С I, имеет существенно более простую по сравнению с Pg структуру и выдает данные из O\ С O. Причем на схожих входных данных программы выдают похожие результаты. Тогда скорость обработки данных Ii программой Pi будет выше, чем скорость обработки данных Ii программой Pg. Если сделать предположение о том, что на входе программы Pg всегда будут данные из множества Ii, то для их обработки можно использовать Pi [11]. (Подобную технику оптимизации компьютерных программ использовали В. Ф. Турчин, С. А. Романенко [12] при создании РЕФАЛ-компилятора.)
В рассматриваемом случае программой Pg является программа, которая может работать с дискретной либо непрерывной системой порядка 2 или 3, заданной на произвольном фазовом пространстве. Программа Pi работает только с конкретной системой, т. е. в ее коде явно заданы размерность и обработка арифметических выражений, описывающих систему. Был разработан и реализован компилятор, позволяющий получить Pi на основе описания динамической системы, заданной в конфигурационном файле программы.
6.2. Алгоритм локализации инвариантного множества. Такой алгоритм реализован на языке C++ с использованием технологии MS Visual Studio 2005 [13, 14]. Благодаря внедрению поддержки смешанных вычислений, производительность реализованного инструмента существенно возросла. Визуализация результата реализована с помощью технологии GNUPLOT [15]. При реализации алгоритма была использована реализация интервальной арифметики из библиотеки BOOST. Инструмент был разработан для платформы Windows. Смешанные вычисления реализованы с помощью технологии Antlr [16] и набора программных инструментов MinGW [17]. Предложенный инструмент работает с дискретными и непрерывными системами размерности 2, 3. При организации процесса обработки непрерывных систем был реализован интервальный метод 4-го порядка типа Рунге-Кутта, созданный Ю. И. Шокиным [3], и использованы приведенные оценки ширины интервала-решения, полученного данным методом.
6.2.1. Пользовательский интерфейс. Взаимодействие между пользователем и приложением может осуществляется либо посредством командной строки, либо путем ввода данных через специальное оконное приложение. Пользовательский графический интерфейс реализован на языке C# [13] с помощью технологий .Net Windows Forms [13], GNUPLOT [15].
6.2.2. Р е а л и з а ц и я смешанных вычислений. При реализации данного алгоритма была сформирована библиотека, содержащая в себе средства для объектной аппроксимации исходной системы в виде некоторого экземпляра динамически генерируемого класса, хранящего основные параметры системы (размерность и границы фазового пространства и т. д.) и позволяющего эффективно ими управлять [11]. При таком подходе из конфигурационного файла извлекается описание динамической системы (ДС) и подается на вход специальному анализатору [18], который строит абстрактное дерево разбора. Далее производится сопоставление между функциями и переменными, входящими в описание системы и их определениями, и строится семантическое дерево разбора. Затем специальный генератор [18] по построенному дереву разбора генерирует файл, содержащий код на языке С++ класса, являющегося абстрактной моделью исходной системы. После этого производятся компиляция и формирование библиотеки (DLL) посредством технологии MinGW [17]. Далее DLL загружается в адресное пространство запущенного инструмента.
6.3. Алгоритм поиска псевдотраекторий. Алгоритм поиска е-траекторий реализован на языке C#. С целью повышения производительности инструмента был
разработан распределенный алгоритм построения множества workset. Для эффективной работы с арифметическими выражениями в реализацию данного алгоритма была внедрена поддержка смешанных вычислений. Визуализация результата реализована с помощью GNUPLOT технологии. Пользовательский интерфейс создан с помощью технологии .Net Windows Forms [13].
Реализация смешанных вычислений. Механизм эффективного взаимодействия с арифметическими выражениями реализован с помощью технологий Antlr [16], .Net reflection [13] и .Net Code Document Object Model [13]. В процессе обработки конфигурационного файла создается класс на языке C#, который является абстрактной моделью исходной динамической системы. В отличие от инструмента локализации инвариантного множества, в данном случае сборка (Assembly [13]) производится сразу в адресном пространстве запущенного инструмента без создания физического DLL-файла^иблиотеки).
7. Численные эксперименты.
7.1. Построение инвариантного множества.
Пример 1. Рассмотрим отображение Икеда [19]
x1 = 2 — 0.9(x cos т(x,y) — y sin т(x,y)) ,
y1 = 0.9 (ycosт(x,y)+хsinт(x,y)) ,
тдет(х,у) = 0.4- 1+J+y2.
Вычисления проводились в области [—2, 4] х [—3, 2]. В таблице приведены значения параметров символических образов, построенных для последовательных шагов подразбиений. Полученное инвариантное множество показано на рис. 1.
У
1.5
Рис. 1. Инвариантное множество отображения Икеда
Параметры символических образов, полученных при последовательном подразбиении
Шаг Разбиение <і <? г
1 50x50 0.12 14.39699 0.0206512260410854
2 100x100 0.06 14.40033 0.0160618211918824
3 200x200 0.03 15.67012 0.00311592509026526
4 600x600 0.01 7.974839 0.00594671759220252
5 1200x1200 0.005 7.36462932 0.00353825551610235
Пример 2. Рассмотрим систему О. Юнге
Х1 = у — цх,
У1 = Лу(1 — х),
г\ = х — ^г,
где Л = 2.35; р = 0.5; 7 = 0.1.
Начальная область М = [—0.38, 0.98] х [0.05,1.45] х [—0.38, 0.98]. Начальное разбиение - [100,100]. Количество подразбиений - 10. Результат показан на рис. 2.
0.89
Рис. 2. Инвариантное множество отображения О. Юнге
Пример 3. Рассмотрим систему Дуффинга
х = У,
у = х — х3 — 0.15у.
При данных значениях параметров у системы существуют два асимтотически устойчивых состояния равновесия ( —1,0) и (1,0). Точка (0,0) является фокусом. Начальное разбиение - [100,100]. Количество подразбиений - 30. Вычисления производились по методу Рунге-Кутта с шагом по Ь, равным 0.03. Количество итераций - 100. На рис. 3, а в верхнем ряду показано положительно инвариантное множество в области М = [—1.5,1.5] х [—2, 2]. На рис. 3, б построено положительно инвариантное множество в окрестности фокуса ( — 1,0). На рис. 3, в построены положительно инвариантные
множества в окрестностях особых точек, в каждой окрестности расчеты производились независимо; на рис. 3, г - отрицательно инвариантные множества в окрестностях состояний равновесия.
1.5 1.5 -1.5 1.5
Рис. 3. Инвариантное множество уравнения Дуффинга в зависимости от начальной области
Пояснения в тексте.
7.2. Построение £-траектории. Рассмотрим отображение Хенона /(х,у) = (1 + у — 1.4х2х, 0.3х). Проверим существование псевдотраектории для начальных данных:
хо = (0.631254,0.189406); х** = (0.660379015,0.287993948); к = 6.
Отрезок точной траектории, проходящий через точку х0:
1. (0.631254, 0.189406);
2. (0.6315317, 0.1893762); £1 = 0.1027503;
3. (0.6310109, 0.1894595); £2 = 0.1028179;
4. (0.6320148, 0.1893033); £3 = 0.1026858;
5. (0.6300834, 0.1896044); £4 = 0.1029481;
6. (0.6337973, 0.1890250); £5 = 0.1024765.
Возьмем наименьшее £ = £5 и построим последовательность интервалов, следуя алгоритму:
1. ( [ 0.621006348297279, 0.641501651702721 ], [ 0.179158348297279, 0.199653651702721 ] );
2. ( [ 0.621284090774879, 0.641779394180321 ], [ 0.186301904489184, 0.192450495510816 ] );
3. ( [ 0.620763269837767, 0.64125857324321 ], [ 0.186385227232464, 0.192533818254096 ] );
4. ( [ 0.621767174695831, 0.642262478101274 ], [ 0.18622898095133, 0.192377571972963 ] );
5. ( [ 0.619835787656796, 0.640331091062239 ], [ 0.186530152408749, 0.192678743430382 ] );
6. ( [ 0.660379015, 0.660379015 ], [ 0.287993948, 0.287993948 ] ).
Выбирая в указанных интервалах произвольные точки, можно получить е-траектории.
Если выбрать в качестве начальных данных координаты седловых неподвижных точек (0.61313544, 0.18940634) и ( — 1.1313545, —0.3394063), то при этих условиях не существует приближенной траектории, проходящей через заданные точки.
8. Заключение. В работе показано применение методов интервальной арифметики к численному исследованию динамических систем. Используется метод аппроксимации системы с помощью символического образа, представляющего собой ориентированный граф, построенный по системе. Ячейки разбиения рассматриваются как интервальные вектора в пространстве соответствующей размерности. Разработан и реализован метод локализации инвариантных множеств с оценкой параметров символического образа, позволяющих определить точность построения. Также реализован способ построения псевдотраектории системы, проходящей через заданные точки. Приведены результаты численных экспериментов.
Литература
1. Dellnitz M., Junge O. Set oriented numerical methods for dynamical systems // Handbook of dynamical systems / ed. by B. Fiedler. Stockholm: Elsiver, 2002. Vol. 2. P. 223—235.
2. Осипенко Г. С. О символическом образе динамической системы // Граничные задачи / под ред. В. А. Алексеева. Пермь, 1983. С. 101—105.
3. Шокин Ю. И. Интервальный анализ. Новосибирск: Наука, 1981. 111 с.
4. Шарый С. П. Конечномерный интервальный анализ. Новосибирск: XYZ, 2009. 572 с.
5. Меньшиков Г. Г. Интервальный анализ и методы вычислений. СПб.: Науч.-исслед. ин-т химии С.-Петерб. ун-та, 1998—2001. Вып. 1. Введение в интервальную организацию вычислений. URL: www.apmath.spbu.ru/ru/education/courses/elective / menshikov.html.
6. Югай С. А. Разработка и применение функционально-интервального анализа для обеспечения доказательных вычислений на ЭВМ: дис. на соискание учен. степени канд. физ.-мат. наук. Фрунзе: Ин-т математики АН Кирг. ССР, 1988. 97 с.
7. BOOSTREF BOOST. URL: http://www.boost.org.
8. Осипенко Г. С., Ампилова Н. Б. Введение в символический анализ динамических систем. СПб.: Изд-во С.-Петерб. ун-та, 2005. 236 c.
9. Tucker W. A rigorous ODE solver and Smale’s 14th problem // Found. Comput. Math. 2002. Vol. 2. P. 53-117.
10. Neumaier A. Interval Methods for systems of equations. Cambridge: Cambridge University Press, 1990. 245 p.
11. Ershov А. Mixed calculation // Journal: In scientific world. 14.02.1984. URL: http://ers-hov.iis.nsk.su/archive/eaindex.asp?did=2596.
12. Романенко С. А., Турчин В. Ф. Рефал-компилятор // Труды 2-й Всесоюз. конференции по программированию. Заседание Б. Новосибирск, 1970, 3-6 февраля. URL: http://www.refal.org/ori-gins/RfcVkp2/index.htm.
13. MSDNREF MSDN. URL: http://msdn.microsoft.com/library/.
14. Standard Template Library. URL: http://msdn2.microsoft.com/en-us/library/c191tb28(VS.80).aspx.
15. GNUPLOTREF Standard Template Library. URL: http://gnuplot.sourceforge.net/.
16. ANTLRREF Antlr. URL: http://www.antlr.org/.
17. MINGWREF MinGW. URL: http://www.mingw.org/.
18. Ахо А., Сети Р., Ульман Дж. Компиляторы: принципы, технологии, инструменты. М.: Изд. дом Вильямс, 2003. 768 с.
19. Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system // Opt. Commun. 1979. Vol. 30. P. 257-261.
Статья рекомендована к печати член-кор. РАН, проф. Г. А. Леоновым.
Статья принята к печати 28 мая 2009 г.
