Алгоритм вычисления скользящего режима для системы с гладкой границей разрыва Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 15, № 2, 2010

Алгоритм вычисления скользящего режима для системы с гладкой границей разрыва

В. В. Коровиным. Ю.В. Фролова Омский государственный университет, Россия e-mail: [email protected]

Описан численный алгоритм, предназначенный для решения систем обыкновенных дифференциальных уравнений, фазовое пространство которых состоит из двух областей непрерывности поля направлений системы с общим участком гладкой границы, где реализуется скользящий режим. Особенность алгоритма состоит в точном определении точки пересечения траектории решения с границей областей непрерывности. Для вычисления решения в скользящем режиме строятся две траектории по разные стороны от границы. Результаты вычислительного эксперимента показали сходимость предлагаемого алгоритма. Получена экспериментальная оценка вычислительных затрат алгоритма.

Ключевые слова: дифференциальные уравнения с разрывной правой частью, методы Рунге — Кутты, скользящий режим.

Введение

Обыкновенные дифференциальные уравнения (ОДУ) с разрывной правой частью активно используются для описания математических моделей разного рода процессов и систем. Исследованию динамических систем с разрывной правой частью, в том числе со скольжением, посвящены работы [?-?]. Такие системы получили различные названия: гибридные, релейные, кусочно-сшитые, системы Филиппова, динамические системы со скольжением и с разрывами, системы с клеточной структурой. Наиболее полное теоретическое описание систем ОДУ с разрывной правой частью дано в [?].

К настоящему времени имеется много работ, посвященных различным аспектам реализации алгоритмов численного решения ОДУ с разрывной правой частью [?-?]. Однако, проведя анализ существующих реализаций численных методов решения ОДУ [?-?], нельзя считать, что существуют готовые реализации методов решения систем Филиппова в общем виде.

В данной статье предлагается численный алгоритм, предназначенный для решения систем ОДУ с фазовым пространством, состоящим из двух областей непрерывности ПОля направлений системы с общим участком гладкой границы, на которой реализуется скользящий режим. В отличие от предыдущих работ авторов [? - ?] здесь реализован алгоритм, имеющий три режима: непрерывный, пересечения границы и скольжения вдоль границы. Особенность алгоритма состоит в том, что для определения точки пересечения траектории решения с границей областей непрерывности правой части уравнения не используются значения функции правой части системы вне области непрерывности.

© ИВТ СО РАН, 2010.

Для вычисления значения функции правой части системы на границе (в скользящем режиме) используются значения функций по разные стороны от границы. Таким образом, вдоль границы движутся две траектории численного решения, которые являются приближением точного решения с двух сторон. После выхода со скольжения две траектории сливаются в одну, и решение системы продолжается в непрерывном режиме.

В работе проведен вычислительный эксперимент с реализованным алгоритмом на примере, имеющем гладкую границу со скользящим режимом. Результаты показали сходимость алгоритма к точному решению. Получена экспериментальная оценка вычислительных затрат алгоритма.

1. Определение системы ОДУ на поверхности разрыва

Будем рассматривать автономную систему обыкновенных дифференциальных уравнений с кусочно-непрерывным полем направлений. Для таких систем характерно разбиение пространства состояний на отдельные области, в которых векторное поле является гладким. Границы между этими областями называют поверхностями разрыва. В общем виде система ОДУ записывается как

Их

^ = I(х), х е (1)

где векторное поле I (х) может быть либо гладким, либо кусочно-гладким. Предположим, что пространство состояний состоит из двух областей Сг и Су, разделенных поверхностью разрыва , которая определяется гладкой скалярной функцией д^ (х)

Бгз = {х е : дгз(х) = 0},

а также

Сг = {х е Ега : дгу(х) > 0}, Су = {х е Ега : дгу(х) < 0}.

Тогда систему (??) можно переопределить как

Их = ( Рг(х), х е Сг, Иг = 1 Р(х), х е Су,

(2)

где Рг и Ру — достаточно гладкие функции. Будем считать, что эти функции определены в соответствующих областях Сг и Су, а также на поверхности Бгу. При этом допускается, что функция Рг может быть не определена в области Су и, наоборот, функция Ру может быть не определена в области С г.

Если векторные поля Рг и Ру вблизи поверхности разрыва Бгу одновременно направлены к ней или от нее, то решение будет двигаться вдоль этой поверхности. Такое решение называют скользящим режимом. Открытое подмножество Бгу поверхности Бгу, где векторные поля одновременно направлены к поверхности разрыва или от нее, обычно называют поверхностью скольжения. Кроме того, если справедливо неравенство

Ир-^ (дгу )(х) < 0, х е ¡гу, то поверхность разрыва Бгу является устойчивой, но если

(дц)(х) > 0, х е БИ,

то — неустойчива. Здесь через Ср (д)(х) обозначена производная Ли от функции д(х) вдоль векторного поля Г (х)

ср= |(х) = Ц(х> = {£ (х>Г(х>

На поверхности разрыва систему (??) можно доопределить как

— = (х), х ^ £>гз,

Fij {Х) = +m^iM ^ {Х),

(3)

2 2

а —1 < ¡¡^ (х) < 1. Поскольку траектория решения должна двигаться по поверхности , то векторное поле Г^ должно быть направлено по касательной к этой поверхности, т.е. (д^)(х) = 0. Отсюда

Vij(x) = -

LFi+Fj Í9ij)(x) LFj-Fi Í9ij)(x) '

(4)

В результате получаем систему вида

Fi(x), x е Gi

dt = 1 Fij (x), x E Sij F (x), x е Gj

dx dt

(5)

Такое доопределение системы (??) согласуется с подходом Филиппова, где система на поверхности разрыва заменяется дифференциальным включением с минимальным вы-

Fi Fj

гладкость функций Fij и g-, решение системы (??) будет единственным справа, кроме точек на Sj-, где оба вектopa Fi и Fj ортогональны к S-. Если оба вектора направлены к поверхности Sто эти точки являются стационарными точками системы, и решение будет единственным. Если же оба вектора направлены от поверхности разрыва, то SSij

ности^не обладает единственностью. Однако если решение не начинается с точки в этой области, то и траектория решения не попадет на нее. Далее решения, начинающиеся с точки на неустойчивой границе, рассматриваться не будут, поэтому все исследуемые решения системы (??) будут обладать единственностью cnpciBci.

2. Точное определение точки входа траектории решения на поверхность разрыва

Для точного вычисления решения со скользящим режимом необходимо точно определить точку входа траектории решения на поверхность разрыва. Но мы будем искать не одну точку, а две, которые находятся достаточно близко от поверхности, но по разные от нее стороны. Обе точки должны быть близки к точному положению траектории решения на поверхности разрыва.

Выше было сделано предположение о том, что функция Рг (х) может быть не определена в области Су, а функция Ру(х) не определена в

Сг

для вычисления

функции Ргу требуется определить значения функций Рг и Ру для одной и той же точки х, где х е Сг или х е Су.

Для точки х е Сг найдем близкую к ней точку х* е Су такую, что \х — х*\ < е\х\, где е — заданный предел локальной относительной погрешности решения. Такую

точку можно найти, когда расстояние от точки х до поверхности Бгу достаточно мало

е

(меньше -|х\). Тогда можно принять точку х* за представление точки х в области Су

и определить значения функций Ру(х) := Ру(х*) и, наоборот, Рг(х*) := Рг(х). х е Сг

будет определяться стандартной схемой численного решения ОДУ с непрерывной правой частью. Однако при приближении к поверхности разрыва для осуществления очередного шага интегрирования может потребоваться значение функции правой части

Су Рг Су

вопрос — как найти длину шага интегрирования и точку пересечения траектории с по-верхностыо разрыва. Предлагаем сделать это в два этапа: сначала определить длину

Сг

статочно близко к границе, а затем, используя экстраполяцию кривой решения, найти две точки хь е Сг и х* е Су, лежащие на продолжении траектории решения, но по разные стороны от поверхности разрыва. Для первого этапа выберем шаг

к = —К Л ), (6)

(дгу)(х)

где К > 0 — константа, регулирующая приближение к границе (чтобы наиденное зна~ чение решения осталось в области Сг, можно задать К = 0.9). Сделаем шаг интегрирования к го точки хк в точку хк+1.

Для экстраполяции кривой решения до границы воспользуемся интерполяционной формулой Ньютона для интерполирования назад, опираясь на три точки хк ~ х(гк — к), хк+1/2 ~ х(гк — к/2) и хк+1 ~ х(гк)■ Точка хк+1/2 будет получена при нахождении точки хк+1, если используется схема Ричардсона для оценки погрешности решения. В точках хк, хк+1/2, хк+г можно вычислить значения функции Рг: ¡к, ¡к+1/2,¡к+1, причем первые два значения были вычислены при нахождении точки хк+1. Обозначим также ¿1 = (хк+1 — хк+1/2)/к, И2 = (хк+1/2 — хк)/к. Тогда интерполяционный многочлен будет иметь

к+1

+

х(гк + вк) и N (в) = хк+1 + вк в-1 /, _ _ в-2

+ вЫ — ¡к+1 + (в — /к+1/2 — 2И1 + 1к+1+

4

(¿2 — 41к+1/2 + 5И1 — 2¡к+1 + • (¡к — 3И2 + 41к+1/2 — 3И1 + ¡к+1)))

(7)

Построенный многочлен имеет четвертый порядок точности.

Чтобы найти точки хь и х*, необходимо определить соответствующие вг и вг такие, что хь = N(вг) и х* = N(вг). Для поиска этих значений воспользуемся итерационным методом Ньютона для решения алгебраического уравнения

дгу (х) = 0, где х = N (в).

Итерационная процедура имеет вид

еа+1 = в3-е- ^(Хз\, ,, 0с = 0.5, ^ = о, 1,... к • (д-)(х8)

Здесь в качестве производной от функции д- используется производная по направлению СЕ1(д-)(х8), умноженная на к (производная сложной функции д—(х) • х'(0)). Значение точки х8 = х(Ьк+08к) вычисляется то формуле (??). Коэффициент С = 1.1 позволяет постоянно "перепрыгивать" через границу, обеспечивая приближение к ней с обеих сторон. Осуществляя итерационную процедуру, необходимо сохранять величину 01 = шах{08},

для которой х(Ьк + 08к) Е Ог, и 0Г = шт{08}, для которой х(Ьк + 08к) Е О-. Таким способом будем определять наиболее близкие точки по разные стороны от поверхности Б-, лежащие на продолжении кривой решения. Если в ходе итерационной процедуры будет выполнено условие |08| > 1, то необходимо остановить процедуру по причине расходимости метода. Успешным завершением процедуры является выполнение условия

101 - 0г I <

£

К • -

П

ГД6 К КОНСТаНТа (экспериментально установлено значение к = 0.0005, позволяющее найти точку пересечения с заданной точностью £, значение п определяется по формуле

П = т~~, гДе т = ^ , ^08 = 08 - 08-1, в = 2, 3,... 1 - т Д08-1

К

ходит за рамки данной статьи.

3. Условия продолжения решения на поверхности разрыва

Имея значения точек хь и х£, можно продолжить решение системы после достижения поверхности разрыва. При этом о дальнейшем поведении решения можно сделать вывод, анализируя значения вектор-функций Рг(хь) и Fj(хЬ) по отношению к поверхности Б-. Если выполнено условие пересечения

СРг(д-)(хь) -Се,(дг-)(хЬ) > 0, (8)

то траектория решения должна покинуть поверхность разрыва. При СЕ1 (д-)(хь) < 0 решение переидет в О г и продолжить его необходимо из точки хь, пр и СЕ1 (д- )(хь) > 0 решение переидет в О.,- и продолжить его необходимо из точки хЬ- Если эти векторы одновременно направлены к поверхности разрыва Б- или от нее:

СЕ1 (дг-)(хь) • СЕ, (дг-)(х*ь) < (9)

то реализуется скользящий режим. Причем если первый множитель СЕ1 (дг-)(хь) > 0 (векторы направлены к поверхности), то скользящий режим является устойчивым (поверхность притягивающей), и наоборот, при СЕ1 (д-)(хь) < 0 (векторы направлены от поверхности) скользящии режим является неустойчивым (поверхность отталкивающей). В последнем случае траектория решения может "соскользнуть" с поверхности

в любой момент, что не дает возможности говорить об единственности решения системы. Однако возможна ситуация, когда выполняется равенство

^ (9ц)(хь) • (9ц)(х*>) = 0. (10)

В этом случае, если один из сомножителей не равен нулю, то на поверхности разрыва одно из векторных полей направлено ортогонально к поверхности, а проекция другого вектора определяет скорость скольжения вдоль поверхности. Вместе с тем, если оба сомножителя равны нулю, то рассматриваются условия более высокого порядка — вторые производные функции дгу (х) [?].

4. Две скользящие траектории вблизи поверхности разрыва

Наиденные выше точки хь и Хь представляют собой пару, которая служит приближением точки решения х(гк+1) на поверхности разрыва Бгу, причем

+ _ + I 0 + гк+1 = гк +--^-

Для продолжения решения рассмотрим случай скользящего режима.

Предлагаем находить не одну траекторию решения, а две — по обе стороны ОТ ПОверхности разрыва. Одна траектория стартует из точки хь, вторая — из точки хЬ■ Траектории будем искать как решение системы

— = (х; хь), — = Еуг(х>; х), х Е Сг, хь Е Су, Ъ > Ък+1 (11)

= (х; х )' ~!г

с начальными данными

х(гк+1) = хь, х(гк+1) = х».

Функция Еу определяется формулами (??), (??) с заме ной аргумента х на х> для функции Ру\

Е(х) + Еу(х>) + Еу(х>) - Е(х) х>)

Е гу (х; х ) — п + О (х; х ),

^гу(х; х ) = I Ад \ / т \ , (12)

(х>),Р, (х>)) -I ^ (х),П(х)

^(х>),Еу(хьЛ -/

причем Егу(х; х>) = Еуг(х>; х), а это означает, что скорости движения вдоль поверхности разрыва точек траекторий с обеих сторон будут одинаковы. При этом х(г) Е Сг, х>(г) Е

Су ДЛЯ ЛК)бЫХ г ^ ¿к+1-

Движение траектории в скользящем режиме может прекратиться, и тогда траектория либо продолжит движение в непрерывном (не скользящем) режиме, либо достигнет устойчивой точки равновесия и останется в ней, либо продолжит движение в скользящем режиме, но по другой поверхности разрыва. Прекращение скользящего режима возможно при достижении траектории решения точки, в которой выполнимо одно из трех условий: 1 — векторы правых частей направлены по одну сторону от поверхности разрыва; 2 — точка является особой устойчивой; 3 — точка лежит на пересечении с другой поверхностью разрыва. В первом случае решение необходимо продолжить из

точки, лежащей по одну сторону от поверхности разрыва с вектором правой части, во втором — траектория решения останется в особой точке, в третьем — поведение решения будет зависеть от векторных полей по разные стороны от обеих поверхностей разрыва и в данной работе не рассматривается.

5. Описание алгоритма

Предлагаем алгоритм вычисления численного решения кусочно-сшитой системы (Piece-wise Sewing System) обыкновенных дифференциальных уравнений (??) с ТОЧНЫМ ВЫчислением скользящего режима.

Алгоритм. Обеспечивается вычисление решения кусочно-сшитой системы обыкновенных дифференциальных уравнений на интервале времени [a, b] для начального значения xq. Решение сохраняется в массивы Т и X. Используется схема Рунге — Кутты для интегрирования уравнения в области непрерывности. Для поиска решения на границе применяется метод экстраполяции полиномом Ньютона (??)■ Для нахождения решения в скользящем режиме используется метод Рунге — Кутты для системы уравнений (??), (??).

Шаг 0 (инициализация). Установим начальные значения x — Xq. t <— a i 0 Сохраним в массив Т[i] — t, X [i] — x, i — i + 1. Выберем начальный шаг интегрирования h. Установим режим mode — cont.

Шаг 1 (выбор режим,а, вычислений). Если mode = cont, то перейдем на шаг 2, если mode = intrs — на шаг 3, если mode = slide — на шаг 5.

h

Кутты. Если произошло пересечение поверхности разрыва, то выберем длину шага hi согласно формуле (??), немного не доходя до пересечения границы, выполним шаг h —— hi и переключим режим "mode —— intrs. Переидем на шаг 6.

Шаг 3 (режим пересечения). Найдем значение двух точек x^ и xl, лежащих на продолжении траектории решения: до и после пересечения границы, расстояние между точками должно быть меньше заданного предела погрешности е (см. разд. 2). Если алгоритм не сходится, то уменьшим шаг h — h/2, переключим режим mode — cont и перейдем на шаг 1, иначе — на шаг 4.

Шаг 4 (скольжение или продолжение?). Если условие скольжения (??) ВЫ ПОЛнено, то сохраним среднюю точку и обе точки используем для скользящего режима mode — slide. Если условие скольжения не выполнено, то сохраним обе точки и продолжим решение в области выхода траектории в непрерывном режиме mode — cont. Перейдем на шаг 6.

Шаг 5 (режим скольжения). Выполним шаг интегрирования методом Рунге — Кутты для скользящего режима. Если произошло пересечение границы области, то выберем длину шага h2 до пересечения границы, выпо л ним шаг h — h2 до границы и переключим режим mode — intrs. Если скользящий режим завершился, mode cont

Шаг 6 (сохранить точку решения и продолжить вычисления?). Если погрешность найденной точки решения меньше заданной точности, то сохраним результат Ь ^ Ь + к, Т [г] ^ Ь, X [г] ^ х, г ^ г + 1. Если Ь < Ь, то вернемся к шагу 1, иначе — закончим вычисления.

Решение системы (??) по представленному алгоритму определяется тремя режимами: непрерывный, пересечение и скольжение. Непрерывный режим реализуется одной из схем Рунге — Кутты. Например, можно использовать схемы Мерсона, Фельберга или Дормана — Принса [?], а также многошаговые методы, но в нашем вычислительном эксперименте была применена схема Рунге — Кутты — Фельберга.

х

в два этапа: определение номера г области Ог, содержащей х, и вычисление соответ-ствующби функции Fi(x). Определение номера области сводится к вычислению логических функций, определяющих условия ограничения пространства области. Однако при вычислении решения вблизи границы по схеме Рунге — Кутты может возникнуть ситуация, когда для одного шага по схеме необходимо вычислить значения функции правой части в разных областях Ог и О-. Это означает, что возможно пересечение границы областей траекторией. Поэтому необходимо более точно определить точку пересечения траектории решения с поверхностью разрыва, для чего реализован режим пересечения. Схема определения точек пересечения траектории решения с поверхностью разрыва приведена выше в разделе 2.

После нахождения точек пересечения необходимо определить дальнейшее поведение решения: продолжить решение в другой области или перейти в скользящий режим (условия определения поведения (??), (??))• Если выполнено условие пересечения, то обе точки хь и хЬ сохраняются в выходные массивы X и Т, а решение продолжает-хь

выходные массивы сохраняется одна точка Ь8 = (Ьь + ЬЬ)/2, х8 = (хь + хЬ)/2. Из обеих то-

х хь системы (??)■

Численное решение в скользящем режиме находится методом Рунге — Кутты для правой части (??)■ Причем при вычислении правой части системы (??) можно сократить количество выполняемых арифметических операций, учитывая условие Fij (х; хь) = Fji(xь; х). Также отметим, что значение производной функции д-(х) можно заменить выражением

^дг- х) ^ ^дг- (хь) ^ дг-(х ) - дг-(х)

гу> Н Гр /V»_ /V»

Завершение скользящего режима так или иначе связано с пересечением одной из траекторий или обеими траекториями одновременно некоторой поверхности разрыва. Если одна из траекторий перескочила на другую сторону поверхности разрыва (или обе поменялись сторонами), то либо траектория должна покинуть поверхность в силу завершения скользящего режима, либо шаг интегрирования необходимо уменьшить и продолжить скольжение. Если же одна из траекторий или обе пересекли другую поверхность разрыва, то это означает, что траектория достигла точки пересечения поверхностей разрыва и необходимо детально изучить дальнейшее направление движения траектории.

Если оценка погрешности найденной точки решения на данном шаге меньше заданной точности, то значение точки сохраняется в массивы X и Т. Если же шаг оказался

неудачным, то итерация алгоритма повторяется снова, стартуя с предыдущей сохраненной точки. Оценка погрешности численного решения производится с помощью метода Ричардсона. Дополнительно вычисляется решение с двумя шагами половинной длины. Эта схема используется как для непрерывного, так и для скользящего режима. Вложенные схемы оценки погрешности показали себя хуже метода Ричардсона [?]. поэтому в настоящей работе не использовались. В режиме пересечения погрешность оценивается условием завершения итерационной процедуры метода Ньютона.

При достижении конечной точки интервала времени Ъ = Ь алгоритм завершает свою работу.

6. Вычислительный эксперимент

Простой пример скользящего режима можно продемонстрировать на решении следующей системы. Состояние системы описывается вектором (х1}х2) € И2. Система имеет две области 01 и 02 непрерывного решения, разделенные кривой разрыва д(х1,х2) := х2 = 0, 01 = {(х\,х2) : х2 > 0}, 02 = {(х\,х2) : х2 < 0^. В области 01 решение системы задается системой дифференциальных уравнений

£ = -°.цх2 -2),

= 0.3(х1 — 1), (х\,х2) Е С\,

СЪ

(13)

В области С2

^ = —°Л(Х2 — 1),

= 0.3(х1 + 1), (х1, х2) Е 02.

Съ

(14)

Фазовый портрет системы (рис. ??, а) склеен из двух областей 01 и 02, решение в каждой из которых представляет собой цикл вокруг особой точки (тип особой точки центр), причем область 01 содержит особую точку (1; 2) для функции а особая точка (—1; 1) для функц ии Е2 в об л ас ть 02 не попадает. Траектории численного решения движутся против часовой стрелки. При склеивании областей образуется участок грани-

(—1; 0) (1; 0)

В остальной части границы пересечение происходит на тех ее участках, где траектории

б

Рис. 1. Фазовый портрет системы

(а); области фазовой плоскости (б)

непрерывны, но не являются гладкими, причем при х1 < —1 траектории пересекают границу х2 = 0 сверху вниз, а при х1 > 1 — снизу вверх.

Условно все фазовое пространство системы можно разделить на три непересекающиеся области А, В, С (рис. ??,б), траектории в которых имеют разное поведение. Область А прсдст^влябт (1; 2)

х2 = 0. Внутри этой области траектории являются циклическими и образуют концен-

ВА

О1 О2

входят в скользящий режим сверху и, сместившись до точки (1; 0), входят в цикл, пред-ст&вляютции АС

но в этом случае траектории входят в скользящий режим снизу и далее переходят на

А

Для вычислительного эксперимента, исследующего поведение предложенного численного алгоритма, используем систему (??), (??)• Нас будет интересовать поведение алгоритма в трех режимах: непрерывный, пересечение и скольжение. Все эти режимы реализуются в предложенной системе ОДУ.

А, В, С

дим начальные точки: (2; 2) из А, (0; — 4) из В, (0; 7) из С. Первая траектория находится в области непрерывности системы и требуется для изучения поведения алгоритма в условиях непрерывной задачи. Вторая (третья) траектория сначала проходит точку пересечения границы, затем входит в скользящий режим сверху (снизу), выходит со скольжения в точке (1; 0) и переходит в циклическое решение на границе области А. Эти траектории позволяют изучить поведение алгоритма в трех ситуациях: пересечение границы разрыва, переход траектории в скользящий режим и выход из скольжения.

О1 О2

ски. Решение в общем виде выглядит как

С1 С2

у/3

С1 = х1(Ь0) — Хъ С2 =--(х2(Ь0) — х2)-, гДе Х1, Х2 — координаты особой точки в

3

О1, О2 найденно го общего решения можно сконструировать решения для

АВС

Для области А с заданной начальной точкой х1(0) = 2, х2(0) = 2 получаем решение

х1(Ь) = С1 ГОБ

+ Х1,

Это решение является циклическим с периодом Ьа = —-—

Для области B с заданной начальной точкой x^0) = 0 x2(0) = —4 получаем решение, склеенное из четырех траекторий (I, II, III, IV):

г, ч i V3\ 5V3 . (V3 .

Xi(t) = со^ —1\ + — sin I — t ) — 1,

x2(t) = v^sln | I3tj — 5cos^j30 + 1.

Эта траектория достигает границы в точке x[(t1) = 2, x2(t1) = ü в момент времени

t1 = —-. Далее решение продолжается в области G1 траекторией 9

X1! (t) = cos | ^(t — h)^ + ^ sln^ß(t — t1) I + 1,

x2(t) = V^sln | l3(t — t1— 2cos(^3(t — t^ ) +2,

достигающей границы в точке х{г(t2) = 0, x1^(t2) = ü в момент времени t2 = t1 +

iüV3 ( + ( 1 \\ _ „

—3— In + arccos I — 7 I I. Затем решение продолжается траекторией на границе в скользящем режиме и описывается функциями

х[п (t) = 0.15(t — t2), xi11 (t) = 0.

Далее траектория достигает точки выхода из скользящего режима x{n(t3) = 1,

x2H(t3) = 0 в момент времени t3 = t2 +--. Затем решение продолжается циклически

15

xf (^"fsin^tt —13))+1,

xf (t) = —2^1^ — h)\+2} и траектория вновь достигает точку x[V(t4) = 1, x2V(t4) = 0 в момент времени t4 =

20V3n t> + —.

B

можно получить решение в области C.

Вычислительный эксперимент проводился для оценки эффективности работы предложенного алгоритма. Фазовая плоскость системы была разделена на три области (см. рис. ??, б), в каждой из которых траектории имеют качественно одинаковое поведение. Поэтому было выбрано по одной тестовой траектории в каждой области, найдены для них аналитические решения (см. выше) и вычислены несколько численных решений с разной заданной точностью для каждой начальной точки. Отслеживание траектории

проведено на интервале времени от заданной начальной точки до завершения первого цикла на границе области А.

Траектория решения в А (%) находится в области непрерывности

векторного поля направления правой части системы. Эта траектория используется для проверки работы алгоритма в стандартной (непрерывной) ситуации с целью определения не происходит ли увеличения времени счета по сравнению с таковым для традиционных алгоритмов и не нарушается ли сходимость метода на непрерьтвно-диффе-ренцируемом решении. Как показывают результаты вычислений (рис. 11.6). алгоритм сохраняет сходимость к точному решению. На рисунке приведены графики изменения порядка абсолютной погрешности ^ Аегг численного решения при разной заданной локальной точности То1 = 10-2,10-4,10-6,10-8. Видно, что погрешность получаемых решений не превышает заданной локальной точности.

Траектория решения в области В (рис. 11. а) стартует из точки (0; —4), пересекает границу областей непрерывности в точке (2;0), где происходит смена направления

(0; 0)

(1; 0)

А

грешности численного решения (рис. 11. б) показывают сходимость к точному решению при уменьшении значения заданной точности То1. Однако в отличие от непрерывного режима (см. рис. 11. а) при переходе через границу при входе в скользящий режим и выходе из него возникает скачок погрешности, что, как следствие, ведет к увеличению количества вычислительных операций для нахождения решения с заданной точностью (см. ниже рис. 11. а). Тем не менее и в этом случае удается найти решение с необходимой ТОЧНОСТЬЮ.

Траектория решения в области С (рис. 11, а) ведет себя аналогично траектории из В

тория приходит снизу (из области С2), а не сверху (из об ласти как в предыдущем

3.00

2.00

1.00

0.00

lg A err -2

-4 -6 -8 -10 -12

0.0

__-ог^

р-

---Т^4 -V-

г

10.0

20.0

б

30.0

Рис. 2. Траектория решения в области A с указанием точек численного решения, найденных с точностью Tol = 10" (а), и порядок абсолютной погрешности lg A err численного решения на траектории одного цикла при разных заданных значениях точности: Tol = 10"2 (!), 10"4 (2), 10-6 (3), 10"8 (4)

t

*2 4.00

2.00

0.00

-2.00

> У

-4.00

-1.00 0.00

1.00

a

2.00

1§ Aerr -2

-10

-12

-L-n—

2-0-—.

-- 4 ^

0.0

20.0

40.0

б

60.0

Рис. 3. Траектория решения в области В с указанием точек численного решения, найденных с точностью То1 = 10-2 (а), и порядок абсолютной погрешности ^ Аегг численного решения на траектории одного цикла при разных заданных значениях точности: То1 = 10-2 (!), 10-4 (2), 10-6 (3), 10-8 (4)

t

Рис. 4. Траектория решения в области С с указанием точек численного решения, найденных с точностью То1 = 10-2 (а), и порядок абсолютной погрешности ^ Аегг численного решения на траектории одного цикла при разных заданных значениях точности: То1 = 10-2 (□), 10-4 (о), 10-6 (•), 10-8 (v)

случае. Тем не менее мы сочли нужным включить эту траекторию в вычислительный эксперимент. Результаты (рис. ??, б) показали сходимость численного решения к точному при уменьшении величины заданной точности Tol.

Поведение алгоритма можно проследить на графиках изменения количества шагов интегрирования вдоль траектории в областях B и C (рис. ??). Заметим, что рост количества шагов происходит вблизи поверхности разрыва: в первый раз при пересечении поверхности, во второй раз при входе в скользящий режим и в третий раз при выходе из него. Увеличение заданной точности требует больших затрат на вычисление, поэтому происходит рост количества шагов интегрирования.

а б

Рис. 5. Изменение количества шагов интегрирования вдоль траектории в областях В (а) и С (б) при разных заданных значениях точности: То1 = 10-2 (□), 10-4 (о), 10-6 (•), 10-8 (v)

Вычислительные затраты алгоритма для нахождения численного решения адекватно оценить довольно сложно. Во-первых, вычислительные затраты численного метода решения ОДУ традиционно оценивают по количеству вычислений правых частей уравнения. Эта величина хорошо коррелирует с реальными вычислительными затратами для одно- и многошаговых методов, но в представленном алгоритме присутствуют итерационные процедуры, не использующие вычисления правых частей, но требующие временных затрат на сходимость. Во-вторых, вычислительные затраты, которые можно определить экспериментально, зависят от конкретной реализации и от компьютера, используемого для расчетов. Поэтому конкретные временные характеристики, выраженные в секундах работы алгоритма, сильно меняются в зависимости от применяемых модели компьютера и программной реализации.

Логично предположить, что удачной оценкой вычислительных затрат алгоритма будет количество шагов метода Эйлера, которое можно выполнить за время работы на данном компьютере в той же программной реализации. Именно такой критерий оценки производительности оыл выоран. Результаты эксперимента приведены в таблице, где представлены следующие параметры: То1 — заданная точность, М — количество принятых шагов метода для нахождения решения на заданном интервале, Nf — количество выполненных вычислении правых частей, Тс время (в секундах), потраченное на вычисление в данной реализации, Ме — количество шагов метода Эйлера, которое можно выполнить за время Тс, Уь — коэффициент увеличения шагов интегрирования (по отношению к предыдущей строке в таблице), Р — оценка порядка метода. Отметим, что для области А результаты эксперимента приведены для 10 циклов траектории, поскольку на одном цикле время точно измерить сложно в этом случае требуемое время составляет меньше 0.1 с. Для областей В ж С приведено значение времени вычисления одной траектории.

Проведенные вычисления показали, что среднее количество вычислений правых ча-

АВ

С

А

В

С

Оценка вычислительных затрат и порядка точности алгоритма

То1 N N1 Тс N Ун Р

Область А

10-2 108 1287 0.261 1569 — —

10-3 174 2079 0.411 3777 1.615 4.801

10-4 282 3375 0.480 4793 1.623 4.752

10-5 459 5499 0.691 7899 1.629 4.717

10-6 733 8787 1.202 15421 1.598 4.913

10-7 1177 14127 1.703 22796 1.608 4.849

10-8 1879 22551 2.593 35896 1.596 4.923

10-9 2992 35895 4.096 58020 1.592 4.954

10-10 4757 57075 6.489 93245 1.590 4.956

Область В

10-2 58 883 0.451 4366 — —

10-3 77 1169 0.461 4513 1.324 8.207

10-4 105 1538 0.551 5838 1.316 8.393

10-5 149 2105 0.641 7163 1.369 7.337

10-6 217 3006 0.791 9371 1.428 6.463

10-7 315 4230 0.902 11005 1.407 6.741

10-8 472 6178 1.222 15715 1.461 6.079

10-9 719 9199 1.572 20876 1.489 5.784

10-10 1114 14023 2.103 28683 1.524 5.461

Область С

10-2 61 1039 0.410 3763 — —

10-3 77 1263 0.491 4955 1.216 11.794

10-4 100 1589 0.610 6707 1.258 10.028

10-5 138 2100 0.601 6574 1.322 8.258

10-6 181 2545 0.621 6869 1.212 11.981

10-7 265 3636 0.812 9680 1.429 6.454

10-8 392 5207 1.061 13345 1.432 6.412

10-9 593 7678 1.332 17334 1.457 5.929

10-10 922 11689 1.833 24709 1.522 5.479

и 5.579. Этот результат показывает, что добавленные процедуры поиска точек пересечения траектории с границей и точек выхода из скользящего режима не понизили точность метода в целом и увеличили временные затраты не более чем на 24%. Кроме того, условие оценки погрешности на границе областей выбрано слишком сильным. Как показывает эксперимент, порядок метода на низкой заданной точности слишком велик и опускается до 5.461 на высокой заданной точности. Это означает, что условие оценки погрешности требует корректировки в сторону ее увеличения, что также позволит уменьшить вычислительные затраты на низкой заданной точности.

Таким образом, в представленной работе описан алгоритм численного решения системы ОДУ с правой частью, которая терпит разрыв на гладкой поверхности. В ходе реализации алгоритма были решены задачи по нахождению точки входа траектории численного решения на поверхность разрыва, определению дальнейшего поведения решения, вычислению решения в скользящем режиме, нахождению точки выхода со скользящего режима. На примере численного решения одной системы показана сходимость предложенного численного метода.

Список литературы

Айзерман М.А., Пятницкий Е.С. Основы теории разрывных систем. I // Автоматика и телемеханика. 1974. № 7. С. 33-47.

Айзерман М.А., Пятницкий Е.С. Основы теории разрывных систем. II // Там же. 1974. № 8. С. 39-61.

Матросов И.В. О единственности справа решений невырожденных алгебро-дифференциальных уравнений с разрывами // Там же. 2007. № 1. С. 11-19.

Уткин В.И. Скользящие режимы в задачах оптимизации и управления. М.: Наука, 1981. 368 с.

Di Bernardo М., Johansson К.Н., Vasca F. Self-oscillations and sliding in relay feedback systems: Symmetry and bifurcations // Intern. J. Bifurcation and Chaos. 2001. Vol. 11, No. 4. P. 1121-1140.

Di Bernardo M., Kowalczyk P., Nordmark A. Bifurcations of dynamical systems with sliding: Derivation of normal-form mappings // Phys. D. 2002. Vol. 170. P. 175-205.

Gear C.W., Osterby O. Solving ordinary differential equations with discontinuities // ACM Trans. Math. Software. 1984. Vol. 10. P. 23-44.

Gouze J.-L., Sari T. A class of piecewise linear differential equations arising in biological models // Dynamical Systems. 2002. Vol. 17. P. 299-319.

Johansson K.H., Barabanov A.E., Astrom K.J. Limit cycles with chattering in relay feedback systems // IEEE Transactions on Automatic Control. 2002. Vol. 247. P. 1414-1423.

Kowalczyk P., di Bernardo M., Champneys A.R. et al. Two-parameter nonsmooth bifurcations of limit cycles: Classification and open problems // Intern. J. Bifurcat. Chaos. 2006. Vol. 16. P. 601-629.

Kuznetsov Yu.A., Rinaldi S., Gragnani A. One-parameter bifurcations in planar Fillipov systems // Ibid. 2003. Vol. 13. P. 2157-2188.

Leine R.I., van Campen D.H., van de Vrande B.L. Bifurcations in nonlinear discontinuous systems // Nonlinear Dynamics. 2000. Vol. 23. P. 105-164.

Malmborg J., Bernhardsson B. Control and simulation of hybrid systems // Commun. in Nonl. Sci. and Numerical Simulat. 1997. Vol. 30. P. 337-347.

Di Bernardo M., Budd C.J., Champneys A.R., Kowalczyk P. Piecewise-smooth Dynamical Systems. Theory and Applications. Applied Mathematical Sciences. Berlin: Springer-Verlag, 2008.

Zhusubaliyev Z., Mosekilde E. Bifurcations and Chaos in Piecewise-Smooth Dynamical Systems. Singapore: World Scientific, 2003.

Филиппов А.Ф. Дифференциальные уравнения с разрывной правой частью. М.: Наука, 1985. 224 с.

Новиков Е.А., Шорников Ю.В. Численное моделирование гибридных систем методом Рунге — Кутты второго порядка точности // Вычисл. технологии. 2008. Т. 13, № 2. С. 99-105.

Cash J.R., Karp А.Н. A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides // ACM Trans. Math. Software. 1990. Vol. 16, No. 3. P. 201— 222.

[19] Died L., Lopez L. Sliding Motion in Filippov Differential Systems: Theoretical Results and A Computational Approach, http://www.math.gatech.edu/dieci/preps/DL-Fili.pdf (01.03.2009).

[20] Piiroinen P.T., Kuznetsov Yu.A. An event-driven method to simulate Filippov systems with accurate computing of sliding motions // ACM Trans. Math. Software. 2008. Vol. 34, No. 13. P. 1-24.

[21] Shampine L.F., Thompson S. Event location for ordinary differential equations // Comput. Math, with Appl. 2000. Vol. 39. P. 43-54.

[22] Guglielmi N., Hairer E. Computing breaking points in implicit delay differential equations // Advances in Comput. Math. 2008. Vol. 29, No. 3. P. 229-247.

[23] Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М : Мир, 1990. 512 с.

[24] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М : Мир, 1999. 685 с.

[25] Gupta G.K., Sacks-Davis R., Tescher P.E. A review of recent developments in solving ODEs // ACM Comput. Surveys. 1985. Vol. 17, No. 1. P. 5-47.

[26] Park Т., Barton P.I. State event location in differential-algebraic models // ACM Trans, on Modeling and Comput. Simulat. 1996. Vol. 6, No. 2. P. 137-165.

[27] Berzins M., Brankin R.W., Gladwell I. The stiff integrators in the NAG library // ACM SIGNUM Newsletter. 1988. Vol. 23, No. 2. P. 16-23.

[28] Enright W.H., Pryce J.D. Two FORTRAN packages for assessing initial value methods // ACM Trans. Math. Software. 1987. Vol. 13, No. 1. P. 1-27.

[29] Kamel M.S., Enright W.H., Ma K.S. ODEXPERT: An expert system to select numerical solvers for initial value ODE systems // Ibid. 1993. Vol. 19, No. 1. P. 44-62.

[30] Shampine L.F., Reichelt M.W. The MATLAB ODE suite // SIAM J. Sci. Comput. 1997. Vol. 18. P. 1-22.

[31] Коробицын В.В., Фролова Ю.В. Алгоритм численного решения дифференциальных уравнений с разрывной правой частью // Математические структуры и моделирование. Омский гос. ун-т, 2005. Вып. 15. С. 46-54.

[32] Коробицын В.В., Маренич В.В., Фролова Ю.В. Исследование поведения явных методов Рунге — Кутты при решении систем обыкновенных дифференциальных уравнений с разрывной правой частью // Математические структуры и моделирование. Омский гос. ун-т, 2007. Вып. 17. С. 19-25.

[33] Коробицын В.В., Фролова Ю.В., Маренич В.Б. Алгоритм численного решения кусочно-сшитых систем // Вычисл. технологии. 2008. Т. 13, № 2. С. 70-81.

Поступила в редакцию 1 апреля 2009 г.