CC BY
30
Вычислительные технологии
Том 19, № 6, 2014
Символьные вычисления в моделировании
>к
и качественном анализе динамических систем*
А. В. Банщиков, Л. А. БУРЛАКОВА, В. Д. Иртегов, Т.Н. Титоренко
Представлены написанные на языке программирования системы компьютерной алгебры Mathematica комплексы программ для моделирования и качественного анализа в символьном виде динамических систем. В качестве динамических систем, в частности, рассматриваются механические системы и электрические цепи. Описаны используемые алгоритмы и их программная реализация в виде пакетов расширения системы компьютерной алгебры Mathematica. Функциональные возможности комплексов продемонстрированы на конкретных примерах.
Ключевые слова: моделирование, качественный анализ, механическая система, электрическая цепь, система компьютерной алгебры, комплекс программ.
Введение
При построении математической модели (дифференциальных уравнений движения) для сложных объектов и её качественном анализе исследователи сталкиваются с необходимостью оперировать громоздкими аналитическими выражениями. Эти трудности бывают настолько велики, что существенно сужают класс исследуемых моделей.
Проблемы достоверности, точности вычислений и некоторые другие вопросы вычислительной математики, а также вопросы ускорения и наглядности исследовательского процесса могут быть частично сняты, если в качестве инструмента решения задач выбрана система компьютерной алгебры (СКА). Производительность СКА и эффективность её использования существенно возрастают в связи с увеличением быстродействия и оперативной памяти современных компьютеров, позволяющих обрабатывать большие объёмы символьной информации. Символьные вычисления на ЭВМ уже давно перешли в разряд рабочих вычислительных методов с разнообразными приложениями.
Как видно из публикаций, иногда используется подход к применению СКА в качестве калькулятора для решения отдельно взятой задачи. Существует и другой подход, когда на базе внутреннего языка выбранной СКА или какого-либо универсального языка программирования разрабатывается программное обеспечение для решения определённого класса задач.
В настоящей работе описываются комплексы программ [1-3] для моделирования и качественного анализа динамических систем, дифференциальные уравнения которых можно получить по некоторой характеристической функции. В качестве таких систем, в частности, рассматриваются механические системы и электрические цепи. Написанные на языке программирования СКА Mathematica комплексы программ предназначе-
1 Институт динамики систем и теории управления СО РАН, Иркутск, Россия
Контактный e-mail: bav@icc.ru
* Работа поддержана Программой фундаментальных исследований Президиума РАН № 17.1 и Советом по грантам Президента РФ для государственной поддержки ведущих научных школ (НШ-5007.2014.9).
ны для решения указанных задач полностью в символьном виде. Символьные вычисления оказались достаточно эффективными на этапах построения математических моделей и анализа их свойств. На решение таких задач и ориентированы представленные комплексы программ. Комплексы позволяют по геометрическому описанию механической системы или графу электрической цепи построить характеристическую функцию системы, которая используется как для вывода уравнений движения (состояния цепи), так и для анализа свойств решений полученных уравнений. По известным характеристическим функциям и уравнениям движения возможно исследование систем, отличных от механических и электрических, например, на алгебрах Ли.
Алгоритмическую основу комплексов образуют классические методы аналитической механики, теории устойчивости движения и их развитие. Сочетание алгоритмов предметной области с символьными алгоритмами систем компьютерной алгебры позволяет создавать достаточно мощные инструменты для исследования динамических систем. В статье приводится краткое описание алгоритмов, используемых в рассматриваемых программных разработках, что даёт возможность представить круг задач, решаемых комплексами, подходы к моделированию и анализу исследуемых систем. Программная реализация используемых алгоритмов описывается в виде пакетов расширения СКА ЫаЬНешаЫса. Для демонстрации функциональных возможностей обсуждаемых комплексов программ выбраны достаточно простые механические и электрические системы. Исследование сложных механических систем конкретизировано в заключении статьи, а ссылки на соответствующие работы представлены в списке литературы.
1. Математическая модель механической системы
В качестве механической рассматривается система твёрдых тел. Тела соединены так, что либо у них существуют общие точки, либо они могут совершать винтовые движения друг относительно друга. Предполагается, что система находится под воздействием сил потенциальной и непотенциальной природы. Для её описания используется формализм Лагранжа, согласно которому движение системы, находящейся под воздействием только потенциальных сил, характеризуется функцией Лагранжа Ь = Т + и, где и = и(д) — силовая функция, Т = Т(д, д, ¿) — кинетическая энергия, д = (д1,. . . , дп) — обобщённые координаты, д = (с[1,... , дп) — обобщённые скорости.
1.1. Описание механической системы
Алгоритм построения функции Лагранжа основывается на геометрическом (рис. 1) и кинематическом описании исследуемой системы тел. На рис. 1 тело В1 является носителем, его положение определяется относительно инерциальной системы осей Е0. С твёрдым телом В1 связана система осей Е1, имеющая начало в точке О1 Е В1, а с телом В к — Ек, имеющая начало в точке О к Е В к. Возможно также использование полуподвижной системы осей [4]. Поворот Ек относительно Е-7 определяется матрицей а-к, элементы которой являются функциями выбранных обобщённых координат.
Для ввода в компьютер в соответствии с конкретной конфигурацией механической системы указывается количество тел N в системе, а для каждого тела Вк (к = 1,... , N) задаются:
• ] — номер тела В-, с которым соединено тело Вк, или номер системы осей Е7, относительно которой определяется положение тела Вк (] < к для к > 1 и ] = 0 для В1);
Рис. 1. Геометрическое описание механической системы
• номера осей вращения и соответствующие им идентификаторы углов поворотов;
• r
• r
Ok k
Ck
радиус-вектор точки Ок соединения тел В3 и Вк в системы осей Т3; радиус-вектор точки Ск (центр масс тела Вк) в осях Т
k
• v
Ok
вектор относительной линейной скорости точки Ок в проекциях на оси Т3 (для тела В1 вектор абсолютной линейной скорости точки Ок);
• масса тела тк и его тензор инерции ©°к относительно точки Ок;
• список обобщённых координат.
Далее номера осей вращения и соответствующие им идентификаторы углов поворотов для каждого тела Вк используются для автоматического построения:
• матрицы направляющих косинусов а3к, определяющей угловое положение системы осей Тк относительно Т3;
• векторов относительной угловой скорости ш3к к-го тела в системе осей Т3 и абсолютной угловой скорости тела шк = + ш3к;
• вектора абсолютной линейной скорости "0к = "о- + [ шк х го ] + "о точки Ок.
1.2. Кинетическая энергия и силовая функция
Кинетическая энергия тела Вк и всей системы находится соответственно по формулам
N
Tk = 2 mk vOk + 2 Шк • SOk • Шк + mk [ vOk x Шк ] • rCk, T = ^ Tk. (1)
k=i
Силовая функция тела Вк в ньютоновском поле тяготения с точностью до второго порядка малости относительно величин, равных отношению характерных размеров системы к расстоянию от точки О1 до неподвижного притягивающего центра О (см. рис. 1), вычисляется следующим образом
Uk =
vmk vmk
r
[(rO, • rOk ) + (rO, • rCk)]
vmk
(r.
1 )2
Ok)
+ (r
Ok
rCk)
3 v
2Г3
(ek )T ©°k ek + 2V3 (0Oik + ©Ok + ©Ok) +
+
3 vmk
O,
■Ok)
+
O,
rM • (r°
Ok
O1 ACk,
+ const.
(2)
1
2
0
0
k
5
r
N
Соответственно силовая функция всей системы определяется как V = ^ Vк. Здесь
к= 1
V — постоянная тяготения; г = | г°х1 — расстояние от точки 01 до притягивающего центра О; в°к (г = 1, 2, 3) — главные моменты инерции тензора ®°к; вк = а1кв1 — матрица 3 х 1 направляющих косинусов углов, образуемых вектором г° с осями Ек; в1 = о>01Х (х — матрица 3 х 1 направляющих косинусов углов, образуемых г° с СО Е0).
Силовая функция всей системы и входящих в неё тел Вк в поле постоянной тяжести имеет вид
N
0
V = ^ Цк, Цк = Шк ^ ■ гСк) + сопв^ (3)
г
к= 1
где g — вектор ускорения постоянной силы тяжести, г^ = г° + г° + гС. 1.3. Уравнения движения
По известной кинетической энергии (1) и силовой функции (2) (или (3)) системы вычисляется функция Лагранжа
1 п п п
Ь = Т2 + Т1 + То + V = 2 ^^ Мз (д, г) д Т + ^ Аг(д, г) д + Ао(д, г), (4)
г=1 з=1 г=1
где Агз = Азг, А0(д,г) = Т0(д,г) + V(д). По функции Лагранжа могут быть получены дифференциальные уравнения движения в форме уравнений Лагранжа второго рода, Рауса, Гамильтона [4]. Например, уравнения движения в форме уравнений Лагранжа второго рода строятся исходя из формулы (4):
^ л ■■ , (дАгг , дАгз дАгД . . ^ (дАг дАу \ .
£Агг"г +1Ц3=11дз + —Ж) "гдз + £ ( Ч - % )д+
А дАгг . дАг дА0 д& , А г 1 ^ ^
+ £ ---"г + --Г — -ддГ = + А (г = 1,...,п). (5)
Г=1
Здесь — дII/ддг + Qг — непотенциальные силы, слагаемое — дИ/ддг соответствует дисси-пативным силам, определяемым функцией Релея II.
Для построения уравнений Гамильтона выполняется переход к гамильтоновым переменным дг, рг (рг = дЬ(д,д)/ддг), а при наличии у системы циклических координат для построения уравнений Рауса осуществляется переход к переменным Рауса.
2. Математическая модель электрической цепи
В качестве электрической рассматривается цепь, в которой произвольным образом соединены элементы сопротивления (Я), катушки индуктивности (Ь), конденсаторы (С), источники тока (I) или напряжения (Е). В случае линейных электрических цепей указанные элементы имеют постоянные (не зависящие от токов и напряжений в ветвях цепи) характеристики. Электромеханические аналогии [5] позволяют использовать методы аналитической механики для описания электрических цепей и тем самым применять методы исследования механических систем к электрическим цепям.
2.1. Функция Лагранжа линейной электрической цепи
Согласно [5], если в качестве независимых переменных выбраны контурные токи, то состояние линейной электрической цепи можно описать при помощи функций Лагранжа Ь и Релея II
-.га 1 ™ 1 п
Ь = т - П, т = 2 £ Ыэ, П = ^ £ С., 1=2 £ Ыэ. (6)
¿,.7=1 г,э = 1 ¿э ¿,.7=1
Здесь Т — кинетическая, П — потенциальная энергии цепи; д1,...,дп — контурные заряды, д!,... ,дп — их производные по Г; Д., , С. — соответственно омическое сопротивление, индуктивность, ёмкость, п — число независимых контуров цепи.
Если в качестве независимых переменных выбраны напряжения в узлах цепи, то кинетическая и потенциальная энергии и функция Релея т +1 узловой электрической цепи принимают вид
т т т . .
т=1 £ с., п = 2 £ ^, 1 = 1 £ ^, (7)
¿,.7 = 1 ¿,. = 1 ¿,. = 1
где м1,... , ит — узловые потенциалы, и 1,... , ит — их производные по Г.
Вычисление функций (6), (7) на компьютере методом контурных токов и узловых потенциалов опирается на представление электрической цепи в виде графа, где множеству его вершин соответствует множество узлов цепи, а множеству рёбер — множество ветвей цепи. Задача определения токов и напряжений в ветвях цепи сводится к нахождению множества фундаментальных циклов графа, построение которых выполняется стандартными алгоритмами на графах. Для каждой а-й ветви цепи по известному в ней току (напряжению) вычисляются кинетическая То и потенциальная По энергии и функция Релея I. Указанные величины всей цепи определяются как суммы величин То, По, I, взятые по всем ветвям цепи.
По известной функции Лагранжа линейной электрической цепи строятся уравнения состояния цепи в форме уравнений Лагранжа второго рода (5), Гамильтона, Рауса.
2.2. Характеристическая функция нелинейной электрической цепи
В качестве нелинейной рассматривается электрическая цепь, принадлежащая классу "полных" в смысле [6]. Дифференциальные уравнения такой цепи имеют вид
Ьк ^ )-£ = (к =1,...,г), С (иг)—1 = — — (1 = г + 1,...,г + т), (8)
ас дгк ас диг
где ¿к(г^) — индуктивности; С(иг) — ёмкости; ¿1,...,гг — токи, протекающие через катушки индуктивности; иг+1,... , иг+т — напряжения на конденсаторах.
Функция Р(г, и) = Р(г) — С(и) + (г, 7«т) описывает часть цепи, содержащую элементы сопротивления и источники напряжения (тока), и называется смешанным потенциалом. Здесь г = (¿1,..., гг), и = (иг+1,..., иг+т) — вектор-строки, 7 — матрица г х т, элементы которой равны 0, ±1.
В работе [6] предложен метод построения функции Р(г, и) по графу цепи. В соответствии с этим методом электрическая цепь разбивается на "элементарные" ветви,
содержащие только один элемент (сопротивление R, катушку индуктивности L, конденсатор C, источник напряжения E или тока I). Каждому ребру графа цепи ставится в соответствие одна "элементарная" ветвь. Вершинам графа соответствуют узлы цепи, которые соединяют "элементарные" ветви.
Посредством стандартных алгоритмов на графах множество ветвей цепи разделяется на множество ветвей Qj, токи в которых определяются переменными i\,... ,ir, и множество ветвей Qu, напряжения в которых определяются переменными ur+1,... ,ur+m. Показано, что такое разделение всегда возможно для "полных" цепей. Соответственно для ветвей Qj и Qu вычисляются функции
F(i) = ^ У up dip, G(u) = ^ У ip dup.
p>r+m, Q г p>r+m, Qi p
Здесь Г — кривая в пространстве токов и напряжений цепи, вдоль которой удовлетворяются законы Кирхгофа; ip, up — ток и напряжение в ^-й ветви цепи.
Для нахождения слагаемого (i,^uT) граф цепи разбивается на r независимых контуров Лр (р = 1,... ,r) c контурными токами i1,... ,ir, тогда
r
(i,Yu) = ip S ±u°.
р=1 ЛрПто
Здесь ip — ток в р-м контуре, ua — напряжение в a-й ветви, тс — дерево графа цепи, образованное ветвями с элементами C.
Функция P(i, u), кроме задачи составления дифференциальных уравнений цепи, используется также для качественного исследования данных уравнений как аналог функции Ляпунова (см. [6]). Для этого вычисляется её полная производная по времени в силу дифференциальных уравнений цепи и проверяются свойства P(i,u) и производной, например их знакоопределённость, в соответствии с теоремами Ляпунова (Четаева и др.).
3. Комплексы программ LinModel и Circuits
На основе приведённых выше алгоритмов на языке программирования СКА Mathematica написаны программы для построения характеристических функций и уравнений движения механических систем и электрических цепей. Программы объединены в пакеты (комплексы программ) LinModel и Circuits. В системе Mathematica пакеты (называемые пакетами расширения) используются для объединения программ, предназначенных для решения задач определённой предметной области, и позволяют расширить функциональные возможности СКА. После загрузки такого пакета в среде Mathematica с его программами можно работать как со встроенными функциями системы.
Пакет LinModel содержит программы построения функции Лагранжа механической системы по её геометрическому описанию, уравнений движения в форме уравнений Лагранжа второго рода, Рауса, Гамильтона. Для вывода уравнений Гамильтона используется как классический метод, так и схема Дирака, применяемая к системам с "вырожденной" кинетической энергией.
Пакет Circuits включает программы построения функции Лагранжа линейной электрической цепи (методами контурных токов и узловых потенциалов) и смешанного потенциала нелинейной цепи по заданному графу, а также программы вывода дифференциальных уравнений состояния цепи по заданному или вычисленному потенциалу.
Специфика применения пакетов демонстрируется на двух примерах. Пример 1. Задача о движении сферического маятника (рис. 2) в поле сил постоянной тяжести. Введём систему осей Охсух (инерциальную), Ож1у1г1 (связанную с телом) и составим следующее описание механической системы для пакета ЫпМоёе1.
Число тел — 1. Тело 1: т — масса тела; г^, = (0, 0, 0) — радиус-вектор точки О; = (0, 0, 0) — вектор линейной скорости точки О; гС = (0,0,1) — радиус-вектор цен/ т12 0 0 \
тра масс тела; ©° = I 0 т 12 0 I — тензор инерции тела относительно точки О;
000
(3,1) — номера осей вращения; (<^,0) — углы вращения, т.е. положение тела определяется последовательностью двух поворотов: вокруг оси 3 (Ог) на угол ^ и вокруг
Рис. 2. Геометрическое описание маятника
Рис. 3. Интерфейс пакета ЫпМоёе!
оси 1 (Ox 1) на угол в. В использованных обозначениях верхний индекс 0 означает, что движение тела описывается относительно системы осей Oxyz, 1 — относительно Ox1y1z1.
На рис. 3 показан фрагмент работы с пакетом LinModel: загрузка пакета в среде СКА Mathematica, вызов функции пакета для построения лагранжиана механической системы, ввод исходных данных для представленного примера. Ввод в программу описания системы выполняется либо в режиме диалога посредством диалоговых панелей, либо через файл, формируемый пользователем. В первом случае вначале вводятся данные, относящиеся ко всей системе в целом, затем — данные для каждого тела. Введённые данные сохраняются в архиве данных (txt-файлы), их можно отредактировать и использовать повторно для ввода через файл.
После ввода описания механической системы вычисляется функция Лагранжа, которая сохраняется в архиве данных и используется другими программами пакета в качестве исходных данных, например, для вывода уравнений движения.
Ниже приведены построенные пакетом LinModel функции Лагранжа и Рауса:
p2coseche
2L = mi2(sin2 вф2 + в)2) — 2mgl cos в, 2R = ш12в2-----2mgl cos в, где p = const;
ml2
уравнения Лагранжа
1в — l cos в sin вф2 — g sin в = 0, ф + 2^вф в = 0, (9)
уравнение Рауса
функция Гамильтона
mi2 в — p2cosec2f — mgl sin в = 0, (10)
ml2
(pfcosec2 в + p2)
2H = ^---^ + 2mgl cos в,
ml2
где p1, p2 — импульсы.
Пример 2. Рассмотрим использование пакета Circuits для построения смешанного потенциала и уравнений состояния нелинейной электрической цепи, эквивалентная схема и граф которой приведены на рис. 4 и 5. Здесь i1, i2 — токи, протекающие через катушки индуктивности; ua, u4, u5, u6 — напряжения на конденсаторах; C1, C2, Ca, C4, L1, L2, I0, I01, I02 — параметры цепи, соответствующие постоянным ёмкостям, индуктивностям, токам. Нелинейными элементами цепи являются сопротивления, вольт-амперные характеристики которых описываются функциями f1(u6), f2(u4), f3(u3), f4(u5). Графики последних предполагаются качественно совпадающими: они находятся в первом и третьем квадрантах координатной плоскости (i, u), проходят через начало координат и имеют по две экстремальные точки в первом квадранте.
Описание графа цепи (см. рис. 5) для ввода в программу выглядит следующим образом:
{{{1, 2},L1}, {{2, 3},М, {{2, 3},C1}, {{3, 4},L2}, {{4, 5},Io2>, {{4, 5},C4}, {{4, 5},f4>, {{3, 5},C2}, {{3, 5},f2}, {{2, 5},Io}, {{1, 5},Ca}, {{1, 5}, fa}, {{1, 5},«}.
Здесь fj (j = 1,... , 4) — обозначения для функций fj (uk).
Рис. 4. Электрическая цепь
Рис. 5. Граф цепи
Как и в предыдущем случае, данные в программу вводятся посредством диалоговой панели либо файла, формируемого пользователем. Результат работы пакета Circuits имеет вид
Пи = {{{2, 3},Ci}, {{3, 5},C2}, {{1, 5},C3}, {{4, 5},C4}, {{2, 3},/i>, {{3, 5},/2}, {{1, 5},/з}, {{4, 5},/4}, {{2, 5},/о}, {{4, 5},/02, {{1, 5},/qi}; = {{{1, 2},Li}, {{3,4},L2}}; Л2 = {{{3, 4}, L2}, {{4, 5}, C4}, {{5, 3},C2}}; Л1 = {{{1, 2},Li}, {{2, 3},Ci}, {{3, 5}, C2}, {{5,1},Сз}};
из И4 И5 И6
G(u) = /oiU3 + /02И5 - /o(u4 + Ue) + J /3(u) du + J /2(u) du + J /^(u) du + J /i(u) du;
0000
F(i) = 0; (i,YuT) = ii(u4 + u6 - u3) + ¿2(u - u4); P(i,u) = -G(u) + (i,YuT).
Отметим, что в программе построения смешанного потенциала цепи ток в ветви {na, ng} в направлении от узла na к узлу ng считается положительным, если na < ng, и отрицательным в противном случае. Аналогично определяется знак напряжения в ветви цепи.
Функция P(i,u) служит исходными данными для программы пакета Circuits, формирующей уравнения состояния цепи вида (8). Ниже приведены построенные уравнения:
di1 di2 L1— = u4 + u6 - ua, L2— = u5 - u4, dt dt
dua . . . du4 . . . .
= I01 + i1 + fa(ua), C2— = i2 - i1 + /2^4) - Io, dt dt
C4^TT = I02 - i2 + /4^5), C1= -i1 - Io + A(ue). (11)
dt dt
Пример построения функции Лагранжа линейной электрической цепи можно найти в работе [3].
4. Качественный анализ уравнений движения
В этом разделе описываются алгоритмы выделения стационарных решений и инвариантных многообразий (ИМ), полученных в символьном виде уравнений движения, и исследования их устойчивости на основе метода Рауса — Ляпунова [7] и его модификаций [8, 9] и теорем Ляпунова об устойчивости по первому приближению.
4.1. Выделение стационарных решений методом Рауса — Ляпунова
Пусть для автономной системы обыкновенных дифференциальных уравнений
Х = X(x), x е Rn, (12)
известно некоторое число алгебраических (или сводящихся к ним) первых интегралов
Vj(x) = Cj (j = 0,1,... ,k), где Cj — постоянные интегралов. В соответствии с методом
Рауса — Ляпунова стационарные решения системы вида (12) могут быть получены как
решения задачи на условный экстремум первых интегралов данных уравнений. Для
этого, например, из интегралов задачи формируется их линейная комбинация K = к
У] ^j Vj (x), где ^j = const, и записываются условия стационарности K по переменным j=0
задачи:
Ж
^ = 0 (i = 1,...,n). (13)
dx.
При
д 2K
_ _ = 0 (i, I = 1,... ,n) уравнения (13) совместно с частью уравнений V0(x) = dxidxi
c0,... , Vk(x) = Ck определяют стационарные решения уравнений (12). Если указанное условие не выполняется, то система уравнений (13) будет вырожденной. Согласно [8], в последнем случае уравнения (13) определяют стационарные инвариантные многообразия уравнений (12). Аналогично стационарные решения системы с функцией Лагранжа L = L(qm+1,... , qn, qi,... , qn) и m циклическими интегралами можно найти из уравнений
dR
7—= 0, qLi = 0 (i = m +l,...,n), qi = pl = const (I = l,...,m). dqi
Здесь R — функция Рауса.
4.2. Исследование устойчивости стационарных решений методом Рауса — Ляпунова
Метод Рауса — Ляпунова позволяет не только выделить стационарные решения и ИМ уравнений вида (12), но и исследовать их устойчивость. Интеграл К в этом случае используется в качестве функции Ляпунова.
Запишем вторую вариацию интеграла К в окрестности найденного ранее стационарного решения ж0 = (ж°,... ):
2 ¿2 K = ££ ^ ^ Vj
j=0 г, i=i
Zi Zi, (14)
где ^г — отклонения значений переменных в возмущённом движении от их значений в невозмущённом движении.
Согласно теореме Рауса — Ляпунова, условия знакоопределённости квадратичной формы (14) будут достаточными условиями устойчивости исследуемого решения. Для получения указанных условий применяются критерии знакоопределённости квадратичной формы. Если система не имеет первых интегралов, то для построения функции Ляпунова используются дифференциальные следствия уравнений движения, например, определяющие закон изменения энергии и т. п.
4.3. Исследование устойчивости стационарных решений по первому приближению
Уравнения (5) после разложения в ряд в окрестности стационарного решения qi = q, const, qi = 0 (i = 1,... , n) могут быть представлены в виде
мС + (я + 2С) С + (К + р) С = д(с, С), (15)
где С — столбцы (п х 1) отклонений обобщённых координат, скоростей и ускоре-
ний от их значений в невозмущённом движении; ) — столбец нелинейных по С
слагаемых выше первого порядка малости, ф(0, 0) = 0; М = Мт, Я = Ят, С = — Ст, К = Кт, Р = — Рт — постоянные матрицы п х п соответственно кинетической энергии, диссипативных, гироскопических, потенциальных и неконсервативных позиционных сил. При д((,С) = 0 получаем уравнения первого приближения
М< + 5 С + ^С = 0, где 5 = Я + 2С, W = К + Р. (16)
Характеристическое уравнение для уравнений (16) запишем в виде
ае1 (МЛ2 + 5Л + W) = йоА2га + М2"-1 + ^Л2га-2 + ... + й^-А + = 0. (17)
Анализ корней уравнения (17) позволяет решить вопрос об асимптотической устойчивости или неустойчивости нулевого решения уравнений (15) и соответствующего решения уравнений (5). Для этого применяются критерии Рауса — Гурвица, Льенара — Шипара и др. отсутствия нулевых корней и корней с положительными вещественными частями. Для консервативных механических систем уравнение (17) содержит Л только в чётных степенях (т. е. все Л.2г-1 = 0, г = 1,... , п). Устойчивость в этом случае возможна, если характеристический полином имеет лишь чисто мнимые корни. Критерий [11] позволяет получить условия, при которых корни рассматриваемого полинома относительно Л2
0
x
являются вещественными отрицательными числами. В случае консервативных механических систем уравнения (16) принимают вид М ( + + К ( = 0 и для анализа корней соответствующего характеристического уравнения, кроме перечисленных выше критериев, можно использовать первые интегралы этих линейных уравнений [8, 12].
Опираясь на классические теоремы Кельвина — Четаева [10] о влиянии структуры сил на устойчивость, можно решить вопрос как об устойчивости, так и о стабилизации нулевого решения системы (15).
5. Комплекс программ Analysis
На основе описанных алгоритмов качественного анализа уравнений движения создан пакет расширения Analysis на базе СКА Mathematica. Пакет содержит программы построения уравнений стационарности для нахождения стационарных решений и ИМ уравнений движения по известным первым интегралам задачи и исследования устойчивости найденных таким способом решений уравнений движения методом Рауса — Ляпунова. Кроме того пакет включает следующие программы: исследование устойчивости стационарных решений на основе теорем Ляпунова об устойчивости по первому приближению и Кельвина — Четаева по структуре сил: построение уравнений первого приближения как по известной характеристической функции задачи (Лагранжа, Рауса, Гамильтона), так и по нелинейным уравнениям движения; выделение линейных квадратичных интегралов линейных уравнений движения для исследования спектра системы; получение условий асимптотической устойчивости линейных систем (критерии Рауса — Гурвица, Льенара — Шипара и т.п.), условий знакоопределённости квадратичной формы (критерий Сильвестра, теорема о приводимости двух квадратичных форм к сумме квадратов одновременно и др.); разбиение матриц линейных уравнений движения на симметричные и кососимметричные части для определения структуры сил, действующих на систему.
Специфику применения пакета Analysis продемонстрируем на примерах 3, 4: нахождение ИМ уравнений движения (9) и исследование на устойчивость положений равновесия электрической цепи (см. рис. 4).
Пример 3. Рассмотрим задачу выделения ИМ уравнений (9) посредством редукции, т. е. вначале решим эту задачу для соответствующей раусовой системы. Для выделения ИМ уравнения (10) применим модифицированный метод Рауса — Ляпунова [13], основанный на использовании "расширенной" функции Рауса R = R + pf (в) в, где f (в) = ¿$(в)/сИ, в(в) — некоторая гладкая функция, p = const. Уравнения Рауса для основной и "расширенной" функции Рауса по позиционным координатам задачи, как известно, совпадают и имеют вид (10).
Используя функцию R в качестве исходных данных для пакета Analysis, построим уравнения стационарности этой функции
дК ,2д дК p2cosec2в cte-в , . .
= ml в + pf (в) = 0, — = pf (в)в + ---—— + mgl sin в = 0. (18)
дв дв J у ' ml2 к J
Найдём условие вырожденности системы (18):
p ^cosec^ ^в — f ^f'^)) + mgl sin в = 0.
ml2
Интегрируя данное уравнение, получим выражение для f через параметры задачи:
2 m/2
f= (d - mg/ cos 0) - ctg20 (d = const). (19)
p2
Таким образом, уравнение Рауса (10) имеет два семейства ИМ m/20 + pf^2(0) = 0, где f 1,2(0) — соответствующие решения уравнения (19). В явном виде уравнения этих семейств ИМ записываются в виде
m/20 2 m/2(d - mg/ cos 0) - p2 ctg20 = 0. (20)
Найденные ИМ раусовой системы могут быть "подняты" как инвариантные в фазовое пространство исходной лагранжевой системы. Для этого к уравнениям (20) достаточно добавить соответствующий циклический интеграл. В рассматриваемой задаче указанный интеграл записывается как dL/дф = m/2 sin2 0 ф = p.
Таким образом, уравнения Лагранжа (9) имеют два семейства ИМ:
m/2 sin2 0 ф = p, m/2 0 ± -у/2 m /2 (d - mg/ cos 0) - p2ctg20 = 0.
Пример 4. Для уравнений (11) рассмотрим задачу выделения положений равновесия (т. е. не зависящих от времени t решений данных уравнений) и исследования их устойчивости. Для нахождения решений приравняем к нулю правые части уравнений (11):
U4 + «в - из = 0, /01 + ii + /зЫ = 0, ¿2 - ¿i + ,/2(u4) - lo = 0,
U5 - «4 = 0, /02 - ¿2 + /4(U5) = 0, ¿1 + /0 - fl («б) = 0. (21)
Здесь /i(«6) = giMe(u6 - 3si)2, f («4) = g2«4(«4 - 3S2)2, /3(«з) = дз«з(«3 - 3^з)2, f4(«5) = g4«5(«5 - 3s4)2 — кубические интерполяционные многочлены Лагранжа, аппроксимирующие соответственно функции fi(w6), f2(«4), /з(«з), f4(«5), Sj — первая экстремальная точка функции fj(«), gj = fj(sj)/(4вз).
Исключив с помощью части уравнений (21) переменные «з, «5, ii, i2 из остальных, получим два полиномиальных алгебраических уравнения относительно «4 и «б:
/oi - /0 + gi мб(мв - 3si)2 + gз(м4 + «в)(«4 + «б - 3 вз)2 = 0,
g2 «4(^4 - 3s2)2 - gi «в («в - 3si)2 + g4 «4(^4 - 3S4)2 + /02 = 0. Считая в последних «4 = «4 = const, «б = «б = const, найдем /0i, /02 как функции от
«4, «0, gj:
/0i = /0 - gi^M - 3si)2 - gз(м0 + мб)(«4 + «0 - 3вз)2,
/02 = gi^« - 3si)2 - w4[g2(«0 - 3S2)2 + g4(«4 - 3S4)2]. (22)
Таким образом, при ограничениях (22) уравнения (11) имеют в качестве решения семейство положений равновесия следующего вида:
¿i = giмб(«0 - 3вз)2 - /0, ¿2 = gi «°(«б - 3si)2 - g2 м4(«0 - 3в2)2,
«з = м0 + «б, «4 = «0, «5 = «0, «б = «0. (23)
Используя теоремы Ляпунова об устойчивости по первому приближению, исследуем устойчивость элементов полученного семейства с помощью пакета Analysis, полагая для компактности получаемых далее выражений
s2 = s3 = s4 = 2sb g2 = ga = g4 = gl, L2 = L1, C2 = Ca = C4 = C1. При этих условиях соотношения (22), (23) запишутся соответственно как
I01 = I0 - 18g1s1, I02 = -62g^,
i1 = 2g1 s1 - I0, i2 = -30g1s1, ua = 4s1, u4 = 2s1, u5 = 2s1, u6 = 2s1. (24)
Исходными данными для пакета Analysis при исследовании устойчивости являются уравнения состояния и исследуемое решение (здесь это уравнения (11) и решение (24)). Как и в предыдущих случаях, данные в программу вводятся посредством диалоговой панели либо файла, формируемого пользователем. По входным данным пакет формирует уравнения возмущенного движения в первом приближении
L1Z1 = <4 + Се - <a, C1Za = <1 - 12g1s1Za, ^<5 = -(2,
L1<2 = <5 - <4, C1<4 = <2 - <1, C1<6 = - <1 - 3 g1s1<6 и соответствующее им характеристическое уравнение
4C4L2 A6 + 15pC?L1 A5 + C?L1(20C1 + 9p%)A4 + 60pC?Lx Aa+
+C1(20C1 + 27p2L1)A2 + 45pC1 A + 9p2 = 0, (25)
где <j (j = 1,6) — отклонения в возмущённом движении, p = f1(s1)/s1.
По коэффициентам уравнения (25) выписываются условия существования у данного уравнения корней с отрицательными вещественными частями в виде условий Льенара — Шипара. Приведём только подлежащие проверке условия:
L1(20C1 + 9p2L1) > 0, pL1 > 0, C1(20 C1 + 27p2L1) > 0, C1P > 0,
Aa = C[p2L4(8C1 + 9p%) > 0, A5 = paC19L1(4C1 + 3p2Lx)2 > 0.
Очевидно, данные неравенства всегда выполняются, так как исходя из постановки задачи C1 > 0, L1 > 0, а в силу выбора значений f1(s1), s1 (в первом квадранте координатной плоскости (i, u)) p > 0. Следовательно, по теореме Ляпунова об устойчивости по первому приближению исследуемое решение асимптотически устойчиво.
Заключение
В работе представлены реализованные в виде пакетов расширения системы компьютерной алгебры Mathematica комплексы программ для моделирования и качественного анализа в символьном виде динамических систем, дифференциальные уравнения которых можно получить по некоторой характеристической функции. В качестве таких систем рассмотрены системы твёрдых тел и линейные и нелинейные (из класса "полных") электрические цепи. Движение (состояние) механических систем и линейных электрических цепей описывается функцией Лагранжа, что позволяет применять методы исследования механических систем к линейным электрическим цепям. Нелинейная
электрическая цепь описывается функцией, называемой смешанным потенциалом. Характеристические функции исследуемых систем можно использовать для получения уравнений движения (состояния) и как функции Ляпунова при качественном анализе этих уравнений. Эффективность описанных комплексов программ подтверждается, в частности, их применением при исследовании в символьном виде динамики искусственных спутников (спутника с гиродинами [14], устойчивости относительных равновесий спутника-гиростата на круговой орбите [15]), а также для получения инвариантных многообразий и исследования их устойчивости в задачах Клебша — Тиссерана — Бруна [16]. В настоящее время представленные комплексы программ продолжают развиваться за счёт расширения алгоритмической базы.
Список литературы
[1] Банщиков А.В., Бурлакова Л.А., Иртегов В.Д., Титоренко Т.Н. Программный комплекс LinModel для анализа динамики механических систем большой размерности. Свидетельство о гос. регистрации программы для ЭВМ № 2008610622. ФГУ ФИПС. 1 февраля 2008 г.
[2] Банщиков А.В., Бурлакова Л.А., Иртегов В.Д., Титоренко Т.Н. Программный комплекс для выделения и исследования устойчивости стационарных множеств механических систем. Свидетельство о гос. регистрации программы для ЭВМ № 2011612429. ФГУ ФИПС. 5 июля 2011 г.
[3] Иртегов В.Д., Титоренко Т.Н. Об использовании электромеханических аналогий // Системы поддержки принятия решений для исследования и управления энергетикой. Новосибирск: Наука, 1997. C. 136-145.
[4] Лурье А.И. Аналитическая механика. М.: ФИЗМАТГИЗ, 1961. 824 с.
[5] Синицкий Л.А. Методы аналитической механики в теории электрических цепей. Львов: Вища школа, 1978. 139 с.
[6] Brayton R.K., Moser J.K. A theory of nonlinear networks — 1 // Quarterly of Appl. Math. 1964. Vol. 22, No. 1. P. 1-33.
[7] Ляпунов А.М. О постоянных винтовых движениях твёрдого тела в жидкости. М.: Изд-во АН СССР, 1954. Т. 1. C. 276-319.
[8] Иртегов В.Д. Инвариантные многообразия стационарных движений и их устойчивость. Новосибирск: Наука, 1985. 142 с.
[9] Иртегов В.Д., Титоренко Т.Н. Об инвариантных многообразиях систем с первыми интегралами // Прикл. математика и механика. 2009. Т. 73, вып. 4. С. 531-537.
[10] Четаев Н.Г. Устойчивость движения. Работы по аналитической механике. М.: Изд-во АН СССР, 1962. 535 с.
[11] Кац А.М. К вопросу о критерии апериодической устойчивости // Прикл. математика и механика. 1951. Т. 15, вып. 1. C. 120.
[12] Банщиков А.В., Бурлакова Л.А. Об алгоритмах символьных вычислений при исследовании устойчивости // Программирование. 1997. № 3. C. 72-80.
[13] Иртегов В.Д., Титоренко Т.Н. Об одном методе выделения инвариантных многообразий лагранжевых систем // Современные технологии. Системный анализ. Моделирование. 2012. № 3 (25). C. 30-34.
[14] Банщиков А.В., Члйкин С.В. Моделирование и анализ устойчивости спутника с гиро-динами с помощью программного комплекса LinModel // Там же. 2010. № 4 (22). C. 37-45.
[15] Члйкин С.В., Блнщиков А.В. О гироскопической стабилизации относительных равновесий орбитального осесимметричного гиростата // Математическое моделирование. 2013. Т. 25, № 5. C. 109-122.
[16] Иртегов В.Д., Титоренко Т.Н. Об инвариантных многообразиях в задаче Клебша — Тиссерана — Бруна // Прикл. математика и механика. 2012. Т. 76, вып. 3. C. 374-382.
Поступила в 'редакцию 14 апреля 2014 г.
Symbolic computation in modelling and qualitative analysis of dynamic systems
Banshchikov Andrei V.1, Burlakova Larisa A.1, Irtegov Valentin D.1, Titorenko Tatiana N.1
Purpose. We propose a technology for investigation of dynamics and stability for motion of mechanical systems and electrical circuits in a symbolic form. The software packages that support the proposed technology are described.
Design/methodology/approach. The algorithms implemented in the packages are based on classical methods of analytical mechanics, stability theory of motion and methods representing their further development.
Findings. Algorithms for modelling and qualitative analysis of mechanical systems, linear and nonlinear electrical circuits have been described. These algorithms were implemented using computer algebra system Mathematica and presented in the form of its add-on packages. The developed software allows constructing of the characteristic function for the above systems on the base of their geometric description. The function can be used both for deriving equations of motion (state equations of circuit) and for qualitative analysis of stationary solutions (invariant manifolds) of these equations. Functional abilities of the packages have been demonstrated by the specific examples.
Research limitations/implications. Software packages are applied for solution of the problems mentioned above entirely in a symbolical form. Such technique is efficient enough for construction of mathematical models (differential equations) and for their qualitative analysis, in particular, for parametric analysis.
Originality/value. Using the developed software, new results have been obtained in investigation of the following problems: dynamics of artificial satellites, motion of a rigid body in an ideal fluid, Euler's equations on Lie algebras, electrical nonlinear circuits and some others. Packages continue to evolve through the expansion of the algorithmic framework. Keywords: modelling, qualitative analysis, mechanical system, electrical circuit, computer algebra system, complex of programs.
Received 14 April 2014
1 Institute for System Dynamics and Control Theory SB RAS, Irkutsk, 664033, Russia Corresponding author: Banshchikov Andrei V., e-mail: bav@icc.ru
Aknowlegements: This research was supported by the Presidium of the Russian Academy of Sciences (basic research program No. 17.1) and by the President Grant for Governmental Support of the Leading Scientific Schools (grant No. 5007.2014.9).
