Численный эксперимент. В рассматриваемой работе Дж. Баумгартабыли сравнены 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 интегралов системы были использованы.
Все решениясравнивались с наиболее точным стандартным решением. Кроме того, была выполненапроверка по результатам прямого и обратного интегрирования. Как показывают результаты,