Алгоритмы построения границ областей устойчивости линейных гамильтоновых систем с помощью пакета Matlab Текст научной статьи по специальности «Математика»

УДК 517.923

DOI 10.25559/SITITO.2017.4.525

Юмагулов М.Г., Белова А.С.

Башкирский государственный университет, г. Уфа, Россия

АЛГОРИТМЫ ПОСТРОЕНИЯ ГРАНИЦ ОБЛАСТЕЙ УСТОЙЧИВОСТИ ЛИНЕЙНЫХ ГАМИЛЬТОНОВЫХ СИСТЕМ С ПОМОЩЬЮ ПАКЕТА МА^АВ

Аннотация

В статье приводятся основные этапы алгоритма построения областей устойчивости динамических систем, описываемых линейной гамильтоновой системой вида х' = А{1,а,р)х. Алгоритм основан на методах теории нелинейных колебаний исследования устойчивости стационарных решений линейных дифференциальных уравнений с периодическими коэффициентами, зависящих от малого параметра. Алгоритм реализован с помощью математического пакета Matlab. В качестве приложения рассмотрена задача построения области устойчивости треугольных точек либрации плоской ограниченной эллиптической задачи трех тел.

Ключевые слова

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

Yumagulov M.G., Belova A.S.

Bashkir State University, Ufa, Russia

THE ALGORITHMS FOR CONSTRUCTING THE BOUNDARIES OF THE STABILITY REGIONS OF LINEAR HAMILTONIAN SYSTEMS BY USING MATLAB

Abstract

The article presents the main stages of the algorithm for constructing the stability regions of dynamical systems described by a linear Hamiltonian system of the form x' = A(t,a,p)x. The algorithm is based on the methods of the theory of nonlinear oscillations of stability studies of stationary solutions of linear differential equations with periodic coefficients depending on a small parameter. The algorithm is implemented by using Matlab CAS. As an application, we have solved the problem of constructing the stability region of libration points of a flat elliptic restricted three-body problem is analyzed in detail.

Keywords

Hamiltonian systems; stability; the stability region; periodic solutions; parameter; the three-body problem.

Введение

В статье рассматривается линейная гамильтонова система

= а, в)х, ХЕК2« (1)

где а, в — скалярные параметры, матрица А(Х а, в) является непрерывной и Т -периодической по t : + Т, а, в) = А(Х а, в) . Уравнение (1) при всех значениях параметров а и в имеет точку равновесия х = 0, которая при одних значениях параметров может быть устойчивой по Ляпунову, а при других — неустойчивой.

Пусть все мультипликаторы системы

^ = Ааа0,£0)х, ХЕК™ (2) по модулю равны 1. Возможны два основных случая:

1) все мультипликаторы являются простыми и при этом ни один из них не равен 1 либо -1;

2) имеются кратные мультипликаторы (кратности 2) равные 1 либо -1, остальные являются простыми.

В первом случае при всех значениях (а, в) близких к (а0,^о) система (1) является устойчивой. Во втором случае происходит изменение характера устойчивости в

окрестности точки (а0,/30) . А именно, как правило, в плоскости параметров (а, в) через точку (а0,р0) проходят одна или несколько гладких кривых, которые образуют границу областей устойчивости и неустойчивости (см. Рис. 1).

Ниже через П будем обозначать плоскость параметров (а, в). Множество ОсП в плоскости параметров ( а, в) будем называть областью устойчивости системы (1), если для любых (а, в) £ О точка равновесия х = 0 системы (1) является устойчивой в линейном приближении, а для любых ( а, в) £ F = П\О эти точки неустойчивы в линейном приближении. При этом множество F будем называть областью неустойчивости системы (1). Точку (а, в) будем называть граничной точкой области устойчивости системы (1) если в каждой ее окрестности имеются точки из О и F. Множество Г граничных точек будем называть границей области устойчивости системы (1). Обратим внимание на то, что в приведённых здесь понятиях термины устойчивость и неустойчивость понимаются как устойчивость и неустойчивость точек равновесия в линейном приближении.

(ао,Ро)

Задачи о построении и изучении границ областей устойчивости, об исследовании поведения динамических систем при переходе параметров через эти границы являются важными и интересными задачами теории управления и регулирования, теории динамических систем, теории нелинейных колебаний и их многочисленных приложений. Исследованию таких задач посвящена обширная литература; здесь предложен ряд эффективных методов исследования, решен ряд важных с теоретической и практической точек зрения задач (см. [1] - [5] и имеющуюся там библиографию). Исследования активно продолжаются в различных направлениях (см., например, [6] - [8]).

Вопросы построения границ областей устойчивости системы (1) являются сложными

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

В настоящей работе предлагается новый алгоритм определения границ областей устойчивости системы (1), основанный на модификации метода М.Розо [8] и формулах теории возмущений линейных операторов. В качестве приложения рассмотрена задача о построении границ областей устойчивости точек либрации плоской ограниченной эллиптической задачи трех тел. Предлагаемая схема использует систему компьютерной математики МаЙаЬ.

Алгоритм построения области устойчивости

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

Приведем основные этапы алгоритма построения границ области устойчивости точки равновесия х = 0 системы (1). Для простоты ограничимся рассмотрением ситуации, когда система (1) имеет вид:

dx

— = (Л0 + (а - а0)В1 (0 + Св -^0)^2С£))х ,

х £ М2М (3)

где А0 - постоянная матрица, В1(Ь) и В2(0 - непрерывные и Т-периодические матрицы. Общий случай может быть рассмотрен по той же схеме, но требует более громоздких построений. Пусть матрица А0 имеет пару

собственных значений — , а остальные ее

т

собственные значения также являются чисто мнимыми, но не кратны числу ^. Ниже для

простоты будем считать, что Т=2п.

Граничную кривую системы (3) будем строить как функцию в = Ро + КЮ , где £ = а — а0 и ^£) = ух£ + Ф(£) и Ф(£) = 0(£2) при £ ^ 0. Приводимая ниже схема построения границ областей устойчивости системы (1) позволяет получить более полные формулы для представления функций /(е), а именно, вида: /(£) = У1£ + Уг£2 + —+ У„£п + 0(£п+1), при этом указывается алгоритм вычисления значений коэффициентов ух, у2 ,...,уп.

Этап 1. Произведем в системе (3) замену в = р0 + ^е) и а = а0 + £ , которая позволит перейти к дифференциальному уравнению, зависящему от малого параметра £:

^ = (Д0 + еА1(ь))Х, ХЕК2

(4)

где Дх(0 - непрерывная 2 п -периодическая матрица.

Свойства устойчивости решения х = 0 уравнения (3) зависит от мультипликаторов матрицы А0 + еА1(1) . К сожалению, явное построение мультипликаторов матрицы Д0 + еА^ (£) возможно лишь в самых простых случаях; исследованию этого вопроса посвящены работы многих авторов [1, 2].

Этап 2. Одним из эффективных способов исследования устойчивости решения х = 0 уравнения (3) является предложения М.Розо [8] схема перехода к равносильному уравнению вида

у' = (До + Е^ + Е^а^у, где матрица - некоторая специально конструируемая постоянная матрица. А именно, матрица является решением матричного уравнения

г2п г2п

I е~А°тЗ1еА°тйт = I е"л°тД0ел°тйт (5) -'о -'о

Известно, что уравнение (5) имеет единственное решение тогда и только тогда, когда для матрицы А0 выполняется условие отсутствия Т - резонанса т.е. любая пара А; и А| различных собственных значений матрицы А0 удовлетворяет соотношению Я^ —

2 па1 .

А) Ф ——, I - мнимая единица, ^ — целое число .

На втором этапе выясняется выполнено ли для матрицы А0 условие отсутствия Т -резонанса. Если оно не выполняется, то предлагается алгоритм перехода к равносильному уравнению, для которого это условие выполнено.

Этап 3. Решается матричное уравнение (5). Его решение зависит от параметра у1 , определяющее искомую границу (точнее, касательную к границе) области устойчивости. Пусть З*^) - это решение уравнения (5).

Этап 4. Найдём собственные значения матрицы З*^) . Тогда искомыми будут те значения у1, при которых эта матрица будет иметь нулевые собственные значения.

Приложение: плоская ограниченная эллиптическая задача трех тел

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

массами т0 и тх : т2 << тх < т0 . Движение третьего тела Р2 в указанной задаче в координатах Нехвилла (^,п) описывается неавтономной системой с 2п -периодической правой частью вида (см., например, [4 — 6]):

¡У-2п' = р(*-ц+ ^

(?2+Л2)3/2 ((?-1)2+Л2)3/2/'

п

2? = р("

П +

м- л

Ц =

(ц-1)л___

(?2+л2)3/2 ((?-1)2+л2)3/2

(6)

эксцентриситет

где р =-, , ,

кеплеровской орбиты (0 < £ < 1), t — истинная аномалия, т0 и тх - массы активно гравитирующих тел (0 < тх < т0), ц - параметр масс (0 < ц < 1) . Штрихами обозначены производные по t . Система (6) является гамильтоновой.

Система (6) имеет пять постоянных решений - точек либрации: прямолинейных L1, L2, Lз и треугольных L4,L5. Треугольные точки

либрации: L4 ; у), L5 ф -у) могут быть

как устойчивыми, так и неустойчивыми. Нас будет интересовать зависимость свойства устойчивости решений от параметров ц и £ в окрестности точки либрации L4.

Вопросам построения областей устойчивости треугольных точек либрации посвящены многочисленные исследования. Известные здесь наиболее полные результаты были получены во второй половине прошлого столетия (см. [6] и имеющуюся там библиографию). Исследования в указанном направлении активно продолжаются (см., например, [6, 7, 10 — 13]). Приведем схему построения границ областей устойчивости треугольных точек либрации на основе вышеприведенного алгоритма с использованием пакета МаЙаЬ. Приводимые ниже формулы получены с помощью аппарата символьных вычислений этого пакета, т.е. являются точными.

Отметим, что граница области устойчивости треугольных точек либрации в плоскости параметров (и, е) образована тремя

непрерывными линиями Гх, Г2 и Г3. Кривые Гх, Г2 начинаются в точке (и0,0), а кривая Г3 - в точке

О*, 0); здесь ^ = ± - ^ = 0,028595 ... , ц* = ± -

— = 0,038520.

18

Переход к основному уравнению

Рассматриваемая задача может быть (см. [12]) сведена к изучению гамильтонова уравнения

^ = А(£,цДЖ h Е К4, (7)

где

A(E,M) =

3V—

—р(1-2ц)

3V—

Р(1 - 2ц) 0

4 Р

-2 0

При этом матрица А(£, цД) представима в виде: А(£, цД) = А0(ц) + (-scost + £2 cos21 — -)А1(ц), где

0 0 10-1

0 0 0 1

3 3V— А„(ц) = т — (1-2ц)

Ах(ц) =

3V—

—(1-2ц)

0 0

3

4

3V—

L-4 (1 -

9

4

0 2

-2 0 0 0

3V—

—(1-2ц)

9

4

0 0

0 0

Таким образом мы перейдем к уравнению вида

h' = A0(^)h + (—scost + £2 cos2 t)Ax(^)h

+ £3А3(£,цД), heM4 (8)

А3(£,цД) = -:

-А1(ц), матрица А3 является

l + £COSt

2п-периодичной и непрерывной по t. Система (8) является системой вида (3) и поэтому для нее применим вышеприведенный алгоритм. Ограничимся приведением лишь основных результатов для кривой Г3.

Построение кривой Г3

Покажем, что уравнение (8) определяет в точности одну гладкую кривую Г3 , начинающуюся в точке (ц*,0) и состоящей из граничных точек области устойчивости системы (8).

Зафиксируем число ф, 0 < ф < п и рассмотрим прямую, заданную уравнениями в параметрической форме: ц = ц* + Scos ф , £ = 5sinф, где S - вспомогательный параметр. В плоскости (ц, £) эта прямая является касательной к кривой Г3 и проходит через точку (ц*,0), образуя с осью ц угол ф. Опишем этапы построения этой касательной.

Шаг 1. Подставим уравнения указанной прямой в (8) и перейдём к зависящему от малого параметра S и угла ф уравнению

h' = (А0 + 5PX (t, ф) + 52P2(t,ф,£2,ц2))h + 53P3(t,5, ф, £2, ц2Ж hel4 (9)

где Px(t, ф) = cosфB0 - sinф cost Ах ,

P2 (t, ф, £х, цх) = —sinф cosф costB0 + sin2 ф cos21 Ах, а матрица P3 (t, S, ф) является 2п-периодичной по t и непрерывной. Подсчет показывает, что верны равенства:

А0=А0(ц*) =

Ах=Ах(ц*) =

B0 =

0 0

3

4

V2—

0 0

3

4

V2—

0 0

0

-3V—

0 0

V2—

9

4 0

0

V2—

9

4 0

0

-3V— 2

0

0

-2 0

0 0 0 0

0 0

0 0

0 0 0 0

0 0 0 0

Q =

Шаг 2. Найдем собственные вектора матрицы А0 , соответствующие ее неполупростыми

собственным значениям: Л12 = -р, А34 = --р.

, V2 , V2

Заметим, что перечисленные собственные значения удовлетворяют соотношению A¡ - A¡ Ф qi при q e Ъ, поэтому 2п-резонанс отсутствует.

Шаг 3. Построим матрицу Q из собственных векторов.

i(-8-iV46) ¿(32iV2 - 18V2—) i(-8 + iV46) ^(-32iV2 - 18V2—)"

iV2 2 -iV2 2

i(4iV2-V2—) ^(-8 + 4iV46) i(-4iV2-V2—) i(-8-4iV46)

1 0 1 0 -Произведём замену h = Qy и перейдем от (9) к равносильному уравнению

y' = (Л0 + 5РХ (t, ф) + 52Р2 (t, ф,£2,ц2) +

53P3(t,5,9,£2^2))y,yel4 (10)

где Px(t, ф) = cosфB0 - sinф cost^ ;

Р2 (t, ф, £х, цх) = —sinф cosф costB0 + sin2 ф cos2 tÁ^

Á0 = Q"1A0Q =

Ai = Q"1A1Q

/i(23V2 + 2iV23) 3(19 + 3iV46) 20 50

-13iV2 + 2V2—

(- "¡/Z 1 0 0^

0 -¡/Z 0 0

0 0 ¡/Z 1

\ 0 0 0

1

i(V2 + 6iV23) -i(-39i + 7V46)\

20 50

2 + iV46 3(19 - 3iV46)

-i(V2 - 6iV23) i(39i + 7V46) -i(23V2 - 2iV23)

20

(2 - iV46)

50

-31iV2 - 6V2—

20 1

= Q B„Q 3 20

í t3iV3(2i + V46)

3V69 8

7(iV6 +V69)

__V-(4 - 7iV46)

50

1

13iV2 + 2V2— 20

3

3 3

—V—(18 + iV46) —(-41W6 + 9V69)

v--l2(-8i + V46) 40V-(36 + 7i^46)

2jV3(18 - iV46) 5-j(41iV6 + 9V69)^

--j2(8i+V46) 40V-(36 - 7i^46)

3 27

-—¡V—(-2i + V46) —(-iV6 + V69)

3V69 8

40V— (4+7iV46) j

Шаг 4. Метод М. Розо позволяет с помощью замены у = (I - с 2п -периодической

матрицей Нх00 перейти от (10) к системе

z' = (Л0 + 5З1(ф) + 52З2(ф,ц2,£2) +

5=%а,ф,ц2,£2)^ ^ЕК4 (11)

с постоянными матрицами и

непрерывной 2п -периодической матрицей 53. Свойства устойчивости системы (11) (а значит, и треугольных точек либрации в линейном приближении исходной системы (8)) можно определить с помощью постоянной матрицы + 5З1(ф) + 5%(ф, ц2,£2). Шаг 5. В соответствии с методом М. Розо

рассмотрим матричное уравнение

г2п г2п

5-SotP1(t,9)eSotdt

г2п г2

I e-s°tS1es°tdt = I

-'о -'о

(12)

По решению этой системы построим матрицу (ненулевые элементы которой совпадают с соответствующими элементами решения системы (12)).

Шаг 6. Найдём собственные значения матрицы S*^):

^1,2,3,4 =

-63 cos2 ф + 32V69 cosф sinф

Шаг 7. Определим, при каком значении ф полученные собственные значения будут равняться 0. Отсюда получим значение ф =

Таким образом, касательная к кривой Г3 проходит через (ц*, 0) и образует угол 90° с осью Оц .

По аналогичной схеме могут быть найдены и последующие приближения к кривой Г3 . В частности, приведенный алгоритм позволяет установить, что кривая кривой Г3 описывается функцией £ = /3 (5) , где S = ^ — и f3 (5) = ZiV5 + ф3(5); здесь

3(2^со5ф-Ы138со5ф) 2 7i-^совф+2 7-^бЭсовф+ 1005Шф 20 50

З-^бЭсовф Зi(4iVЗcosф+7V138cosф)

SK9) =

8 о

40

о о

?1 =

\

621

— = 3,529863 ..., 4

а нелинейность Фз(5) удовлетворяет соотношению: ф3(5) = 0(5) при 5^0.

На Рис. 2 приведен вид кривых, образующих границы области устойчивости треугольных точек либрации системы (8). Изображенные кривые получены в соответствии с предложенным алгоритмом с использованием формул, содержащих вторые приближения.

Рис.2. Граница области устойчивости треугольных точек либрации

Заключение

В статье были приведены основные этапы нового алгоритма построения границ областей устойчивости для линейных гамильтоновых систем вида х' = А(Ь,а,р)х . Разработанный алгоритм реализован с помощью системы компьютерной математики МаЙаЬ. В качестве приложения рассмотрена задача о построении областей устойчивости в линейном приближении треугольных точек либрации плоской ограниченной эллиптической задачи трех тел.

0

0

о о

3(2^со5ф+Ы138со5ф) — 2 7i-^совф+2 7-^бЭсовф+1005Шф

20

З-^бЭсовф

31(—41-^со5ф+7-^138со5ф)

(12)

Литература

Арнольд В.И., Козлов В.В., Нейштадт А.И. Геометрические методы в теории обыкновенных дифференциальных уравнений. - Ижевск: НИЦ Регулярная и хаотическая динамика, 2000. 400 с.

Зигель К.Л., Мозер Ю.К. Лекции по небесной механике - Ижевск: НИЦ Регулярная и хаотическая динамика, 2001. - 384 с.

Ибрагимова Л.С., Мустафина И.Ж., Юмагулов М.Г. Асимптотические формулы в задаче построения областей

гиперболичности и устойчивости динамических систем. // УМЖ. 2016. - №3. - С.59 - 81.

Като Т. Теория возмущений линейных операторов. - М.: Мир, 1975. 740 с.

Каток А.Б., Хасселблат Б. Введение в теорию динамических систем. М.: МЦНМО. 2005. 464 с.

Маркеев А.П. Точки либраций в небесной механике и космодинамике. М.: Наука, 1978. 312 с.

Маршал К. Задача трех тел. М.-Ижевск: Институт компьютерных исследований. 2004. 640 с.

Розо М. Нелинейные колебания и теория устойчивости. М.: Наука. Глав. ред. физ.-мат. лит., 1971. 288 с.

Чезари Л. Асимптотическое поведение и устойчивость решений обыкновенных дифференциальных уравнений. - М.: Мир,

1964. - 477 с.

10. Юмагулов М.Г., Беликова О.Н. Бифуркации периодических решений в окрестностях треугольных точек либрации задачи трех тел. // Известия высших учебных заведений. Математика. 2010. № 6. С. 82-89.

11. Юмагулов М.Г., Беликова О.Н. Бифуркация 4-п периодических решений плоской ограниченной эллиптической задачи трех тел. // Астрономический журнал. 2009. Т. 86, № 2. C. 170-174.

12. Юмагулов М.Г., Сухоруков А. В. Алгоритмы построения границ областей устойчивости точек либрации в задаче трех тел с помощью пакета Maple.// Современные информационные технологии и ИТ-образование, 2016. Т. 12, № 4. C. 181-188.

13. Kovacs T. Stability chart of the triangular points in the elliptic restricted problem of three bodies. - Mon. Not. R. Astron. Soc. 2013. V. 430. Issue 4. Pp. 2755-2760.

14. Chiang H.:D., Alberto L.:F. C.,. Stability regions of nonlinear dynamical systems : theory, estimation, and applications. - Cambridge University Press. 2015. 484 p.

References

1. Arnol'd V.:I., Kozlov V.:V., Nejshtadt A.:I. Geometricheskie metody v teorii obyknovennyh differencial'nyh uravnenij. - Izhevsk: NIC Reguljarnaja i haoticheskaja dinamika, 2000. 400 s.

2. Zigel' K.:L., Mozer Ju.:K. Lekcii po nebesnoj mehanike - Izhevsk: NIC Reguljarnaja i haoticheskaja dinamika, 2001. - 384 s. Associate-

3. Ibragimova L.:S., Mustafina I.:Zh., Jumagulov M.:G. Asimptoticheskie formuly v zadache postroenija oblastej giperbolichnosti i ustojchivosti dinamicheskih sistem. // UMZh. 2016. - №3.- S.59 - 81.

4. Kato T. Teorija vozmushhenij linejnyh operatorov. - M.: Mir, 1975. 740 c.

5. Katok A.:B., Hasselblat B. Vvedenie v teoriju dinamicheskih sistem. M.: MCNMO. 2005. 464 c.

6. Markeev A.:P. Tochki libracij v nebesnoj mehanike i kosmodinamike. M.: Nauka, 1978. 312 c.

7. Marshal K. Zadacha treh tel. M.-Izhevsk: Institut komp'juternyh issledovanij. 2004. 640 s.

8. Rozo M. Nelinejnye kolebanija i teorija ustojchivosti. M.: Nauka. Glav. red. fiz.-mat. lit., 1971. 288 c.

9. Chezari L. Asimptoticheskoe povedenie i ustojchivost' reshenij obyknovennyh differencial'nyh uravnenij. - M.: Mir, 1964. - 477 s.

10. Yumagulov M.G., Belikova O.N. Bifurkatsiya 4n-periodicheskikh resheniy ploskoy ogranichennoy ellipticheskoy zadachi trekh tel. // Astronomicheskiy zhurnal. 2009. T. 86, № 2. C. 170-174.

11. Yumagulov M.G., Belikova O.N. Bifurkatsii periodicheskikh resheniy v okrestnostyakh treugol'nykh tochek libratsii zadachi trekh tel. // Izvestiya vysshikh uchebnykh zavedeniy. Matematika. 2010. № 6. S. 82-89.

12. Yumagulov M.G., Su^orukov A. V. Algoritmy postroenija granic oblastej ustojchivosti tochek libracii v zadache treh tel s pomoshh'ju paketa Maple.// Sovremennye informacionnye tehnologii i IT-obrazovanie, 2016. T. 12, № 4. S. 181-188.

13. Kovacs T. Stability chart of the triangular points in the elliptic restricted problem of three bodies. - Mon. Not. R. Astron. Soc. 2013. V. 430. Issue 4. Pp. 2755-2760.

14. Chiang H.:D., Alberto L.:F. C.,. Stability regions of nonlinear dynamical systems : theory, estimation, and applications. - Cambridge University Press. 2015. 484 p.

Поступила: 15.10.2017

Об авторах:

Юмагулов Марат Гаязович, доктор физико-математических наук, профессор, заведующий кафедрой дифференциальных уравнений факультета математики и информационных технологиий, Башкирский государственный университет, yum mg@mail.ru

Белова Анна Сергеевна, магистрант факультета математики и информационных технологий, Башкирский государственный университет, 89177662488@yandex.ru

Note on the authors:

Yumagulov Marat G., Doctor of Physical and Mathematical Sciences, Professor, Head of the Department of Differential Equations of the Faculty of Mathematics and Information Technology, Bashkir State University, yum mg@mail.ru

Belova Anna S., graduate of the Faculty of Mathematics and Information Technology, Bashkir State University, 89177662488@yandex.ru