Алгоритм метода дополнительных переменных Текст научной статьи по специальности «Математика»

Сер. 10. 2012. Вып. 2

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

ПРИКЛАДНАЯ МАТЕМАТИКА

УДК 517.9:519.6

Л. К. Бабаджанянц, К. М. Брэгман

АЛГОРИТМ МЕТОДА ДОПОЛНИТЕЛЬНЫХ ПЕРЕМЕННЫХ

Предисловие. Рассматриваются полные системы дифференциальных уравнений в частных производных (и среди них системы обыкновенных дифференциальных уравнений - ОДУ), правые части которых можно записать при помощи четырех действий алгебры и любых допустимых суперпозиций функций конечного числа аргументов, принадлежащих неограниченно пополняемому набору функций, называемому библиотекой. Предлагается основанный на методе дополнительных переменных (МДП) автоматизированный алгоритм их сведения к полиномиальным автономным системам, т. е. к системам с полиномиальными по неизвестным правыми частями. Библиотеку можно пополнять любыми функциями, которые удовлетворяют полным, не обязательно автономным, полиномиальным системам. Например, библиотека может содержать все элементарные функции и очень многие специальные функции, используемые в приложениях. Алгоритм разработан на основе результатов статьи [1] с учетом возможностей пакета «Mathematica»[2], которые позволили наделить его следующими дополнительными свойствами:

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

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

- конечные уравнения и данные получаются в аналитической форме.

Алгоритм и пример его применения описываются в п. 2, а в п. 1 приведены необходимые предварительные сведения.

1. Введение в метод дополнительных переменных. Кратко изложим необходимые понятия и результаты из статьи [1].

Бабаджанянц Левон Константинович — доктор физико-математических наук, профессор кафедры механики управляемого движения факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 88. Научные направления: дифференциальные уравнения, классическая и небесная механика. E-mail: [email protected].

Брэгман Константин Михайлович — аспирант кафедры механики управляемого движения факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Научный руководитель: доктор физико-математических наук, проф. Л. К. Бабаджанянц. Количество опубликованных работ: 1. Научные направления: дифференциальные уравнения, вычислительная математика. E-mail: [email protected].

© Л. К. Бабаджанянц, К. М. Брэгман, 2012

1.1. Дифференциальные уравнения. Полиномиальные системы,. Полагая, что

X = (хи...,хт) е Ят, Ь е Я, а = («1 ,...«) € / = /), Д е Я, рассмотрим уравнения

дхг/дЬ^ = Д(х,а), г е [1 : т], ] е [1 : в], дх/дЬ = /(х,а), (1)

которые в случае в = 1 можно записать в виде системы ОДУ.

Полиномиальной назовем систему (1), в которой все функции Д - полиномы по х 1, ..., х^т с коэффициентами, зависящими от параметра а. Будем говорить, что скалярная функция р аргумента х = (х1,...,ха) удовлетворяет полиномиальной системе, если она является одной из компонент вектор-функции (того же аргумента х = (х1,...,ха)) - решения некоторой полиномиальной системы. Класс скалярных функций аргумента х = (х1, ...,ха), удовлетворяющих полиномиальной системе, обозначим £а. Так как любую функцию аргумента х = (х1,...,ха1), принадлежащую классу £а1, можно считать и функцией аргумента х = (х1 ,...,ха2) класса £а2 при <71 < а2, то примем, что а1 < а2 ^ £а1 С £а2, т. е. £ С £2 С .... Большое число специальных функций, представленных в справочниках и компьютерных системах, принадлежит £1 (и тем более £а при а > 1).

1.2. Метод дополнительных переменных. Этот метод сводит уравнения (1) к полиномиальной системе. Он состоит в том, что находят дополнительные переменные, функции хт+1(х1,..., хт),..., хт+к (х1,..., хт), которые удовлетворяют следующим условиям:

а) правые части уравнений (1) - полиномы по х1,..., хт, хт+1,..., хт+и;

б) дхт+\/дхг - полиномы по х1,...,хт+^, I = 1,...,к, г = 1,...,т.

1.3. Классы Тт. Классом Т^,а, т е [1 : называют множество скалярных функций аргумента х = (х1,... ,хт), которые можно получить из х1,... ,хт при помощи конечного числа операций +, —, х,/ и конечного числа функций р1, . ..,р1 е £а и их суперпозиций. Если в число этих функций включить р такую, что р(а) = 1/а, то в определении Т^ можно не использовать операцию /, так как величина а/Ь равна а х р(Ь).

Как показано в [1], для применимости к системе (1) МДП необходимо и достаточно, чтобы все / принадлежали какому-нибудь из классов . Предлагаемый в п. 2 для таких систем алгоритм МДП основан на понятии библиотеки.

1.4. Расширения и библиотеки. Если р е £а, то существует вектор-функция (р1,..., рп) - решение полиномиальной системы, где р1 = р. Множество {р1,..., рп} называем расширением р. Объединение расширений нескольких функций - библиотека. Фактически, библиотека содержит набор полиномиальных систем, а не задач Коши. Всякое подмножество библиотеки называем подбиблиотекой, если оно само - библиотека. Разделом называем подбиблиотеку, которая не пересекается ни с одной другой подбиблиотекой, не являющейся ее частью. Объединение разделов является разделом. Раздел будем считать простым, если он не содержит других разделов. Библиотека объединяет свои разделы, и можно сказать, что она делится на них и, в частности, сама может быть своим единственным (простым) разделом.

2. Алгоритм сведения к полиномиальной системе. Здесь последовательно рассмотрим структуру библиотеки как таблицы, алгоритм сведения и пример.

2.1. Библиотека как таблица. Библиотеку будем называть автономной или неавтономной соответственно тому, какими уравнениями она представлена - автономными или неавтономными. Ради удобства и большей универсальности предлагаемый

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

Для пользователя библиотека - пополняемая таблица, состоящая из потенциально неограниченного количества строк и одиннадцати столбцов:

«IDFun» - идентификатор функции (совпадает с ее именем в пакете «Mathematica», если такая функция в нем определена);

«FName» - второе имя функции;

«CNo» - текущий (current) порядковый номер функции в библиотеке;

«ONo» - исходный (original) порядковый номер функции в главной библиотеке;

«ENo» - порядковый номер функции в ее расширении (extention); так как расширение функции не единственно, то количество функций расширения и его состав зависят от составителя библиотеки;

«FArgN» - количество аргументов функции; у различных функций, входящих в расширение, может быть разное количество аргументов, однако можно (и мы будем) считать, что аргументы (а значит, и их количество) у функций расширения одинаковы, приняв, что производные функций по недостающим аргументам равны нулю;

«jENo» - порядковый номер аргумента, производной по которому соответствует данная строка: если, например, в текущей строке ENo=5, jENo=3, то она содержит информацию о производной пятой функции по третьему аргументу;

«FENoList» - список номеров «ONo» функций расширения данной функции;

«RHSPoly» - полином, правая часть дифференциального уравнения для данной функции у: число правых частей (а значит, и строк с одним ENo и разными jENo) равно количеству FArgN аргументов pi,... этой функции (переменные полинома - функции у = у1 ,у2,... из ее списка расширения и их аргументы pi,...);

«LPartNo» - номер раздела библиотеки, содержащего функцию;

«IVFun» - имя программы, вычисляющей значение функции (такой программы может и не быть, в алгоритме этот столбец не используется и носит чисто информационный характер для пользователей).

2.1.1. Зависимость коэффициентов полиномов от параметр о в. Коэффициенты полиномов (правых частей дифференциальных уравнений) в графе «RHSPoly» могут зависеть от параметров. Относительно них мы предполагаем выполненным единственное условие: во всех уравнениях для функций одного раздела библиотеки параметры, обозначенные одинаково, должны иметь один и тот же смысл. При этом сохраняется возможность использовать одни и те же символы для разных или одинаковых по смыслу параметров в различных разделах. Из сказанного следует, что библиотеку целесообразно представлять разбитой на простые разделы.

2.1.2. Разбиение библиотеки на разделы. Для такого разбиения можно предложить различные алгоритмы, основанные на очевидном утверждении: если в списке расширения функции есть функция из некоторого раздела, то и сама функция принадлежит тому же разделу. Кроме того, пользуясь данным утверждением, при пополнении библиотеки новыми функциями имевшуюся нумерацию разделов можно сохранить, а новым разделам (если появятся) присвоить иные номера (старые разделы могут пополниться новыми функциями).

2.1.3. Зависимость функций от аргументов и параметров. Имя в графе «FName» задается в формате Name(^b...;p1,...), где p1,... - список аргументов функции, а fi1,... - список параметров.

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

а) находим все различные номера О№ из всех списков в графе «^Е^ЫэЬ» для функций ф1,...,фк,

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

2.2. Алгоритм сведения. Вначале опишем элементарное преобразование уравнений вида (1), составляющее основу алгоритма, а затем схему алгоритма.

2.2.1. Элементарное преобразование системы. Рассмотрим систему (1) при условии, что все принадлежат Тт. Алгоритм преобразования следующий:

а) ищем в правых частях функцию вида ук (Р1 (х\,..., хт), . ..,Р, (х\,..., хт)), где Ри - алгебраические полиномы, у = у(р1,... ,р,) - функция из библиотеки, а -натуральные числа;

б) пусть у = у1,...,ул - все функции расширения функции у и пусть

д<Уг/др„ = (У1 ...у^^ г € [1 : ^ где ОЧ - полиномы. Вводим новые дополнительные переменные

Хт+1 = У1 (Р1(Х1, ..., Хт), . . .,Рт)(х-1, . .., Хт)), хт+л у л (Р1 (х1 ■} ..., хт), . .., Р, (х1 ■} ... хт));

в) заменяя все функции у1(Р1,. ..,Р,),..., ул(Р1,... ,Р,) в правых частях уравнений (1) на новые переменные хт+1,..., хт+л соответственно, получаем

дх^/дЬ^ = д%(х1,...,хт,хт+1 ,...,хт+л), г € [1 : т], 3 € [1 : в];

г) выписываем уравнения для введенных переменных (г € [1 : -], 3 € [1 : в]):

, т

дх"т+1/дЬз ^ ^ (хт+1 ^ ... хт+л) ^ ^ дРь> /дхк ' дк (х1... хт ^ хт+1 ^ ... хт+л) .

к=1

2.2.2. Схема сведения.

Шаг 1. Заносим правые части (1) и начальные данные в пакете «МаШетайса»:

дх^/дЬз = Ц (хь ..., хп), г € [1 : п], 3 € [1 : в], (2)

х1(Ьо) = х1,о,..., хп(Ь0) = хп,о.

Введенные правые части и начальные данные могут содержать параметры. Если для обозначения каких-то из этих параметров использованы символы вида а^, где г -натуральное число, то полагаем, что П равно максимальному из таких г, а если таких параметров нет, то считаем, что П = 0.

Образуем подбиблиотеку, состоящую из объединения расширений всех функций, участвующих в написании правых частей системы (2) (см. п. 2.1.4). Параметры, от которых зависят функции подбиблиотеки, переобозначим символами 1, 2,... и заполним таблицу соответствий вида

NewParam OldParam LPartNo

1 «2 2

"0+ 2 d 7

Шаг 2. Преобразуем (2) шаг за шагом, пока не сведем систему к полиномиальной форме, применяя на каждом шаге элементарное преобразование к системе, полученной на предыдущем шаге (см. п. 2.2.1). На последнем шаге для переменных xi,... ,xn+M выводим систему

dxi/dtj = Rj(Xi,...,Xn+M), i e [1 : n + M], j e [1 : s],

где Rj - полиномы.

Шаг 3. Вычисляем начальные значения для переменных xn+i,... ,xn+M.

Таким образом, исходная задача Коши сведена к полиномиальной.

Приведенная на рисунке блок-схема отражает описанный выше алгоритм сведения системы (1) к полиномиальной форме и соответствует тексту программы, электронный вариант которой можно получить по электронной почте: [email protected].

2.2.3. Пример сведения. Рассмотрим сведение к полиномиальной форме полной системы из шести уравнений в частных производных для трех функций от двух аргументов, причем правые части зависят от шести элементарных функций одного аргумента, одной специальной функции двух аргументов и трехпараметрического семейства функций Вебера. Как можно заметить, эти уравнения выбраны так, что при s = 1 образуют систему ОДУ. Данный модельный, но отнюдь не тривиальный пример не только иллюстрирует алгоритм сведения п. 2.2.2, но и предназначен для отладки соответствующей программы. Вначале опишем используемую библиотеку.

Библиотека. Она является объединением шести подбиблиотек Пб1, ... , Пб6, которые рассмотрим вместе с соответствующими дифференциальными уравнениями. Функции будем обозначать p,pi,..., а их аргументы - p,p1,... .

Пб1. Состоит из одной функции одного аргумента - p = inv(p) = p-i, которая удовлетворяет уравнению dp/dp = — p2, а расширение этой функции включает только ее одну.

Пб2. Входят две функции одного аргумента pi = lnp, р2 = inv(p), которые удовлетворяют системе dpi/dp = р2, dp2/dp = —р2, причем расширение функции pi состоит из функций p1, p2, а расширение функции p2 - из нее одной (т. е. совпадает с Пб1).

Пб3. Содержит две функции одного аргумента - pi = sinp, p2 = cosp, которые удовлетворяют системе dp i/dp = p2, dp2/dp = —pi, а расширение каждой из функций pi, p2 включает функции pi, p2.

Пб4. Состоит из двух функций одного аргумента - pi = sh(p), p2 = ch(p), которые удовлетворяют системе dpi/dp = p2, dp2/dp = pi, а расширение каждой из функций pi, p2 содержит функции pi, p2.

Пб5. Состоит из четырех функций двух аргументов: pi = EK(pi,p2) = E(e, M), где EK = E - эксцентрическая аномалия, pi = e - эксцентриситет, p2 = M — средняя аномалия [3]:

p2 = EKs(pi,p2) = sin pi, p3 = EKc(pi,p2) = cos pi,

Вход Сообщение

и заведение системы об ошибке

Нет

Является ли полиномиальной полученная система?

Нет

Да

Заведение и проверка данных

Подготовка

Найти в системе

функцию, аргументы

которой - полиномы

Преобразо-

вания

Вывод таблицы параметров, переменных, начальных данных и уравнений

Вывод данных

Блок-схема алгоритма и программы сведения

y4 = EKi(pi,p2) = (1 - Р1 cos yi) 1, которые удовлетворяют системе

dyi/dpi = y 4, dyi/dp2 = y 4,

dy2/dpi = У2У3 У 4, dy2/dp2 = У3У4,

дуз/др1 = -у2у4, дуз/dp2 = -У2У4,

dy4/dpi = узу2 -piy"2y4, dy4/dp2 = -piy\y2,

причем расширение функции yi состоит из yi,y2,у3,у4, а расширение каждой из функций у2,уз,у4 - из всех этих трех функций.

Пб6. Состоит из функций yi, у2 аргумента p, причем первая из них - функция Вебера Dv [4], а вторая - ее производная, которую обозначим DV. Как известно, они удовлетворяют системе уравнений (a, b, c - параметры)

dyi/dp = y2,

dy2/dp = - (ap2 + bp + c)yi.

Расширение каждой из функций yi,y2 есть {yi,y2}. Перейдем к самому примеру.

Шаг 1. Пользователь заносит систему и начальные условия (j = 1, 2):

dxi/dtj = x3(sincos(aln2 x2 + bx3) + bln4 x2) + sin5(aln2 x2 + bx3) +

+ Dv(g2,1, -1; xi), dx2/dtj = x2 cossin(aln2 x2 + bx3) + cos4(aln2 x2 + bx3) +

+EKj(sinsin(aln2 x2 + bx3),xi), dx3/dtj = xi(ch2(aln2 x2 + bx3) + sin(aln2 x2 + bx3)) + + sh5 (a ln2 x2 + bx3),

(3)

a = (ai, a2) = (a, b) - вещественные постоянные (параметры),

xi(to) = xi,o, x2 (to) = x2,0, x3(to ) = x3,0.

Подбиблиотека состоит из объединения Пб1,... ,Пб6, т. е. из функций

lnp, p-i, sinp, cosp, shp, chp;

EK(pi,p2), EKs(pi,p2), EKc(pi,p2), EKi(pi,p2), Dv(a,b,c; p).

Шаг 2. Преобразуем систему (3) шаг за шагом (шаг 2.1, шаг 2.2,...), пока не будем иметь ее в полиномиальной форме, применяя на каждом шаге элементарное преобразование к системе, полученной на предыдущем шаге.

Шаг 2.1:

а) в качестве первой функции y(P) возьмем, например, ln x2;

б) функциями расширения lnp будут функции yi = lnp,y2 = p-i, которые удовлетворяют системе

dyi/dp = y2, dy2/dp = -y2>.

Введем дополнительные переменные x4 = ln x2, x5 = x-1 ;

в) заменяя lnx2, x- 1 во всех их вхождениях в правые части уравнений (3) на Х4, соответственно, получаем новую запись этих исходных уравнений (j = 1, 2):

dx1/dtj = x|(sincos(ax2 + bx3) + bx\) + sin5 (ax\ + bx3) + Dv(g2,1, —1; x1),

dx2/dtj = x2 cossin(ax2 + bx3) + cos4(ax2 + bx3)++EKj(sinsin(ax2 + bx3),x1),

Bxs/Btj = x"f(ch2(ax4 + bx¡) + sin(ax2 + bx¡)) + sh5(ax2 + bx¡);

г) уравнения для введенных дополнительных переменных следующие:

dx4/dtj = д ln x2/dtj = x-1 • dx2/dtj = x5(x2 cossin(ax2 + bx3) + + cos4(ax2 + bx3) + EKj (sinsin(ax2 + bx3),x1)),

dx5/dtj = dx-1 /dtj = —x-2 • dx2/dtj = — x^^ cossin(ax2 + bx3) + + cos4(ax2 + bx3) + EKj(sinsin(ax2 + bx3), xjJ).

Итак, в результате шага 2.1 будем иметь систему (j = 1, 2)

dx1/dtj = x|(sincos(ax4 + bx3) + bx4) + sin5 (ax^ + bx3) + Dv(g2,1, —1; x1),

dx2/dtj = x¡ cossin(ax2 + bx3) + cos4(ax2 + bx3) + +EKj (sinsin(ax2 + bx3 ),x1),

dx^/dtj = x"f(ch2(ax4 + bx¡) + sin(ax2 + bx¡)) + sh5(ax4 + bx¡),

dx4/dtj = d ln x2/dtj = x-1 • dx2/dtj = x5(x¡ cossin(ax2 + bx3) + + cos4(ax2 + bx3) + EKj (sinsin(ax2 + bx3),x1)),

dx5/dtj = dx-1 /dtj = —x-2 • dx2/dtj = — x^^ cossin(ax2 + bx3) + + cos4(ax2 + bx3) + EKj (sinsin(ax2 + bx3),x1)).

Шаг 2.2 - шаг 2.7. Действуя аналогично, вводим переменные

x6 = sin(ax2 + bxs), x7 = cos(ax4 + bx¿),

xg = sh(ax4 + bxs), xg = ch(ax4 + bx%), x10 = cos x6, xn = sin x6, x12 = sin x7, x13 = cos x7, x14 = EK(xn,x1), x15 = EKs(xn ,x1), x16 = EKc(xn,x1), x17 = EKi(xn,x1), x18 = Dv(g2,1, —1; x1), x1g = DV (g2, 1, —1; x1), и в результате шага 2.7 получаем итоговую систему (j = 1, 2)

dx1/dtj = x|(x12 + bx4) + x6 + xJ18, dx2/dtj = x2 x10 + xíJ + x"[4, dxs/dtj = x{(xl + x6) + xg, dx4/dtj = x5(x^ x10 + xíJ + x^), dx5/dtj = —x5(x2 x10 + x7 + xj14),

дхб/dtj = xr(2ax4x5 (x^xio + x4 + ) + b(xi(x2 + x^) + xf)),

dxi/dtj = -xe,(2ax4x5(x2 xio + x7 + x"4) + 6(xl(xg + x6) + x|)),

dxg/dtj = xg(2ax4 x5(x2 xio + x4 + x"4) + b(x\(x1 + x6) + x|)),

dxg/dtj = xg(2ax4 x5(x2 xio + x4 + x"^) + b(x\(x1 + xe) + x|)),

dxio/dtj = —xnxr(2ax4 x5(x^xio + xTj + x"4) + b(xJ1(x'9 + x6) + x5)),

dxio/dtj = —xnxr(2ax4 x5(x^xio + xTj + x"4) + b(x\(x\ + x6) + xf)),

dxii/dtj = xio x7(2ax4x5 (x\x\o + x7 + x"4) + b(x{(x9 + x6) + xf)),

dxyi/dtj = —xi3xe(2ax4 x5(x2>xw + x| + x"4) + b(x{(xl + xe) + xf)),

dxis/dtj = xi2 xe(2ax4x5 (x^xio + x7 + x"4) + b(x'i(x2g + x6) + xf)),

dxi4/dtj = xi5 xi^xioxr (2ax4x5(x'¡2 xio + x4 + x"4) + b(x\(x2 + xe) + xf)) + + xi7(xl (xi2 + bx;|) + x'6 + j),

dxi5/dtj = x 15x16x17x1o x7(2ax4x5 (x2xio + x4 + x"4) + b(xJi(xg) + x6) + xf)) + + xi6 xi7 (x3(xi2 + bx4) + xl + x"f),

dxi6/dtj = —x^x^ xio x7(2ax4x5(x2 xio + x4 + x"4) + b(x" (x2 + x6) + xf)) — — xi5 xi7 (x3(xi2 + bx4) + x5 + x\f),

dxi7/dtj = (xi6 x27 — xiixÍ5x37)xiox7 (2ax4x5(x2 xio + x7 + x"4) +

+ b( x\ (xg + x6) + x^f )) — xiix25x37 (xj| (xi2 + bx4) + x5 + x{f),

dxif /dtj = xig (xl(xi2 + bx4) + x5 + xjif) — x2f xig(x2, xw + x7 + x""4 ) — — x\f xi9 (x"1 (x9 + x6) + xf),

dxig/dtj = — xif(g2 xi + xi — 1)(x|(xi2 + bx4) + x5 + j ) — — xifxig(x2 xio + xi 7 + x^) — x^ xig (xi (x2 + x6 )+ xf ).

Исходные и дополнительные переменные и начальные данные для них следующие:

xi, x2, x3; xi,o, x2,o, x3¡o; x4 = ln x2, x4,o = ln x2,o; x5 = x2 , x5,o = x2,o; x6 = sin(ax2 + bx3), x6,o = sin(ax2 o + bx3,o); x7 = cos(ax2 + bx3), x7,o = cos(ax2 o + bx3,o); xf = sh(ax4 + bx3), xf ,o — sh( ax 4 o + bx3¡o); xg = ch(ax^ + bx3), xg,o = ch(ax4 o + bx3,o); xio = cos x6, xioto = cos x6to; xii = sin x6, xiio = sin x6to; xi2 = sin x7, xi2,o = sin x7,o;

xi3 = cos X7, xi3,o = cos X7,o; xu = EK(xn,xi), X14,0 = EK(xii,o, xi,o); xi5 = EKs(xii,xi), xi5,o = EKs(xii,o ,xi,o); xie = EKc(xii,xi), xi6,o = EKc(xii,o, xi,o); xir = EKi(xii,xi), xi7,o = EKi(xii,o, xi,o); xis = Dv(g2,1, -1; xi), xi8,o = Dv(g2,1, -1; xi,o); xiq = DV(g2,1, -1; xi), xig,o = DV(g2,1, -1; xi,o).

Литература

1. Бабаджанянц Л. К. Метод дополнительных переменных // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2010. Вып. 1. С. 3—11 .

2. Wolfram Mathematica Documentation Center // URL: http://reference.wolfram.com/mathematica/ guide/Mathematica.html.

3. Холшевников К. В., Титов В. Б. Задача двух тел. СПб.: Изд-во С.-Петерб. ун-та, 2007. 180 с.

4. Абрамовиц М., Стиган И. Справочник по элементарным функциям / пер. с англ.; под ред. В. А. Диткина, Л. Н. Кармазиной. М.: Наука, 1968. 800 с. (Abramowitz M., Stegun I. A. Handbook of mathematical functions.)

Статья рекомендована к печати проф. Л. А. Петросяном. Статья принята к печати 28 февраля 2012 г.