textarchive.ru

Главная > Задача


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

Постановка задачи

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

Аналитические методы позволяют получить основные неравенства в движении объектов. В XVIII и XIX веках с их помощью были созданы модели движения больших планет Солнечной системы, предсказано существование Нептуна. Построенные аналитические теории в XX веке перестали соответствовать точности, качеству и количеству наблюдений положений небесных тел (Абалакин и др., 1976). В настоящее время аналитические модели движения в их классическом виде актуальны для двух спутников Марса и главных спутников планет-гигантов. Эти спутники имеют очень малые численные значения эксцентриситетов и углов наклонений орбит к плоскости экватора основной планеты (Емельянов, 2005).

Качественные исследования дают общий взгляд на картину движения объектов выбранного класса. При таком подходе начальная постановка задачи, близкая к реальному положению вещей, претерпевает существенные упрощения. Изучается, например, движение материальной точки в поле притяжения центрального тела (Солнца) и Юпитера, обращающегося вокруг Солнца по эллиптической орбите (Дубошин, 1975). Для малых планет, находящихся в орбитальном резонансе с Юпитером, были получены интересные результаты, проливающие свет на эволюционные характеристики резонансных орбит, но качественные выводы были сделаны в предположении достаточно малых эксцентриситетов и наклонностей астероидов (Герасимов, 1992).

Задача исследования эффектов прохождения астероида в непосредственной близости от большой планеты и задача о движении объектов в случае, когда их орбиты пересекаются в плоскости эклиптики, представляют большие и почти непреодолимые трудности для качественных и аналитических методов. Следует отметить важные результаты по эволюции угла наклонения, долготы восходящего узла, эксцентриситета и аргумента перигелия для малых планет группы Аполлона и Амура, полученные на основе численно-аналитического метода (Вашковьяк, 1986).

Численное интегрирование дифференциальных уравнений движения в настоящее время является самым распространённым способом исследования эволюции элементов орбит и прогноза положений небесных тел Солнечной системы вообще и астероидов в частности (Баканас и др., 2003).

Как и аналитические методы, методы численного интегрирования развивались последовательно вместе с возникновением новых задач науки и техники. У всех методов есть общие черты (аппроксимация полиномами по времени, оценка точности вычислений, шаг численного интегрирования, интервал надёжности расчёта) и существенные отличия (применение разделённых разностей, использование производных, количество шагов, наличие неизвестных величин в правой части нелинейных уравнений). Для решения конкретной задачи и достижения необходимой точности вычислений могут оказаться полезными лишь несколько методов. Если произойдут изменения в постановке задачи, то часть применяемых методов могут оказаться менее точными, причём заранее, теоретически, такое предсказать невозможно. Каждая система дифференциальных уравнений с заданными начальными условиями (задача Коши) требует дополнительных исследований.

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

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

Уравнения движения

В алгоритмах используются следующие единицы измерений:

Единица расстояния – астрономическая единица (AU), 1 AU = 149597870.691 км. Единица времени – эфемеридные сутки. Единица массы – масса Солнца . Квадрат постоянной Гаусса равен гелиоцентрической гравитационной постоянной .

Через обозначим вектор положения малой планеты относительно центра Солнца в системе экватора и эклиптики на стандартную эпоху J2000.0 (календарная дата 2000, 1.5 январь, юлианская дата 2451545.0).

Три дифференциальных уравнения движения в векторной форме имеют вид

,

где

– ускорение, вызываемое Солнцем,

– ускорение, обусловленное притяжением больших планет Солнечной системы и Луны.

Ускорения вычисляются по формулам:

где

– вектор положения большой планеты или Луны относительно Солнца,

–масса возмущающего тела в единицах массы Солнца.

Положения планет как функции эфемеридного времени известны и заданы численными эфемеридами DE405/LE405 (Standish at all, 1998). В этом же отчёте даны числовые значения отношений массы Солнца к массам планет.

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

Одношаговые методы

Большой класс методов численного интегрирования составляют одношаговые методы. Дифференциальное уравнение заменяют на нелинейное разностное уравнение (Бронштейн и Семендяев, 1980), а решение уравнения представляют в виде полинома по степеням шага численного интегрирования:

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

В одношаговом методе Рунге-Кутта четвёртого порядка относительно вектора положения и скорости на следующем шаге определены формулами (Handbook of Mathematical Functions, by M. Abramowitz and I. A. Stegun, 1965):

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

Многошаговые методы

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

Многошаговые метод Адамса, метод Штёрмера, метод Буллирша, например, являются явными, иногда их называют методами экстраполяционного типа. Метод Коуэлла, с другой стороны, является неявным, то есть методом интерполяционного типа. Вначале вектор положения определяют на основе численных значений вектора скорости и вектора ускорения в момент (предсказание), а затем предсказанное значение используют для вычисления конечных разностей и нахождения нового значения (исправление). Это обстоятельство определяет преимущества того или иного способа вычислений. Методы экстраполяционного типа работают надёжно в случае гладкости правых частей уравнений и быстрого убывания численных значений конечных разностей. Неявные методы содержат на несколько операций больше за счёт процедуры предсказания положения, но с помощью этой процедуры удаётся подправить конечные разности в случае резких изменений величины ускорения.

Другие методы

Среди новых способов численного интегрирования, появившихся в последние десятилетия, следует отметить метод Эверхарта (Everhart, 1974). Профессор Эверхарт в своей статье назвал его неявным одношаговым методом и применил для исследования движения короткопериодических комет, испытывающих сближения с Юпитером. Подробный разбор алгоритма выполнен в работе (Татевян и др., 1996). Оказалось, что алгоритм вычислений, не требуя информации о положениях объекта в предыдущие моменты времени, объединяет все преимущества многошаговых методов (конечные разности высокого порядка в численном виде) и неявных методов (предсказание и коррекция). Необходимость дополнительных обращений к процедуре вычисления вектора правых частей полностью компенсируется специальным неравномерным разбиением шага интегрирования. Более того, поскольку результатом очередного шага являются многочлены, аппроксимирующие решение на интервале этого шага, то вектор положения и вектор скорости могут быть вычислены в любой точке внутри этого интервала. В методах Рунге-Кутта, Адамса, Коуэлла, например, результатом очередного шага интегрирования являются численные значения независимых переменных в момент времени, отстоящий от предыдущего момента на длину шага. Для нахождения решения в произвольный момент времени приходится усложнять алгоритм и использовать интерполяцию значений в узлах интегрирования.

Хорошей внутренней точностью отличаются такие одношаговые методы интегрирования, в которых на каждом шаге используются аналитические формулы для вычисления производных высокого порядка от решения дифференциальных уравнений (Яров-Яровой, 1974). У такого подхода есть свои трудности. В случае близкого прохождения двух тел, например, резко уменьшается взаимное расстояние между объектами, а эта величина в возрастающей степени входит в знаменатели всех производных, начиная со второй. На данном, достаточно малом участке орбиты, ряд Тейлора , составленный на основе производных, окажется расходящимся (Штифель и Шейфеле, 1975).

Были разработаны и другие модификации, реализующие идеи нескольких способов численного интегрирования. Один из таких алгоритмов улучшает сходимость и точность метода Эйлера за счёт использования степенных разложений более высокого порядка, конечных разностей и схемы предсказание-исправление (PECE algorithm, predict-evaluate-correct-evaluate) . Алгоритм был назван именем Эрмита.

Производные по независимой переменной от решения имеют вид

На момент известны вектор положения и вектор скорости , вычисляем вектор ускорения и вектор . С помощью степенного разложения выполняем предсказание вектора положения и вектора скорости на момент :

Вычисляем вектора и . На основе этих промежуточных величин определяем исправленные значения векторов положения и скорости на один шаг вперёд:

В третьих слагаемых этих соотношениях использованы конечные разности.

В последней декаде двадцатого века одношаговый метод Эйлера второго порядка был использован в схеме симплектического интегратора (Wisdom, 1991), когда правые части уравнений движения зависят только от координат.

Три строки алгоритма Эйлера можно интерпретировать следующим образом:

1) - скорость объекта изменяется на интервале в полшага под действием гравитационных сил, вычисленных в начальный момент ;

2) - объект движется по орбите на интервале один шаг;

3) - скорость объекта изменяется на интервале в полшага под действием гравитационных сил, вычисленных в момент .

На основе этих простейших соотношений автор статьи (Construction ofhigherordersymplecticintegrators, by Haruo Yoshida, 1990, Phys. Lett. A150, 262) разработал оригинальный метод интегрирования, получивший название метода Йошида.

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

Симплектический метод с успехом применяют для интегрирования уравнений движения N-тел, когда по условиям задачи не нужна высокая точность расчётов (Duncan at all, 1998). Варианты сближения двух или нескольких материальных точек обрабатываются с помощью других алгоритмов (Fabbio Migliorini, 1998).

Результаты вычислений

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

В качестве первого примера в таблице представлены расчёты, выполненные для малой планеты Тутатис на интервале времени 3 года (2001, 2002, 2003 годы). На этом промежутке астероид находится далеко от больших планет и, в частности, от Земли. Вычисления проводились для методов четвёртого порядка относительно малого параметра – шага интегрирования. Для каждого метода выполнялись предварительные расчёты для выбора наилучшего начального шага интегрирования. В последней строке таблицы приводятся оценки для метода Эверхарта 15 порядка. Время вычислений в соответствующей колонке содержит относительные величины. За единицу принято время вычислений с помощью интегратора Эверхарта. Метод Йошида имеет достаточную точность и большую величину шага интегрирования.

Таблица. Астероид Тутатис, интервал интегрирования 3 года, сравнение методов.

метод

шаг (сутки)

ΔR (AU)

время вычислений

Рунге-Кутта (4)

0.1

2.32e-06

10.0

Адамса (4)

0.5

8.82e-08

0.1

Эрмита (4)

1.0

1.54e-07

0.1

Йошида (4)

10.0

6.86e-08

0.005

Эверхарта (15)

1.0

2.01e-15

1.0

В следующей таблице выполнены оценки методом вложения для той же малой планеты Тутатис на интервале времени 8 лет, с 2000 по 2007 годы. В этом промежутке малая планета приближалась к Земле.

Таблица. Астероид Тутатис, интервал интегрирования 8 лет, сравнение методов.

метод

шаг (сутки)

ΔR (AU)

время вычислений

Рунге-Кутта (4)

0.1

1.3e-07

40.0

Адамса (4)

0.5

1.5e-06

2.0

Эрмита (4)

1.0

4.4e-05

0.9

Йошида (4)

10.0

1.6e+00

0.005

Эверхарта (15)

1.0

6.9e-14

1.0

Увеличение времени интегрирования методом Рунге-Кутта обусловлено автоматическим уменьшением шага интегрирования при сближении объекта с Землёй. Метод Йошида даёт неправильные качественные и количественные результаты. При уменьшении длины шага результат не меняется. Теоретически это обстоятельство обусловлено тем, что симплектический алгоритм был получен при условии малости возмущающих сил по сравнению с притяжением центрального тела.

Выполним сравнение численных методов интегрирования высоких порядков. В таблице использованы сокращения: А-Б (12) - метода Адамса-Башфорта 12 порядка; К (12) – метод Коуэлла 11 порядка; Б-Ш – экстраполяционный метод Буллирша-Штёра; Йо (8) – симплектический интегратор Йошида 8 порядка; Э (15) – метод Эверхарта 15 порядка

Вычислительный код многошагового метода Адамса-Башфорта представлен в книге «Астрономия на персональном компьютере» (Монтенбрук и Пфлегер, 2002). Авторы отмечают, что метод был разработан (Shampine and Gordon, 1975) и обладает рядом полезных свойств: количество используемых предшествующих точек не задаётся, а самостоятельно выбирается программой; величина шага является переменной; код написан как самостартующая процедура и не требует задания нескольких начальных точек; значение вектора решения можно получить для любого момента времени.

Алгоритм метода Коуэлла 11 порядка дан, например, в книге «Линейная и регулярная небесная механика» (Штифель и Шейфеле, 1975).

Вычислительный код экстраполяционного метода Буллирша-Штёра заимствован из свободно распространяемого пакета программ “Mercury” (John E. Chambers and Fabbio Migliorini).

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



Скачать документ

Похожие документы:

  1. ТЕЗИСЫ ДОКЛАДОВ ХХXVII САМАРСКОЙ ОБЛАСТНОЙ СТУДЕНЧЕСКОЙ НАУЧНОЙ КОНФЕРЕНЦИИ

    Тезисы
    ... на основе 4-азагомоадамантана проявляют антивирусную активность. В продолжение проведенного исследованияпо ... метода Адамса на больших интервалах времени. В качестве примера абсолютной неустойчивости задачи Коши рассмотрено дифференциальноеуравнение ...
  2. ТЕЗИСЫ ДОКЛАДОВ ХХXVII САМАРСКОЙ ОБЛАСТНОЙ СТУДЕНЧЕСКОЙ НАУЧНОЙ КОНФЕРЕНЦИИ

    Тезисы
    ... на основе 4-азагомоадамантана проявляют антивирусную активность. В продолжение проведенного исследованияпо ... метода Адамса на больших интервалах времени. В качестве примера абсолютной неустойчивости задачи Коши рассмотрено дифференциальноеуравнение ...
  3. Проекты российской академии наук для участия в реализации направлений технологического прорыва

    Документ
    ... массив по небольшому элементу числу его элементов. На модельных примерах показана большая эффективность метода, назадачах ... прогресс в исследованиях КГ реального времени и увеличение мощности графических процессоров позволил технологиям виртуальной и ...
  4. Отчетный доклад президиума российской академии наук

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

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

Другие похожие документы..