Численная стабилизация уравнений движения небесных тел

Министерство образования и науки Российской федерации

Томский государственный университет

Физический факультет

Реферат

Численная стабилизация

уравнений движения

небесных тел

Томск-2004

Содержание

1. Введение 3

2. Численная стабилизация с применением всех законов сохранения в задаче многих тел 4

1. Асимптотическая стабилизация с использованием интеграла

движения 4

1. Асимптотическая стабилизация с помощью энергетического

соотношения 5

2.3 Асимптотическая стабилизация с помощью законов сохранения

в задаче N-тел 6

3. Численная стабилизация дифференциальных уравнений

кеплеровского движения 7

3.1 Неустойчивость классических уравнений 7

3.2 Стабилизация кругового кеплеровского движения 8

3.3 Численный эксперимент 9

4. Метод Накози 10

4.1 Вводные замечания 10

4.2 Уравнения связей 12

4.3 Оценка метода 14

5. Об эффективности и точности метода Накози при интегрировании

ограниченной задачи трех тел 15

5.1 Вводные замечания 15

5.2 Оценка метода при использовании различных интеграторов 16

5.2.1 Поправки по многообразию 16

5.2.2 Тестовые орбиты 16

5.3 Выводы 17

6. Литература 18

1. Введение

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

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

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

В настоящем реферате мы рассматриваем работы, в которых даны два разных подхода к проблеме стабилизации решений динамических уравнений, это – различные способы асимптотической стабилизации уравнений движения, предложенные Дж. Баумгартом [1,2] , и метод П. Е. Накози [3], основанный на удержании решения около некоторой интегральной поверхности. Кроме того, мы рассмотрим кратко работу М.А. Марисона [4], в которой приводятся результаты анализа точности и эффективности метода Накози при использовании его с различными интеграторами.

2. Численная стабилизация с применением всех законов сохранения

в задаче многих тел.

2.1 . Асимптотическая стабилизация с использованием интеграла движения

В своей статье [1] Баумгарт рассматривает механическую систему, составленную из N точек, где уравнения движения относительно инерциальной декартовой системы координат и времени t задаются в следующем виде

, (2.1)

где F_i внешние силы. Предполагается, что соотношение в форме

(2.2)

можно получить из этой системы. Поэтому система имеет первый интеграл

(2.3)

где постоянная k определяется начальными условиями. Уравнение (2.3) принимаются за внутреннюю неголономную связь

В течение численного интегрирования уравнений движения (2.1) значения функции , в котором k должна быть константой, используются как проверка. В силу неустранимой, вычислительной ошибки, которую компьютер достигает после n шагов на момент времени t₀ , имеет значение

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

которое является свойством наличной дифференциальной системы

Автор желает улучшить эту ситуацию так, чтобы значение уменьшилось до k в течение дальнейшего интегрирования. Чтобы выполнить это он вводит связь f, которая уменьшит ошибку элемента k. Таким образом, связь представляется в форме

(2.4а)

где f линейна, относительно ускорения.

Если система подвергается связям, которые содержат ускорения, то можно использовать принцип Гаусса.

Гауссова связь выбирается в форме дифференциального уравнения

(2.5)

где положительная функция.

Гауссова связь обладает свойством понижения ошибки элемента, которая в данном случае должна иметь значение k. Начальная связь предполагает выбор , который приводит к N=const, и не уменьшает ошибку элемента k. Поэтому необходим выбор . Выбор подходит для многих задач, при этом неверное значение уменьшается экспоненциально.

Правило Гаусса дает уравнения движения, которые удовлетворяют связи [1]

(2.6)

где множитель Лагранжа.

В случае одной связи очень легко решить систему (2.6) для ускорений и исключить множитель Лагранжа. Результатом является следующая дифференциальная система, подходящая для численного интегрирования:

(2.7)

Правая сторона уравнения (2.7) есть фиктивная сила реакции или стабилизирующая сила. Эта сила равняется нулю в аналитическом вычислении, но не для компьютера.

2.2 Асимптотическая стабилизация с помощью энергетического соотношения

На примере стабилизирующего метода (2.6) в качестве (2.5) Баумгарт рассматривает энергетическое соотношение. В этом случае внешние силы зависят от потенциала U, который зависит только от вектора положения x_i:

(2.8)

Это приводит к хорошо известному соотношению энергии

(2.9а)

где E постоянная полная энергия.

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

(2.9b)

Чтобы применить свой метод для уменьшения ошибок Баумгарт вводит Гауссову связь

(2.10)

и получает энергостабилизирующие уравнения движения:

(2.11а)

или после сокращения:

(2.11b)

Этот метод уменьшения ошибок можно легко расширить для возмущенного случая. Получается

(2.12а)

вместе с

(2.12b)

где малый возмущенный параметр и возмущенные силы.

2.3 Асимптотическая стабилизация с помощью законов сохранения в задаче N-тел

Обозначим вектор положения точки с массой m_i как r_i. Уравнение движения есть

(2.13)

где U хорошо известный гравитационный потенциал.

Автор предполагает, что центр масс неподвижен и

(2.14)

Это соотношение легко выполняется путем использования соответствующей системы координат.

Для уравнения движения (2.14) существует следующие известные законы сохранения:

Угловой момент

(2.15)

Соотношение энергии

(2.16)

где c постоянный вектор углового момента и E постоянная общая энергия.

Стабилизирующие Гауссовы связи есть:

(2.17)

и

(2.18)

Выражения (2.17) и (2.18) линейны относительно ускорения.

Уравнения движения, которое содержит фиктивные стабилизирующие силы реакции

связи имеют теперь форму

(2.19)

где векторный множитель, а скалярный множитель Лагранжа.

Из уравнения (2.19) видно, что суммированием по всем i можно избавиться не только от внутренних сил, но также от фиктивных стабилизирующих сил

(2.20)

которые сочетаются с уравнением (2.14):

Однако фиктивные силы реакции не разрушают соотношение (2.14), поскольку теперь вводятся подогнанные координаты, которые удовлетворяют соотношению (2.14) автоматически.

Но перед этим должна быть разрешена система уравнений движения (2.17)-(2.19) относительно ускорения r , и исключены 4 скалярных множителя Лагранжа.

Затем автор вводит исправленные координаты. Первая возможность есть:

(2.21)

Другая возможность есть введение координат Якоби..

С помощью введения исправленных координат понижается степень свободы системы.

Некоторые замечания:

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

- если динамическая задача N-тел является слабо возмущенной, то вектор углового момента c и общая энергия E незначительно отличаются от постоянных величин и метод стабилизации еще может быть применен при этом дополнительные уравнения

(2.22)

должны быть приняты во внимание.

3 Численная стабилизация дифференциальных уравнений кеплеровского движения

3.1 Неустойчивость классических уравнений

В полярных координатах уравнения движения есть

(3.1)

K² гравитационный параметр притягивающего центра масс и точка обозначает дифференцирование по времени t. Уравнение энергии

(3.2)

есть первый интеграл уравнения (3.1) и h отрицательная общая энергия частицы. Частное решение уравнения (3.1) есть круговые движения радиусов . Их полярные углы соответственно

(3.3)

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

Теорема 1: Дифференциальная система (3.1) является неустойчивой в смысле Ляпунова, так как изменение (3.3) со временем в кругом движении зависит от радиусов.

В статье Баумгарта [2] теорема 1 подтверждается обсуждением вариационных уравнений системы (3.1). За опорное круговое движение принимается r=a с

(3.4)

Результирующее вариационное уравнение имеет вид

(3.5)

(3.6)

А частное решение

, (3.7)

соответствует ограниченному круговому решению. Общее решение системы (3.5), (3.6) получается путем добавления чисто периодических членов, которые не существенны для исследования стабилизации. Из правой части соотношения (3.7) следует теорема 2.

Теорема 2: В случае кругового движения нестабильность теоремы 1 линейна; угловой коэффициент есть .

Далее обсуждаются особенности показателей степени вариационных уравнений.

Результат записывается в виде матрицы:

Собственные значения матрицы есть . Ранг матрицы равен 3, следовательно существует только один собственный вектор соответствующий повторяющемся собственному значению, равному нулю. Эта ситуация обеспечивает линейную неустойчивость. Задача устойчивости означает уменьшение ранга до значения 2.

3.2 Стабилизация кругового кеплеровского движения

Много попыток было сделано, чтобы довести до конца эту задачу. Они были неудачны, так как физическое время t использовали как независимую переменную. Баумгарт также продолжил изучать эту задачу, но за новую независимую переменную взял s, которая удовлетворяет следующим требованиям. Соотношение определяется с помощью устойчивого дифференциального уравнения, которое не связано с уравнением орбиты. Это достигается преобразованием уравнения (3.1) к виду

(3.8)

Использовалось уравнение энергии (3.2). Новая независимая переменная s определяется дифференциальным соотношением.

(3.9)

Уравнение (3.8) становится

(3.10)

или

(3.11)

Это дифференциальное уравнение для вычисления физического времени t устойчиво в смысле Ляпунова при условии h>0.

Замечание: Общее решение уравнения (3.11) есть

(3.12)

В возмущенном случае может быть использовано любое дифференциальное уравнение подобное уравнению (3.11) или формула (3.12) где не постоянные величины, но медленно меняющиеся функции. Причем первый член в (3.12) есть временной элемент в смысле Штифеля и Шайфеля.

Основные уравнения (3.1) и (3.2) перепишутся после введения (3.9.) как:

(3.1а)

(3.2а)

(3.13)

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

Теорема 3: Дифференциальная система (3.1а) неустойчива в смысле Ляпунова.

Систему (3.1а) дополняется регулирующим членом, который используется в левой части уравнения (3.2а) и он фактически равен нулю. Таким образом, r – уравнение (3.1а) перепишется как

(3.14)

или

(3.15)

– весовой множитель. Правая часть – уравнения (3.1а) не меняется. Контрольный член эквивалентен фиктивной центральной силе.

Независимость угловой скорости увеличивает надежду, что дифференциальная система

(3.16)

(3.17)

более устойчива и лучше ведет себя, чем система (3.1).

Теорема 4: В случае кругового движения неустойчивость системы из теоремы 3 линейна; угловой коэффициент .

Введение фиктивно времени s . ослабляет нестабильность посредством множителя (Теорема 2).

Баумгарт рассматривает стабилизирующую систему (3.16) при ().

(3.18)

Общее решение стабилизирующей системы (3.18) есть

(3.19)

Следовательно, вариационные уравнения устойчивы. При введении

(3.20)

уравнение (3.18) преобразуется как

(3.21)

Собственные значения матрицы есть 0, 0, , ее ранг равен двум, и. существуют два линейно независимых собственных вектора. Эта ситуация дает стабильность вариационных уравнений.

1. Численный эксперимент

В рассматриваемой работе Дж. Баумгарта были сравнены 3 метода для вычисления кеплеровской орбиты. Так как орбиты лежат в плоскости, точка на орбите определятся двумя координатами, которые мы запишем в векторной форме . Были рассмотрены следующие дифференциальные уравнения

(1): Классические уравнения задачи двух тел

(2): Стабилизированные уравнения

(3): Уравнение осциллятора

Эти дифференциальные уравнения интегрировались численно методом Рунге-Кутта 4 порядка, используя фиксированный размер шага независимой переменной, t для (1) и s для (2) и (3). За 1 оборот было выполнено 100 шагов интегрирования.

За все вычисления большая полуось эллиптической орбиты имела постоянную длину, равную одной единице, пока эксцентриситет менялся.

Ошибка в расстоянии вследствие численного интегрирования была определена как

(3.22)

где

(3.23)

В (3.2) и (3.3) это определение было интерпретировано следующим образом. Пусть s задано в настоящий момент, для которого оценена ошибка и пусть t соответствующее значение времени t вычисленного посредством численного интегрирования. Точные координаты были определены для этого значения времени t путем решения кеплеровского уравнения. Дальнейшее вычисление координат получено путем приведения численного интегрирования к аргументу s.

• Было проведено два краткосрочных эксперимента, в которых использовались небольшие эксцентриситеты. Вычислены два оборота тела с началом в перицентре.

Как видно из результатов самый точный метод под номером (3), но метод под номером (2) не намного хуже. Метод под номером (1) имеет очень плохую точность и ухудшается для больших эксцентриситетов.

• Так же был проведен один краткосрочный эксперимент, в котором использовались большие эксцентриситеты. Вычислялся один оборот тела, с началом в апоцентре. Метод под номером (3) очень точен благодаря регуляризации. Метод (2) ухудшается для слишком больших эксцентриситетов из-за сингулярности в начале системы координат.

• Был проведен один долгосрочный эксперимент, с использованием больших эксцентриситетов. Ошибка в расстоянии после 50 оборотов показана в следующей таблице:

0.9 0.95 0.99	5.3*10-4 7.6*10-4 5.3*10-3	8.2*10-5 1.2*10-4 2.7*10-4

Метод под номером (3) работает даже при e=0,95.

4. Метод Накози

4.1 Вводные замечания

П. Накози [3] представил метод, который эффективно использует интегралы в численном интегрировании гравитационной системы. Метод приводит к решениям высокой точности, в то время как он использует меньше времени вычисления, чем обычные процедуры численного интегрирования, которые используют интегралы непосредственно.

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

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

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

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

Ответ, кажется, зависит от способа использования интегралов. Интегралы могут использоваться для понижения порядка системы до (6n-p), но интегралы движения пониженной системы часто не сохраняются во время вычисления.

Результирующие уравнения движения пониженной системы могут быть на много сложнее чем первоначальная система порядка 6n (например, появляется нелинейность, введенная интегралом энергии). Кроме того, уравнения порядка (6n-p) могут потерять симметричность первоначальной системы.

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

Идея использования поправок, найденных методом наименьших квадратов, чтобы сохранить интегралы, имеет геометрическую интерпретацию. В течение интегрирования, ошибки в вычислении могут послужить причиной тому, что решение покинет (6n-p) – мерную гиперповерхность, определенную интегралами. Поправки наименьших квадратов к 6n переменным возвращают решение на поверхность вдоль нормальной поверхности. С помощью исправления переменных, решение остается на первоначальной гиперповерхности в течение численного интегрирования.

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

4.2 Уравнения связей

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

(4.1)

два уравнения

(4.2)

должны быть решены с помощью уравнения (4.1), чтобы определить величины x, y, и . Где, множитель Лагранжа.

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

Пусть (4.3)

интеграл системы. Уравнение (4.3) определяет 3 – мерную гиперповерхность, помещенную в фазовое 4 – мерное пространство.

В течение процесса численного интегрирования системы, вычисленное решение, полученное на время t, выглядит следующим образом:

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

(4.4)

где малая величина. Решение покидает интегральную поверхность, определенную уравнением (4.3) и остается на поверхности определенное уравнением (4.4). Это используется для того, чтобы получить поправки и вычислить вектора , , так чтобы

(4.5)

Квадрат величины вектора поправок можно записать как

(4.6)

Поправки выбираются так чтобы функция уравнения (4.6) минимизировалась, при условии связей (4.5). Для 4 – мерного случая уравнения (4.2) имеют вид

(4.7)

Уравнения (4.5) и (4.7) решаются для пяти неизвестных и . Уравнение (4.5) линеаризуется обычным способом

(4.8)

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

Решение уравнений (4.7) и (4.8) для поправок при условии (4.3) и (4.4) имеет вид

(4.9)

Вектор поправок добавляется к вычисленному вектору состояния , чтобы получить новый вектор состояния который удовлетворяет интегралу (4.3), с ошибкой порядка . Уравнение плоскости, заданное (4.8) исключает члены второго и выше порядков.

Результат (4.9) можно обобщить на динамическую систему порядка 6n , имея p интегралов. Вектор состояния системы , где вектор столбец в фазовом пространстве с компонентами . Обозначим вектор положения как , а вектор скорости как . Вектор состояния можно записать в следующем виде

(4.10)

Следовательно, уравнения движения системы есть

(4.11)

где вектор F есть функция вектора x и времени, p интегралов системы записываются как

(4.12)

Уравнение (4.11) решается численным интегрированием. Частные производные интегралов энергии (4.12) относительно компонент вектора состояния есть элементы матрицы . Так что,

На момент времени t, из-за ошибок вычисления, некоторые или все p компонент вектора E не равны нулю. Так что

где это вектор ошибок, чьи элементы малые величины.

Желательно вычислить вектор поправок так чтобы вектор

удовлетворял уравнению

Вектор выбирается так, чтобы величина была минимальной. Здесь, W весовая матрица и верхний индекс T означает операцию транспонирования матрицы.

Как в уравнении (4.8), каждый элемент вектора E разлагается по степеням вектора .

.

Члены второго порядка и выше отбрасываются, при , разложение уменьшается до:

(4.14)

Решение, заданное уравнением (4.7) принимает вид

(4.15)

Уравнения (4.14) и (4.15) решаются относительно (6n+p) неизвестных: компоненты двух векторов . Уравнение (4.15) разрешается относительно , результат подставляется в уравнение (4.14) и получаем

Последнее уравнение разрешается относительно и результат подставляется в уравнение (4.15). Решение для вектора поправок будет иметь вид

(4.16)

Для гравитационной системы, вектор F уравнения (11) задается как

(4.17)

где обозначает вектор столбец, чьи 3n компоненты есть .

U отрицательная функция потенциальной энергии системы и определяется как

Численное интегрирование системы уравнений (4.11) относительно F определенное уравнением (4.17), дает вектор решения на время t. Поправки могут быть вычислены с помощью уравнения (4.16).

4.3 Оценка метода

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

Метод был применен к гармоническому осциллятору и к системе двух тел с помощью следующей процедуры. Два множества решений получены путем численного интегрирования с различными начальными условиями. Первое множество решений получено без использования интегралов, в то время как другое множество получено с использованием поправок определенных рассматриваемым методом. Поправки применялись к вектору состояния на каждом шаге интегрирования. Все результаты интегрирования сравнивались с истинным решением системы, чтобы определить относительную точность неисправленного и исправленного решений. Решения гармонического осциллятора получены путем использования метода Рунге-Кутты 4 порядка с постоянным шагом. Оба неисправленное и исправленное решения гармонического осциллятора использовали такой же размер шага и такое же число шагов интегрирования. Решения системы двух тел были получены путем использования процедуры предиктора-коректора с переменным шагом. Оба решения, неисправленное и исправленное, для системы двух тел получены интегрированием с одним и тем же шагом и одним и тем же числом шагов интегрирования.

Применение метода к гармоническому осциллятору в фазовом 2 – мерном пространстве показывают не значительную разность в точности между исправленным и неисправленным решениями.

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

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

Метод был применен к гравитационной системе 25-тел, используя стандартные начальные условия. Сначала было выполнено высокоточное, неисправленное численное интегрирование системы. Ошибка усечения интегрирования была понижена до предельной компьютерной мощности. Систему интегрировали вперед и назад по времени. Это решение было получено как точное, стандартное решение с которым другие, менее точные сравнивались.

Программа численного интегрирования, которую использовали для оценки исправленного метода Рунге-Кутты-Филберга седьмого порядка переменного шага, была применена к задаче 25 тел. Два множества решений получены с помощью численного интегрирования. Первая система порядка 6n интегрировалась без использования интегралов. Затем систему порядка 6n интегрировали и все или различные комбинации 10 интегралов системы были использованы. Все решения сравнивались с наиболее точным стандартным решением. Кроме того, была выполнена проверка по результатам прямого и обратного интегрирования.

Как показывают результаты, метод Накози дает более эффективный численный процесс интегрирования, так как получается наибольшая точность при использовании одинакового времени вычисления и такая же точность при меньшем времени вычисления. Использование только интеграла энергии дает более эффективное численное интегрирование, чем интегрирование, использующее все 10 интегралов. Причина этого в следующем: (1) Ошибка в энергии существенно больше, чем ошибки в интегралах углового момента и центра масс. (2) Использование всех 10 интегралов требует обращения матрицы большой размерности, что неудобно и приводит к большим затратам времени

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

5. Об эффективности и точности метода Накози при интегрировании ограниченной задачи трех тел

5.1 Водные замечания

В работе [4] автор предпринял попытку применить метод стабилизации П. Накози к ограниченной задаче трех тел в сочетании с численным методом интегрирования Булирша-Штера.

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

(5.1)

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

В круговом случае , существует интеграл движения:

(5.2)

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

В этой статье рассматривается численные методы и необходимые условия. Первое, что рассматривается, это частный метод которым интегрируют уравнения, показано, что он является очень важным. Затем, пользуются константой Якоби, чтобы откорректировать численные ошибки усечения и округления в круговой ограниченной задаче. Этот метод можно применить к любой задаче, в которой существуют константы движения. Эффекты различных методов и техники сравниваются, используя орбиты в круговой и эллиптической ограниченных задачах трех тел..

5.2 Оценка метода при использовании различных интеграторов

5.2.1 Поправки по многообразию

П. Накози (1971) открыл способ использовать константы движения для проверки и поправки ошибок в положении и скорости. Когда существует интеграл уравнения движения, орбиту ограничивает поверхность в фазовом пространстве. Когда численная ошибка происходит, движение прыгает по различным поверхностям, соответствующим разным орбитам. Коррекционная схема П. Накози использует множители Лагранжа, чтобы найти методом наименьших квадратов наикротчайший путь назад на исходное многообразие. Мюресон (1988) применил эту поправку многообразия к обоим нерегулязированным и регулязированным уравнениям ограниченной задачи трех тел.

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

Если отклонения очень малы, тогда вычисленная орбита будет близка истиной орбите («истинная орбита» равна «вычисленной орбите с бесконечной точностью»). Следовательно, если отклонения систематические, тогда вычисленная орбита будет отклоняться от истиной орбиты и получение самого короткого обратного пути на многообразии будет давать систематическую компоненту в ошибках и вычисленная орбита будет значительно отклоняться от истинной.

В круговой ограниченной задачи трех поправки делались всякий раз, когда различия между текущей и начальной постоянной Якоби были не нулевые. Практически, это означает, что C была сохранена постоянной с точностью приблизительно 10^-16.

5.2.2 Тестовые орбиты

Чтобы показать полезность этих методов их проинтегрировали разными способами. Эти орбиты трудны для интегрирования, так как они близки к . Интегрирование продолжалось в течение 40 орбитальных периодов планеты.

В течение интегрирования, ошибка в постоянной Якоби контролировалась. Измерение накопленной ошибки также полезно в характеризования точности на орбите. Т.о. величина C определяется как

.

Экспериментально было показано, что стабилизация уменьшает величину C более чем на 5 порядков.

1. Выводы

Приведенные в статье результаты показывают, что экстраполяционный метод Булирши – Штера является очень быстрым и довольно точным. При ошибках такого же порядка как и ошибки, получаемые популярным методом Рунге-Кутты, интегрирование происходит на порядок быстрее. Многообразная коррекционная схема Накози (1971) может уменьшить ошибки фактически до нуля, и с удивительно малыми затратами машинного времени, зависящими от орбиты. К сожалению, реализация возможна при существовании как минимум одного интеграла движения. Окончательно, в особом случае ограниченной задачи трех тел, было обнаружено, что ненулевой эксцентриситет планеты значительно уменьшает точность интегрирования, в сравнении с круговым случаем.

6. Литература

1. J. Baumgarte, Numerical stabilization of all laws of conservation in the many body problem. Celestial Mechanics 8 (1973). P. 223-228.

2. J. Baumgarte, Numerical stabilization of the differential equation of Keplerian motion. Celestial Mechanics 5 (1972). P. 490-501.

3. P. E. Nacozy, The use of integrals in numerical integrations of the N-body problem. Astrophysics and Space Science 14 (1971) 40-51.

4. M. A. Murison, On an efficient and accurate method to integrate restricted three-body orbits. Astron. J, 97 (5), May 1989. P. 1496-1509.