CC BY
16ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ГИБРИДНЫХ СИСТЕМ ЯВНЫМ МЕТОДОМ ТРЕТЬЕГО ПОРЯДКА В ИНСТРУМЕНТАЛЬНОЙ СРЕДЕ ИСМА
А. Е. Новиков, Е. А. Новиков*, Ю. В. Шорников**, Д. Н. Достовалов**
Сибирский федеральный университет, 660041, Красноярск, Россия *Институт вычислительного моделирования СО РАН, 660036, Красноярск, Россия **Новосибирский государственный технический университет, 630092, Новосибирск, Россия
УДК 519.622
На основе трехстадийной схемы типа схемы Рунге - Кутты третьего порядка точности разработан алгоритм переменного шага. Построены неравенства для контроля точности вычислений и устойчивости численной схемы. Для гибридных систем разработан алгоритм выбора шага интегрирования с учетом событийной функции.
Ключевые слова: метод Рунге - Кутты, контроль точности и устойчивости, переменный шаг, гибридные задачи.
The algorithm on the basis of the three-phasic scheme Runge-Kutta of the third order of accuracy and variable step is developed. Inequalities for the control of accuracy and stability of the numerical scheme are constructed. The algorithm of choice a step of integration is developed for hybrid systems with the guard condition.
Key words: method of Runge-Kutta, the accuracy and stability control, variable step, hybrid problems.
Введение. Для большого количества технических задач характерно наличие событий, которые означают наличие точек разрыва в первых производных от фазовых переменных. Такие комбинированные дискретно-непрерывные системы называются гибридными, или системами с переключением [1-3]. В них события происходят в моменты времени, соответствующие нулям некоторой алгебраической функции. Во многих случаях проблема усугубляется большой размерностью и жесткостью рассматриваемой задачи. При решении жестких задач большой размерности основные вычислительные затраты приходятся на декомпозицию матрицы Якоби [4-6], поэтому эффективность алгоритма интегрирования обычно повышается за счет замораживания матрицы Якоби [7]. В случае умеренно жестких задач можно применять алгоритмы интегрирования на основе явных методов с контролем точности вычислений и устойчивости численной схемы [8].
В данной работе создан алгоритм интегрирования переменного шага на основе трехстадийной схемы типа схемы Рунге - Кутты третьего порядка точности. Построены неравенства для контроля точности вычислений и устойчивости численной схемы. Алгоритм интегрирования может быть использован при решении гибридных задач, в том числе умеренно жестких. Разработан также алгоритм обнаружения событий гибридной системы. Рассмотренные
Работа выполнена в рамках АВЦП “Развитие научного потенциала высшей школы (2009-2010 гг.)” (код проекта РНП 2.1.2/4751), ФЦП “Научные и научно-педагогические кадры инновационной России” на 2009-2011 гг. (государственный контракт № П-297), а также при финансовой поддержке РФФИ (код проекта 08-01-00621) и Совета по грантам Президента РФ для государственной поддержки ведущих научных школ РФ (грант № НШ-3431.2008.9).
методы реализованы в среде моделирования гибридных систем ИСМА (инструментальные средства машинного анализа) [9].
1. Трехстадийный метод типа метода Рунге — Кутты. При численном решении задачи Коши для системы обыкновенных дифференциальных уравнений
у' = /(У), У(ьо) = Уо, ¿о < Ь < гк (1)
рассмотрим явный трехстадийный метод типа метода Рунге - Кутты вида
Уп+1 = Уп + рк + Р2^2 + рзкз, кг = к/(уп), к2 = к/(у,п + ^21^1), (2)
кз = к/(уп + взг к 1 + вз2 к2),
где у и / — вещественные Ж-мерные вектор-функции; Ь — независимая переменная; к — шаг интегрирования; к1, к2 и к3 — стадии метода; р1, р2, р3, в21, в31 и в32 — числовые
коэффициенты, определяющие точность и устойчивость (2). В случае неавтономной системы
у = /(t,У), У(ьо) = Уо, ¿о < Ь < ¿к, схема (2) записывается в виде
Уп+1 = Уп + рк + Р2к2 + Рзкз,
к1 = к/(Ьп,Уп), к2 = к/(Ьп + в21к,Уп + в2Л), кз = к/(Ьп + [вз1 + вз2] к Уп + вз1 к1 + Рз2к2)-
Ниже для упрощения выкладок будем рассматривать задачу (1). Однако разработанные методы можно применять для решения неавтономных задач.
Получим соотношения, определяющие коэффициенты метода (2) третьего порядка точности. Для этого разложим стадии к1, к2 и кз в ряды Тейлора по степеням к до членов с к4 включительно и подставим в первую формулу (2). В результате получим
Уп+1 = Уп + (р1 + р2 + рз)к/п + [в21р2 + (вз1 + вз2)рз]к2/п /п+
+ к [в21 вз2рз/п /п + 0,5(в21 р2 + (вз1 + вз2)2рз)/п/п] +
+ к4[0,5в21вз2рз/п/п /п + в21(вз1 + вз2)вз2рз/!^/п/п +
+ ([в21р2 + (вз1 + вз2 ) з рз] /6 } /"'] + 0(hЪ),
где элементарные дифференциалы вычислены на приближенном решении У .
Точное решение у (Ьп+1) в окрестности точки Ьп можно представить в виде
у(Ьп+1) = у(Ьп) + к/ + 0, 5к2/'/ + (к3/6)[/,2/ + /"/2] +
+ (к4/24)[/'з/ + 3/''/'/2 + /'/''/2 + /'''/з] + 0(к5),
где элементарные дифференциалы вычислены на точном решении у(Ьп). Сравнивая ряды для приближенного уп+1 и точного у(Ьп+1) решений до членов с кз включительно при условии Уп = у(Ьп), запишем условия третьего порядка точности схемы (2):
1) р1 + р2 + рз = 1;
2) в21Р2 + (^31 + Рз2)Рз — 1/2;
3) в21р2 + (вз1 + вз2)2рз — 1/3;
4) в21вз2Рз — 1/6.
В предположении уп — у(Ьп) локальную ошибку 8п схемы (2) можно вычислить по формуле 8п — у(1п+-\) — уп+1. Учитывая представления у^п+1) и уп+1 в виде рядов Тейлора, получаем
5п — К
24//з/ + ^24 — 2в21^з2р^ /'/"/2+ ^ — ^21^з2(^з1 + ,$з2)рз^ /''/У2 +
+ 0(К5).
В нелинейной системе условий третьего порядка точности численной схемы имеется два свободных коэффициента. Положим /321 — 0,5 и /Зз1 + вз2 — 1. Тогда на каждом шаге к1, к2 и кз вычисляются в точках Ьп, Ьп + 0,5К и Ьп + К соответственно. Условия третьего порядка запишем в следующем виде:
1) р1 + р2 + рз — 1;
2) 0,5р2 + рз — 1/2;
3) 0,25р2 + рз — 1/3;
4) Рз2рз — 1/3.
Из второго и третьего уравнений данной системы находим р2 — 2/3 и рз — 1/6. Из первого и последнего уравнений получаем р1 — 1/6 и /Зз2 — 2. Из равенства /Зз1 + /Зз2 — 1 следует вз1 — —1. В результате коэффициенты метода (2) имеют вид
1 2 1
р1 — рз — 6 > р2 — 3 > в21 — 2 > вз1 — —1 Д}2 — 2. (3)
При данных соотношениях локальную ошибку 8п схемы (2) можно записать следующим
образом:
¡>п — 24 К1 (/а/ — /"/'/2 — 1 /"'/з) — 0(Н5).
Построим неравенство для контроля точности вычислений метода третьего порядка. Для этого рассмотрим вспомогательную схему
Уп+1, 1 — Уп + гк + Г2к2,
где к1 и к2 определены в (2). Потребуем, чтобы данный метод имел второй порядок точности. Сравнение рядов для у(1п+-\) и уп+1; 1 показывает, что требование второго порядка будет выполнено, если г1 + г2 — 1 и в21г2 — 0,5. Отсюда и из (3) получаем г1 — 0 и г2 — 1. С использованием идеи вложенных методов ошибку еп, з метода третьего порядка можно оценить по формуле
£п, з — Уп+1 — Уп+1, 1 — (р1 — г1)к1 + (р2 — г2)к2 + рзкз.
Тогда неравенство для контроля точности вычислений имеет вид
\\к\ — 2^2 + к3\\< 6е,
где \ \ • \ \ — некоторая норма в Лм; е — требуемая точность расчетов.
Построим неравенство для контроля устойчивости (2) предложенным в [8] способом. Запишем стадии кг, к2 и к3 применительно к задаче у' = Ау, где А — матрица с постоянными коэффициентами. Используя соотношения (3), в результате получим
кг = Хуп, к2 = (Х + 0,5Х 2)уп, к3 = [X + X2 + X 3]у„,
где X = НА. Найдем коэффициенты ¿г, ¿2 и ¿3 из условия выполнения равенства
¿1 кг + ¿2к2 + йзкз = X ъуп.
Данное требование будет выполнено, если = ¿3 = 1 и ¿2 = -2. Нетрудно заметить также,
что к2 — к\ = X2уп. Тогда согласно [8] оценку максимального собственного числа уп, 3 = НХтах матрицы Якоби системы (1) можно вычислить по формуле
Vп, 3 = 0,5 тах (\к1 — 2к2 + к}3\/\кг2 — к\\). (4)
’ 1<г<М
Функция устойчивости Q(x) метода (2) с коэффициентами (3) имеет вид Q(x) = 1 + х + х2/2 + х3/6. Область устойчивости схемы (2), построенная в среде ИСМА [9, 10], приведена на рис. 1. Построены линии уровня ^(х)\ = д при д = 0,05; 0,30; 0,70; 1,00. Интервал устойчивости численной формулы (2) приблизительно равен 2,5, поэтому для контроля ее устойчивости можно использовать неравенство уп> 3 < 2,5.
Полученная оценка (4) является грубой, поскольку максимальное собственное число необязательно значительно отделено от остальных, в степенном методе применяется небольшое количество итераций, дополнительные искажения вносит нелинейность задачи (1). Поэтому в данном случае контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг Нп+г вычисляется следующим образом.
Новый шаг по точности Нас определим по формуле Нас = дгНп, где Нп — последний успешный шаг интегрирования; задается уравнением д3\\еп, 3\\ = е с учетом соотношения еп, 3 = О(Нп). Шаг по устойчивости Н8* зададим формулой Н81 = д2Нп, где д2 определяется из уравнения д2уп, 3 = 2,5 с учетом равенства уп, 3 = 0(Нп). В результате прогнозируемый шаг Нп+1 вычисляется по формуле
Нп+1 = тах[Нп, тт(Нас, Н8*)].
Данная формула позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Именно наличие данного участка существенно ограничивает применение явных методов при решении жестких задач.
2. Гибридная система. Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений вида
1т (х)
Рис. 1. Область устойчивости трехстадийного метода (2) третьего порядка точности
у'— /(у), у(Ьо) — уо, д(у,ь) < 0, (5)
где д (у,Ь) — событийная функция или нелинейный предохранитель. Поскольку многие модели, представляющие интерес, линейны, будем рассматривать их как наиболее важный класс событийных функций. Заметим, что любой нелинейный предохранитель можно привести к линейному виду добавлением дополнительной фазовой переменной х — д (у, Ь). В результате задачу (5) можно записать в виде
у' — /(у), х — |/(у) + х < 0.
Здесь принимается, что событийная функция линейна. Также для упрощения рассуждений будем предполагать, что исходная задача скалярная. Однако все приведенные ниже выводы будут распространяться и на системы уравнений. Особое внимание следует уделить выбору метода интегрирования. Полностью неявный метод использовать нельзя, поскольку он требует вычисления / (у) в потенциально опасной области, где модель не определена. В то же время известно, что явные методы слабоустойчивы. Поэтому будем использовать построенный выше метод (2) с контролем точности вычислений и устойчивости численной схемы. В этом случае событийная динамика описывается соотношением
дп+1 — д(уп+1, Ьп + Кп+1),
где уп+1 вычисляется по формуле (2). Отметим, что стадии к2 и кз в формуле (2) вычисляются в потенциально опасных точках Ьп + 0,5К и Ьп + К. Поэтому рассмотрим событийную функцию вида
дп+1 д(уп + Кп+1 /п, Ьп + Кп+1),
т. е. на решении, полученном с использованием явного метода Эйлера. Разлагая дп+1 в ряд Тейлора и учитывая линейность функции д(у,Ь), имеем
_ I ЪР ( ддп £ , ддп\ ^
дп+1 — дп + Кп+1 ^ ду /п + дЬ ) , (6)
где
— ( . ) ддп — дд(уп, Ьп) ддп — дд(уп, Ьп)
дп — д(уа,гп), 9у ду , дь дь .
В результате получаем зависимость дп+1 от прогнозируемого шага КРп+1.
Теорема [5]. Выбор шага по формуле
кп+1—(т—1)дп! (дп /п+^) > (7)
где 7 Е [0, 1), обеспечивает поведение событийной динамики как устойчивой линейной
системы, которая приближается к поверхности д(у,Ь) — 0. Кроме того, если д(у0, Ь0) < 0, то
д(уп, Ьп) < 0 для всех п.
Доказательство. Подставляя (7) в (6), имеем
дп+1 = 1дп.
Преобразовав рекуррентно данное соотношение, получим дп+1 — 7п+1д0. Так как ^ < 1, то дп ^ 0 при п ^ ж. Кроме того, из соотношения 7 > 0 следует, что функция дп не меняет знак. Следовательно, при до < 0 будет выполняться неравенство дп < 0 для всех п. Тогда событийная функция никогда не пересечет потенциально опасную область д(уп, Ьп) — 0. Теорема доказана.
Сформулируем алгоритм интегрирования с учетом прогноза шага через событийную функцию. Пусть решение уп в точке Ьп вычислено с шагом Кп. Кроме того, известны значения стадий к1} к2 и кз метода (2). Тогда можно выполнить следующий шаг:
Ш!аг 1. Вычисляется /п — /(уп, Ьп).
Ш!аг 2. Вычисляются
дп — д(уп, Ьп), ддп/ду — дд(уп, Ьп)/ду, ддп/дЬ — дд(уп, Ьп)/дЬ.
Ш!аг 3. Вычисляется шаг КРп+1 по формуле (7).
Ш!аг 4. Вычисляется новый шаг Кп+1 по формуле
Кп+1 — тт(кп+1, К+Л,
где Кп+1 задается соотношением
Кп+1 — тах[Кп, т1п(Кас, К5*)].
Ш!аг 5. Выполняется следующий шаг интегрирования.
3. Результаты расчетов. С целью тестирования работы алгоритма рассмотрена типичная гибридная система двух осциллирующих масс на пружинах [11]. Система может находиться в одном из двух локальных состояний: “Раздельно” или “Вместе”. Поведение системы в каждом из состояний описывается системой дифференциально-алгебраических уравнений.
При ^ < |к1п1 — к2п2 — х1(к1 — к2)| имеем
х\ — у1; У1 — к1(п1 — х1)/т1, а1 — к1(п1 — х1)/т1,
х'2 — У2, у'2 — к2(п2 — х2)/т2, а2 — к2(п2 — х2)/т2, (8)
при х1 — х2 и у1 > у2
а1 — [к1п1 + к2п2 — х1(к1 + к2)]/(т1 + т2),
У1 — [к1п1 + к2п2 — х1(к1 + к2)]/(т1 + т2),
х1 — уь (9)
а2 — [к1п1 + к2п2 — х1(к1 + к2)]/(т1 + т2), у2 — [к1п1 + к2п2 — х1(к1 + к2)]/(т1 + т2), х'2 — У2, в — —в,
где т1 и т2 — массы грузов; к1 и к2 — жесткости пружин; п1 и п2 — нейтральные координаты грузов; х1 и х2 — координаты грузов; у1 и у2 — скорости грузов; а1 и а2 — ускорения грузов; в — общая жесткость пружин в состоянии “Вместе”.
Ниже приведен текст компьютерной модели гибридной системы двух осциллирующих масс на языке LISMA [12]. k1=1; k2=2; n1=1; n2=2; m1=1; m2=1; x1=0; x2=3;
separate [s<abs(k1*n1-k2*n2-x1*(k1-k2))] is
s~=10;
x1'=v1;
v1'=k1*(n1-x1)/m1;
a1~=k1*(n1-x1)/m1;
x2'=v2;
v2'=k2*(n2-x2)/m2;
a2~=k2*(n2-x2)/m2;
from;
together [(x1>=x2) and (v1>=v2)] is s=10;
v1=(m1*v1+m2*v2)/(m1+m2);
v2=v1;
v1'=(k1*n1+k2*n2-x1*(k1+k2))/(m1+m2);
a1~=(k1*n1+k2*n2-x1*(k1+k2))/(m1+m2);
x1'=v1;
v2'=(k1*n1+k2*n2-x2*(k1+k2))/(m1+m2);
a2~=(k1*n1+k2*n2-x2*(k1+k2))/(m1+m2);
x2'=v2;
s'=-s;
from separate;
Результаты анализа (8) и (9) в инструментальной среде ИСМА [9] с разработанным алгоритмом обнаружения (рис. 2,a) совпадают с результатами расчета эталонной модели в системе HyVisual [13]. Традиционный анализ системы без алгоритма обнаружения приводит к некачественным результатам (рис. 2,6).
н----------1-----1-------1----------1------1-------1-------1-------1-------1-------1— t —і-------1--------1-------1-------1------1--------1-------1-------1-------1--------н-
0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
Рис. 2. Результаты компьютерного анализа гибридной системы двух масс, полученные с использованием алгоритма обнаружения (а) и без него (б)
Заключение. Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, поскольку оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f. Такая оценка получается грубой, однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий.
Предложенный способ прогноза шага через событийную функцию распространяется на все явные методы типа метода Рунге - Кутты, поскольку при выборе шага используется только явный метод Эйлера.
Список литературы
1. Esposito J., Kümar V., Pappas G. J. Accurate event detection for simulating hybrid systems // Proc. of the 4th Intern. workshop “Hybrid systems: computation and control (HSCC 2001)”, Rome (Italy), March 2001. Berlin: Springer-Verlag, 2001. V. 2034 P. 204-217.
2. Новиков Е. А., Шорников Ю. В. Численное моделирование гибридных систем методом Рунге -Кутты второго порядка точности // Вычисл. технологии. 2008. Т. 13, № 2. С. 98-104.
3. Новиков A. E., Новиков E. A., Шорников Ю. В. Моделирование гибридных систем явными методами // Вестн. КрасГАУ. Ресурсосберегающие технологии механизации сел. хоз-ва. 2009. № 5 С. 88-92.
4. ХаЙрер Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциальноалгебраические задачи / Э. Хайрер, Г. Ваннер. М.: Мир, 1999. 685 c.
5. Новиков Е. А., Шитов Ю. А., Шокин Ю. И. Одношаговые безытерационные методы решения жестких систем // Докл. СССР. 1988. Т. 301, № 6. С. 1310-1314.
6. Новиков Е. А., Шитов Ю. А., Шокин Ю. И. О классе (m, &)-методов решения жестких систем // Журн. вычисл. математики и мат. физики. 1989. Т. 29, № 2. С. 194-201.
7. Новиков Е. А., Шитов Ю. А. Алгоритм интегрирования жестких систем на основе (m, &)-метода второго порядка точности с численным вычислением матрицы Якоби. Красноярск, 1988. (Препр. ВЦ СО АН СССР; № 20).
8. Новиков Е. А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 197 с.
9. Свидетельство об офиц. регистрации программы для ЭВМ № 2005610126. Инструментальные средства машинного анализа / Ю. В. Шорников, В. С. Дружинин, Н. А. Макаров, К. В. Омельченко, И. Н. То-милов. М.: Роспатент, 2005.
10. Шорников Ю. В., Новиков Е. А., ДостовАлов Д. Н. Анализ устойчивости явных методов Рунге - Кутты в инструментальной среде ИСМА // Пробл. информатики. 2009. № 3. С. 42-51.
11. Шорников Ю. В. Моделирование сложных динамических и гибридных систем в ИСМА // Науч. вестн. НГТУ. 2007. № 1. С. 79-88.
12. Свидетельство об офиц. регистрации программы для ЭВМ № 2007611024. Программа языкового процессора с языка LISMA (Language of ISMA) / Ю. В. Шорников, И. Н. Томилов. М.: Роспатент, 2007.
13. Lee E. A., Zheng H. Operational semantics of hybrid systems // Proc. of the 8th Intern. workshop “Hybrid systems: computation and control (HSCC 2005)”, Zurich (Switzerland), March 2005. Berlin: SpringerVerlag, 2005. V. 3414. P. 25-53.
Новиков Антон Евгеньевич — ассист. Института фундаментальной подготовки Сибирского федерального университета; e-mail: jibuti88@mail.ru;
Новиков Евгений Александрович — д-р физ.-мат. наук,
гл. науч. сотр. Института вычислительного моделирования СО РАН; e-mail: novikov@icm.krasn.ru;
Шорников Юрий Владимирович — д-р техн. наук, проф. Новосибирского государственного технического университета; тел. (383) 346-11-00, e-mail: shornikov@asu.sc.nstu.ru, shornikov@inbox.ru;
Достова,лов Дмитрий Николаевич — магистр техники и технологий,
e-mail: dostovalov.dmitr@mail.ru.
Дата поступления — 24.05.10
