Численное моделирование гибридных систем методом Рунге-Кутты второго порядка точности Текст научной статьи по специальности «Математика»

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

Том 13, № 2, 2008

Численное моделирование гибридных систем методом Рунге—Кутты второго порядка точности*

Е. А. Новиков

Институт вычислительного моделирования СО РАН, Красноярск, Россия

e-mail: novikov@icm.krasn.ru

Ю. В. Шорников Новосибирский государственный технический университет, Россия

e-mail: shornikov@inbox.ru

The control inequality for stability of the explicit two-phase Runge—Kutta method is obtained. For hybrid systems with the event function the control algorithm for modeling step is developed.

Введение

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

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

* Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 08-01-00621) и Президентской программы "Ведущие научные школы РФ" (грант № НШ-3431.2008.9).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

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

В настоящее время выделяют два подхода к контролю устойчивости [1, 2]. Первый связан с оценкой максимального собственного числа матрицы Якоби /у через ее норму с последующим контролем (наряду с точностью) неравенства \\Л/У|| < Д [1], где положительная постоянная Д зависит от размера области устойчивости. Очевидно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к значительному увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числа Лтах матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства |ЛЛтах| < Д [2]. Во всех рассмотренных ситуациях такая оценка фактически не вызывает увеличения вычислительных затрат [3]. Следует отметить, что в [2, 3] предполагается наличие в явном методе не менее трех стадий.

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

1. Метод второго порядка точности

При численном решении задачи Коши для системы обыкновенных дифференциальных уравнений

у' = /(у), у^о) = уо, ¿о < г < 4, (1)

рассмотрим явный двухстадийный метод типа Рунге—Кутты вида [4]:

Уп+1 = Уп + Р\к\ + ^2^2, кг = h/(уп), к2 = h/(у„ + вкг), (2)

где у и / — вещественные Ж-мерные вектор-функции; г — независимая переменная; Л — шаг интегрирования; кг и к2 — стадии метода; р1,р2, в — числовые коэффициенты, определяющие свойства точности и устойчивости (2). В случае неавтономной системы

у' = /(t,У), у(го) = Уо, ¿о < г < 4, схема (2) записывается в виде

уп+г = уп + Ргкг + Р2к2, кг = Л/(1п,уп), к2 = Л/(1п + вЛ,у,п + вкг).

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

Получим соотношения на коэффициенты метода (2) второго порядка точности. Для этого разложим стадии к1 и к2 в ряды Тейлора по степеням к до членов с к3 включительно и подставим в первую формулу (2). В результате получим

Уп+1 = Уп + (Р1 + Р2)/ + вР2 к2//п + 0.5в2Р2к3/П / + 0(к4),

где элементарные дифференциалы вычислены на приближенном решении уп, /'п =

д/(уп)/ду, = д2/(уп)/ду2. Точное решение у(1п+1) в окрестности точки Ьп имеет вид 3

у(<п+1) = у(^) + к/ + 0.5к2/'/ + к3 (/''/2 + /'2/) + 0(к4),

6

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

Р1 + Р2 = 1, вР2 = 0.5.

При данных соотношениях локальную ошибку 8п схемы (2) можно записать следующим образом:

¿п =2-—гк3/''/2 + 6 к3/'2/ + 0(к4).

Построим неравенство для контроля точности вычислений метода второго порядка. Для этого рассмотрим вспомогательную схему

уп+1,1 = уп + к1

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

£п,2 = уп+1 - уп+1,1 = Р2(к2 - Ь).

Для повышения надежности данной оценки выберем в = 1. Тогда стадия к1 вычисляется в точке ¿п, а к2 — в точке ¿п+1. Как показывают расчеты, использование производных решения в крайних точках шага приводит к более надежному неравенству для контроля точности вычислений. При в = 1 коэффициенты метода второго порядка определяются однозначно: р1 = р2 = 0.5, а локальная ошибка 8п и неравенство для контроля точности вычислений имеют, соответственно, вид

¿п = -1 к3// + 1 к3/'2/ + 0(к4), 0.5\\к2 - М < е, 1- 6

где \\ ■ \\ — некоторая норма в Ям, е - требуемая точность интегрирования.

Теперь построим неравенство для контроля устойчивости (2) предложенным в [2] способом. Для этого рассмотрим вспомогательную стадию к3 = к/(уп+1). Заметим, что к3 совпадает со стадией к1 , которая применяется на следующем шаге интегрирования, и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем стадии к1, к2 и к3 применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами. В результате получим

к1 = Хуп, к2 = (Х + X 2)уп, к3 = (Х + X2 + 0.5Х3 )уп,

Область устойчивости метода второго порядка точности

где X = ЛА. Легко видеть, что

к2 - кг = X2уп, 2(кз - к2) = X3у

Тогда согласно [3] оценку максимального собственного числа уп,2 = ЛЛп тах матрицы Якоби системы (1) можно вычислить по формуле

Уп,2 = 2шахг<г<^(|к3 - к2|/|к2 - М|). (3)

Область устойчивости метода второго порядка, полученная в инструментальной среде ИСМА [5], показана на рисунке.

Интервал устойчивости численной схемы (2) приблизительно равен двум. Поэтому для контроля ее устойчивости можно применять неравенство уп,2 < 2.

Полученную оценку (3) можно считать грубой потому что: 1) вовсе не обязательно максимальное собственное число сильно отделено от остальных; 2) в степенном методе применяется мало итераций; 3) дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг Лп+1 вычислим следующим образом. Новый шаг Лас по точности определим по формуле Лас = q1hn, где Лп есть последний успешный шаг интегрирования, а q1, с учетом соотношения к2-кг = 0(Ып), задается уравнением q2\\k2 - к1\ = е. Шаг Л3* по устойчивости зададим формулой Л3* = q2hn, где q2, с учетом равенства уп,2 = 0(Лп), определяется из уравнения q2Vh,2 = 2. В результате прогнозируемый шаг Л^н вычисляется по следующей формуле:

Лп+1 = шах[Лп, ш1п(Ласс, Л3*)]. (4)

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

2. Гибридная система

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений вида

у' = /(у), у(Ьо) = уо, д(у,Ь) < 0, (5)

где д(у,Ь) — событийная функция, или нелинейный предохранитель. Поскольку многие модели, представляющие интерес, линейны, будем рассматривать их как наиболее важный класс событийных функций. Заметим, что любой нелинейный предохранитель можно привести к линейному виду добавлением дополнительной фазовой переменной х = д(у,Ь). В результате задачу (5) можно переписать в виде

у = /(у),

х=^ (у)+д~д

х ду1 (У)+ дГ х < 0.

(6)

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

дп+1 = д(уп + 05кРп+1[/п + /п+в + К+1), (7)

где /п = /(уп) и /п+в = /(уп + вк/п). Отметим, что функция /,п+@ вычисляется в потенциально опасной точке. Поэтому рассмотрим событийную функцию

дп+1 = д(уп + К+1 /п, и + кРп+1), (8)

записанную на решении, полученном явным методом Эйлера. Разлагая дп+1 в ряд Тейлора и учитывая линейность д(у,Ь), имеем

. (ддп , . ддп\ ( ,

дп+1 = дп + кп+1 -ду/п + (9)

где

= ( . ) ддп = дд(уп^п) ддп = дд(у,п,ьп) дп =д(уп ду = ду ' дЬ = дЬ .

В итоге получили зависимость дп+1 от прогнозируемого шага крп+]_.

Теорема [6]. Выбор шага по формуле

К+1 = (7 - 1)^п/

д9п . + ддп

!п +

ду дЬ

(10)

где 7 € [0,1), обеспечивает поведение событийной динамики как устойчивой линейной системы, которая приближается к поверхности д(у,Ь) = 0. Кроме того, если д(у0,Ь0) < 0, то д(уп,Ьп) < 0 для всех п.

Доказательство. Подставив (10) в (9), имеем

дп+1 = тдп- (11)

Преобразовав рекуррентно (11), получим

дп+1 = тп+1до-

Учитывая, что 7 < 1, имеет место дп ^ 0 при п ^ ж. Кроме того, из соотношения 7 > 0 следует, что функция дп не меняет знака. Следовательно, при д0 < 0 будет выполняться дп < 0 для всех п. Тогда событийная функция никогда не пересечет потенциально опасную область д(уп,Ьп) = 0, что завершает доказательство. □

Теперь сформулируем алгоритм интегрирования с учетом прогноза шага через событийную функцию. Пусть решение уп в точке Ьп вычислено с шагом Лп. Кроме того, известны значения стадий к1 и к2 метода (1), а также вспомогательной стадии к3 с предыдущего шага. Тогда следующий шаг можно выполнить по определенному правилу. Шаг 1. Вычисляется fn = f (уп,Ьп).

Шаг 2. Вычисляются дп = д(уп,Ьп), ддп/ду = дд(уп,Ьп)/ду, ддп/дЬ = дд(уп,Ьп)/дЬ. Шаг 3. Вычисляется шаг Лрп+1 по формуле (10). Шаг 4. Вычисляется новый шаг Лп+1 по формуле

Лп+1 = шiп(hPn+l,hn+l-l),

где Лп+°1 задан соотношением (4).

Шаг 5. Выполняется следующий шаг интегрирования.

Заключение

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

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

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

[1] Shampine L.M. Implementation of Rosenbrok method // ACM Transaction on Mathematical Software. 1982. Vol. 8, № 5. P. 93-113.

[2] Новиков Е.А., Новиков В.А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Докл. АН СССР. 1984. Т. 277, № 5. С. 1058-1062.

[3] Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 197 с.

[4] Кнауб Л.В., ЛАЕвский Ю.М., Новиков Е.А. Алгоритм интегрирования переменного порядка и шага на основе явного двухстадийного метода Рунге—Кутты // СибЖВМ. 2007. Т. 10, № 2. C. 177-185.

[5] Шорников Ю.В., Новиков Е.А., СолодовниковА М.В. Программа исследования областей устойчивости численных методов «RStable» // Инновации в науке и образовании. 2007. № 3(26). С. 36.

[6] EsposiTO J., Kumar V., Pappas, G.J. Accurate event detection for simulating hybrid systems // Hybrid Systems: Computation and Control (HSCC). Vol. LNCS 2034. Springer-Verlag, 1998.

Поступила в редакцию 28 ноября 2007 г.