научная статья по теме СИМВОЛЬНО-ЧИСЛЕННЫЙ АНАЛИЗ ОГРАНИЧЕННОЙ ЗАДАЧИ ПЯТИ ТЕЛ СРЕДСТВАМИ КОМПЬЮТЕРНОЙ АЛГЕБРЫ Математика

Текст научной статьи на тему «СИМВОЛЬНО-ЧИСЛЕННЫЙ АНАЛИЗ ОГРАНИЧЕННОЙ ЗАДАЧИ ПЯТИ ТЕЛ СРЕДСТВАМИ КОМПЬЮТЕРНОЙ АЛГЕБРЫ»

ПРОГРАММИРОВАНИЕ, 2014, No 3, с. 44-48

КОМПЬЮТЕРНАЯ АЛГЕБРА -

УДК 517.9+519.6+531.39

СИМВОЛЬНО-ЧИСЛЕННЫЙ АНАЛИЗ ОГРАНИЧЕННОЙ ЗАДАЧИ ПЯТИ ТЕЛ СРЕДСТВАМИ КОМПЬЮТЕРНОЙ

АЛГЕБРЫ

© 2014 г. Д. А. Будько, С. А. Щерба

Брестский государственный университет им. A.C. Пушкина, Беларусь 224016 Брест, бул. Космонавтов, 21 E-mail: Dzmitry.Budzko@gmail.com Поступила в редакцию 05.09.2013

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

1. ВВЕДЕНИЕ

Со стремительным развитием систем компьютерной алгебры возрастает и количество ее приложений в самых разнообразных областях науки [1, 2]. Удачное комбинирование символьных и численных методов исследований, а также применение новейших технологий динамической визуализации позволяет добиться новых результатов, затрачивая на это меньшее время. В данной работе использование системы компьютерной алгебры Mathematica позволяет упростить исследование нелинейной системы алгебраических уравнений, содержащей параметры.

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

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

Г 2 1 е<?1 = Г ==(ТТо^ + Т + + ^2

(1 + (а - b)2)3/2'

2__а(1 + 2^i) + (а - b)ß2 +

аШ (1 + а2)3/2 + (1 + (а - b)2)3/2 + b2 ,

ЬШ 2 == 2а^

(1 + а2)3/2 2(а - b)ßi + 1 + ß2

(1 + (а - Ь)2)3/2 Ь2 где параметры ^1,^2 суть отношения между массами тел, выраженные формулами

Ш1 = ^1Ш0, Ш2 = ^2^0-

Геометрическая модель дельтоида представлена на рисунке 1, где роль параметров играют а, Ь, соответственно равные ординатам Ш1 и

Ш2- Параметр ш2 связан с гравитационной постоянной, массой шо и расстоянием между телами с массами Ш1, половина которого принята за единицу длины. Он определяет угловую скорость вращения дельтоида вокруг оси, проходящей через центр масс системы.

Анализ возможных центральных конфигураций для четырех тел, а также доказательство существования центральной конфигурации в виде рассматриваемого дельтоида приведено в [3,

41-

Вторая подсистема состоит из двух уравнений относительно неизвестных ж, у, содержит те же пять параметров ^1, ^2, а, Ь, ш и имеет вид

eq2

+

+

(x2 + y2)3/2

-1 + x

+

((-1+ x)2 + (а - y)2)3/2

+

1 + x

((1+ x)2 + (а - y)2)3/2

^1 +

+

x^2

- xw2 == 0,

(x2 + (b - y)2)3/2

+

(x2 + y2)3/2 V(1 + а2)3/2

+

+

-а + y

((-1 + x)2 + (а - y)2)3/2

+

+

-а + y

((1+ x)2 + (а - y)2)3/2

^1 +

1

To +

-b + y

b2 ((x2 + (b - y)2)3/2

^2 - yw == 0.

Вторая подсистема е^2 определяет относительные положения равновесия ж, у в ограниченной задаче пяти тел, сформулированной на основе дельтоида. Стоит отметить, что сами уравнения ед1 и е^2 выведены с помощью системы Ма(-ке.таЫса путём серии подстановок и правил замены. За основу взяты уравнения движения в пространстве Нехвила в относительных координатах для п тел [5].

Целью данной работы является общий анализ системы: определение области существования решений на плоскости ОаЬ, выявление количества решений в зависимости от значений

Рис. 1. Конфигурация выпуклого (а < b) дельтоида при а = 1.3, b = 2.1

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

2. ОБЩИЙ АНАЛИЗ СИСТЕМЫ УРАВНЕНИЙ

Заметим, что параметры ^i,^2,w2 входят линейно в уравнения eq1, поэтому выразим их а, b

so11 = Simplify [Solve [eq31/. temp ,

temp}] [[1]]/. temp ^ w2] .

Далее, принимая во внимание то, что массовые параметры ^1,^2 должны быть положительными, имеем два неравенства

{^1 > 0,^2 > 0}/. so11.

Решения этих неравенств мы не выписываем здесь ввиду их громоздкости, а лишь проиллюстрируем их решение (область ABCF на рисунке 2) с помощью следующей команды:

RegionPlot [а < b && 0 < && 0 <

/.so11, {а, 0,2.5}, {b, 1,3.5}, MаxRecursion ^ 5, PlotPoints ^ 40, PlotStyle ^ White,

x

y

Boundary Style ^ {Thick, Gray}].

Также на рисунке 2 аналогичным образом изобразим более узкую область ADF серым цветом, где выполняются неравенства 0 < Ц1 < 1,0 < Ц2 < 1 поскольку с точки зрения приложений представляют интерес именно небольшие значения параметров

Рис. 2. Области существования решений на плоскости Oab щи а < Ь. Бифуркационные кривые GE и AE разделяют области с разным количеством решений.

Подставим полученные выражения для

,ц.2,ш2 в уравнения eq2. В результате получится новая система двух уравнений eq3 относительно x,y, которая содержит только два параметра а, Ъ:

eq3 = eq2/.sol1.

Далее представим один из вариантов программы, которая помогает исследовать зависимость количества решений системы

eq3 а, b

области, изображенной на рисунке 2.

Manipulate[ContourPlot[eq3/.a ^ аа

/Ьь ^ bb // Evaluate, {x, -2.5,2.5}, {y, -2,3}, LabelStyle ^ Bold, PlotPoints ^ 20, ContourStyle ^ {{Thick, Black, Dashed},

{Thick, Black}} , MaxRecursion —> 3, Epilog ^ {{PointSize[0.025], (Point[x,y]

/.(setO = {}; Do[soltemp = FindRoot[eq3/.a ^ aa/.b ^ bb // Rationalize[#, 0]&, {{x,x0}, {y,y0}}] // Chop If [!MemberQ [setO, SetPrecision[soltemp, 10]]

AppendTo setO, SetPrecision[soltemp, 10] {x0, -1.5, 1.5, step}, {y0, -0.8, 2.2, step}]; temp = setO))}, Line[{{0, 0}, {1,aa}, {0,bb}, {—1,aa}, {0, 0}}]}], Style ["Given geometric parameters", Bold, Medium], Delimiter,

I {aa, 1.3," a"}, 0.6, 0.8}, {{bb, 2.1, "b"}, 1.3,2.9},

Delimiter, Style["Calculated mass

parameters", Bold, Medium], Delimiter,

Row[{Panel[Dynamic[^1/.sol1/.a ^ aa

/.b ^ bb]," ц1 "], Panel[Dynamic[^2/.sol1

/.a ^ aa/.b ^ bb]," ц "]}, " " ],

Delimiter, Style ["Number of solutions",

Bold, Medium], Delimiter,

Panel[Dynamic[Length[temp]]],

ControlPlacement ^ Left].

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

eq3

виде Epilog функции ContourPlot. Нужно учесть, что система уравнений нелинейная, и заранее неизвестно, сколько решений будет иметь система при разных значениях параметров a, b

Oxy

step

начальные приближения для функции FindRoot. Чтобы исключить попадание одного и того же решения (с разных начальных приближений) в setO

принадлежности решения множеству решений. Оставшийся код настраивает внешний вид объекта Manipulate. Результат работы представленной программы изображен на рисунке 3.

Рис. 3. Зависимость количества решений от параметров a, Ь.

Заметим, что каждое уравнение системы eq3 определяет некоторую кривую на плоскости Oab. Используя встроенные функции Manipulate и ContourPlot системы Mathematica, можно изменять положение ползунков a, b и наблюдать видоизменение кривых, изображенных соответственно сплошной и штриховой линиями на рисунке 3. При этом координаты x и y каждой точки пересечения этих кривых определяют решение системы eq3. Как видно из рисунка 3, имеется 13 точек пересечения кривых для значений параметров a = 1.3, b = 2.1. Отметим, a, b

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

Для придания наглядности код программы был немного упрощен. Удалена проверка, когда метод решения нелинейных уравнений, используемый в FindRoot не сходится, но все равно выдает последнее приближение. Также удалены настройки контроля точности решений Working-Precision, а также максимального количества итераций Maxlterations. При нахождении eq3

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

При этом первостепенной представляется

Oab

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

Jacob = D [eq3 [[1,1]], x] D [eq3 [[2,1]], y]

-D[eq3[[1,1]],y] D[eq3[[2,1]],x] ==0.

Jacob

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

xy

Jacob eq3

Сама задача о нахождении таких параметров a, b

кривым, реализовано в данной работе, используя алгоритмы, приведенные в [6]. Проведенный

eq3

могут иметь 11, 13 и 15 решений. Результаты приведены на рисунке 2, где на плоскости параметров построены две бифуркационные кривые, которые отделяют области с разным

a, b

из области GDE на рисунке 2 имеется ровно

13 решений, для параметров a, b из области AGE - 11 решений, для параметров a,b из области AEF - 15 решений.

СПИСОК ЛИТЕРАТУРЫ

1. Акритас А. Основы компьютерной алгебры с приложениями. М.: Мир, 1994. 544 с.

2. Таранчук В.Б. Основные функции систем компьютерной алгебры. Минск: БГУ, 2013. 59 с.

3. Grebenikov Е.А., Ikhsanov E.V., Prokopenya A.N. Numeric-Symbolic Computations in the Study of Central Configurations in the Planar Newtonian Four-Body Problem. In: Computer Algebra in Scientific Computing, Ganzha V.G., Mayr E.W.,

Vorozhtsov E.V. (eds.), LNCS. V. 4194. Springer, Heidelberg, 2006. Р. 192-204.

4. Прокопе

Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.

Показать целиком