Методы интерполяции, используемые для получения координат и элементов орбит больших планет и малых тел Cолнечной системы Текст научной статьи по специальности «Механика и машиностроение»

Небесная механика и астрономия

УДК 521.182

МЕТОДЫ ИНТЕРПОЛЯЦИИ, ИСПОЛЬЗУЕМЫЕ ДЛЯ ПОЛУЧЕНИЯ КООРДИНАТ И ЭЛЕМЕНТОВ ОРБИТ БОЛЬШИХ ПЛАНЕТ И МАЛЫХ ТЕЛ СОЛНЕЧНОЙ СИСТЕМЫ

А. Ф. Заусаев, Д. А. Заусаев

Самарский государственный технический университет,

443100, Самара, ул. Молодогвардейская, 244

E-mail: Zausaev_AFamail.ru

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

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

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

Исследование эволюции орбит малых тел Солнечной системы является актуальной задачей в связи с решением проблемы астероидной опасности. В настоящее время путём численного интегрирования уравнений движения малых тел Солнечной системы создан банк данных координат, скоростей и элементов орбит астероидов групп Аполлона, Амура, Атона, содержащий около четырёх тысяч объектов и 190 короткопериодических комет на интервале времени с 1800 по 2204 гг. [1,2]. Данные о координатах, скоростях и элементах орбит вышеуказанных объектов приведены в банке данных с шагом сто дней. Для получения информации на любой момент времени в промежутке с 1800 по 2204 гг. необходимо использовать интерполяционные методы, например, простейшие методы численного интегрирования: метод трапеций, метод парабол, интерполирование с помощью полиномов Чебышева и т.д.

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

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

Заусаев Анатолий Фёдорович — профессор кафедры прикладной математики и информатики Самарского государственного технического университета; д.ф.-м.н., профессор. Заусаев Дмитрий Анатольевич — студент Самарского государственного технического университета, специальность «Прикладная математика и информатика».

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

Численное интегрирование уравнений движения методом трапеций. Рассмотрим численное интегрирование уравнений движения методом трапеций на примере системы дифференциальных уравнений в барицентрических координатах. Дифференциальные уравнения движения нами выбраны в следующем виде [3]:

где Д2 = (Хг — X)2 + (У — У)2 + — 2)2; гог — эффективный радиус ¿-того

тела; аог — соответствующее ускорение для ¿-того тела на расстоянии Гог от центра массы; X, У, 2 — барицентрические координаты возмущаемого тела; Хг, Уг, —барицентрические координаты возмущающих тел.

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

Предлагается следующий метод численного интегрирования системы (1). Построение метода численного интегрирования рассмотрим на примере решения системы дифференциальных уравнений вида

(1)

(2)

где X = (х, у, г).

Запишем систему (2) в виде

(3)

Интегрируя её в пределах от ¿о до і

(4)

получим

По формуле (5) вычисляется значение скорости объекта на следующем шаге интегрирования. Для получения координат запишем формулу (5) в виде

(X =

Интегрируя (6), имеем

(Х\ М )

+ Р(х,у, г,і)(і

¿0

Но

(і.

(6)

(X =

'¿о

'¿о

(Х\

(і ) ¿0

+

Р(х, у, г, і)(і

'¿0

(і,

отсюда

X = Х*о +

—Х^ (і — іо) + / [ Р(х,у, г,і)(1і(1і.

(і ) ¿о ііоі¿0

(7)

(8)

В зависимости от вида аппроксимации функции ¥(х,у,г,Ь) получаем различные методы численного интегрирования.

Рассмотрим первый метод, предполагая, что функция ¥(х,у,г,Ь) на каждом шаге изменяется линейно. В этом случае метод численного интегрирования аналогичен алгоритму вычисления интегралов по методу трапеций.

Для нахождения координат на следующем шаге вычисляем правую часть уравнения (2). По формуле

Ъ2 X(£о + Ъ) = X(¿о) + X(¿о) ■ Ъ + Х(£о) ■ —

для каждой координаты определяем положение тела на следующем шаге, где Ъ — шаг интегрирования. Затем находим правые части уравнения (2) в точке ¿о + Ъ. Далее вычисляем средние значения правых частей дифференциальных уравнений (2) по формуле

Хер (і) =

Х(іо) + Х(іо + Ь) 2 '

И, наконец, уточняем значения координат в точке іо + Ь по формуле

Ь2

2

Х(і0 + Ь) = Х(і0) + Х(і0) • Ь + Хср(і) •

Для нахождения скоростей используются формулы

Х(іо + Ь) = Х(іо) + ХСр(і) • Ь,

где X = (X,у,г). Таким образом, происходит пошаговый процесс интегрирования системы дифференциальных уравнений (2).

Численное интегрирование уравнений движения методом парабол. Для

численного интегрирования системы (2) методом парабол необходима информация о координатах и скоростях в трёх точках на каждом шаге интегрирования. Для этого по методу трапеций с шагом ^ находим положение координат

і

I

I

I

и скоростей на моменты времени іо, іо + §, іо + Ь. Затем вычисляем значения функции Х(і) по формуле

= Х^(іо) + 4Х (іо + |) + Х(іо + Ь) (9)

6

Координаты и скорости на следующем шаге находим по следующим формулам:

Н2

X(іо + Ь)= X(іо) + X(іо) • Ь + Х(і) • —, (10)

X(іо + Ь) = X(іо)+ X(і) • ь. (11)

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

Опишем основную идею работы алгоритма данного метода. Основной

вклад в возмущающее действие небесных тел оказывает Солнце. Следова-

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

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

По формулам

X * = X — X0, X * = X — X0 (12)

переходим из барицентрической экваториальной системы координат в гелиоцентрическую экваториальную систему координат, где X* — гелиоцентрические экваториальные координаты планет, Луны и исследуемого объекта; X — барицентрические экваториальные координаты планет, Луны и исследуемого объекта, Xо — координаты Солнца в барицентрической экваториальной системе координат.

По формулам (1) вычисляем возмущающие ускорения исследуемого объекта от больших планет на момент времени іо.

По координатам и скоростям на момент времени іо находим элементы орбит больших планет, Луны и исследуемого объекта по известным формулам [5]. Затем смещаем положение всех планет и исследуемого объекта по орбите на шаг интегрирования, т. е. на момент іо + Ь, по формуле

М = Мо + п(і0 + Ь), (13)

где Мо —средние аномалии планет, Луны и объекта на момент времени іо; М — средняя аномалия вышеупомянутых тел на момент іо + Ь; п — среднесуточное движение этих же тел.

По элементам орбит находим координаты и скорости X* +§, Xt*о+h больших планет, Луны, исследуемого объекта на момент іо + Ь по формулам, приведённым в работе [5].

Определяем возмущающие ускорения для каждой планеты, Луны и исследуемого объекта на момент времени іо + Ь по формуле (1).

Находим среднее значение возмущающих ускорений для планет, Луны и объекта на момент времени іо + §.

Вычисляем возмущённые координаты планет, Луны и объекта на момент времени іо + Ь путём учёта взаимного возмущения планет друг на друга и на исследуемый объект по формулам:

Переходим от гелиоцентрической экваториальной системы координат к барицентрической экваториальной системе координат по формулам:

X = X * + Xo, X = X * + Xo. (15)

Таким образом, получены барицентрические экваториальные координаты больших планет, Луны, Солнца и исследуемого объекта на момент времени to + h. Далее процесс вычислений повторяется. Так осуществляется интегрирование дифференциальных уравнений методом оскулирующих элементов.

Исследование эффективности методов трапеций, парабол и оскулирующих элементов. Эффективность методов трапеций, парабол и оскулирующих элементов исследовалась на примере двух астероидов Aten 2004 FU162 и Apollo 2004 TD18. Рассматривались интервалы времени, на которых происходили тесные сближения астероида с Землёй, а также на которых они отсутствовали. Численное интегрирование проводилось для каждого объекта на интервале времени 100 дней, что соответствует периоду, с которым заданы координаты, скорости и элементы орбит для каждого индивидуального астероида. Результаты вычислений координат и скоростей астероида сопоставлялись с данными, полученными методом Эверхарта 27 порядка. Абсолютной погрешностью в координатах и скоростях применяемого метода считалось максимальное отклонение от результатов, полученных методом Эверхарта.

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

За приемлемую погрешность принималось отклонение в координатах X и в скоростях X по отношению к точным значениям на величину, равную 0,000005 а. е. и 0,0000005 а. е./сут. соответственно. Вышеуказанные значения погрешности соответствуют ошибкам телескопических наблюдений.

V*

Xt0 +h

= Xt*o +

d2X * dt2

В качестве начальной точки отсчёта, где астероиды не имеют тесного сближения с Землёй, был выбран момент времени 2415200,5 J.D. (юлианская дата), что соответствует 30 июня 1900 года.

В табл. 1 приведены расхождения в координатах при различных шагах между исследуемыми методами и методом Эверхарта 27 порядка для астероида Aten 2004 FU162 на момент времени 2415300,5 J.D.

Т аблица 1

Абсолютные погрешности в координатах астероида Aten 2004 FU162

Метод расчёта Шаг (д.) Ах (а.е.) Ау (а.е.) е. (а. А

Метод трапеций 0,01 0,00000013 0,00000036 0,00000007

Метод трапеций 0,1 0,00000235 0,00000163 0,00000071

Метод трапеций 0,5 0,00006241 0,00003222 0,00001623

Метод парабол 0,01 0,00000013 0,00000037 0,00000007

Метод парабол 0,1 0,00000166 0,00000176 0,00000076

Метод парабол 0,5 0,00004544 0,00003525 0,00001651

Метод парабол 1 0,00018380 0,00013973 0,00006977

Метод оск. элемент. 1 0,00000017 0,00000034 0,00000006

Метод оск. элемент. 5 0,00000028 0,00000033 0,00000005

Метод оск. элемент. 10 0,00000065 0,00000031 0,00000003

Метод оск. элемент. 50 0,00001092 0,00000033 0,00000042

Как следует из результатов, приведённых в табл. 1, метод трапеций даёт приемлемую точность при шаге 0,1 день, метод парабол также при шаге 0,1 день, а метод оскулирующих элементов обеспечивает необходимую точность при шаге 10 дней.

В табл. 2 приведены расхождения в координатах при различных шагах между исследуемыми методами и методом Эверхарта 27 порядка для астероида Apollo 2004 TD18 на момент времени 2415300,5 J.D.

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

Т аб л и ца 2

Абсолютные погрешности в координатах астероида Apollo 2004 TD18

Метод расчёта Шаг (д.) Ах (а.е.) е. (а. у А Az (а.е.)

Метод трапеций 0,1 0,00000007 0,00000019 0,00000003

Метод трапеций 0,5 0,00000279 0,00000196 0,00000083

Метод трапеций 1 0,00001170 0,00000871 0,00000334

Метод трапеций 5 0,00029257 0,00022635 0,00008440

Метод парабол 0,1 0,00000010 0,00000023 0,00000001

Метод парабол 0,5 0,00000186 0,00000096 0,00000046

Метод парабол 1 0,00000801 0,00000470 0,00000187

Метод парабол 5 0,00020360 0,00012477 0,00004706

Метод оск. элемент. 5 0,00000018 0,00000026 0,00000000

Метод оск. элемент. 10 0,00000018 0,00000021 0,00000002

Метод оск. элемент. 20 0,00000019 0,00000002 0,00000012

Метод оск. элемент. 25 0,00000005 0,00000011 0,00000013

Метод оск. элемент. 50 0,00000103 0,00000163 0,00000093

В табл. 3 приведены значения погрешностей в координатах для астероида Aten 2004 FU162 на момент времени 2453100,5 J.D. Вследствие того, что внутри интервала интегрирования с 2453000,5 J.D. по 2453100,5 J.D. имело место тесное сближение астероида с Землёй на расстоянии 0,000086 а. е., ни один из рассматриваемых методов не даёт удовлетворительной точности при заданных шагах интегрирования. При дальнейшем уменьшении шага результаты несколько улучшаются, однако время расчёта несоразмерно увеличивается.

Таблица 3

Абсолютные погрешности в координатах астероида Aten 2004 FU162

Метод расчёта Шаг (д.) Ах (а.е.) е. (а. А Az (а.е.)

Метод трапеций 0,01 0,00031019 0,00001726 0,00017959

Метод трапеций 0,005 0,00008566 0,00001487 0,00003469

Метод трапеций 0,001 0,00002342 0,00000799 0,00003469

Метод трапеций 0,0005 0,00002153 0,00000738 0,00003232

Метод парабол 0,0001 0,00002130 0,00000729 0,00003200

Метод парабол 0,01 0,00011984 0,00002983 0,00013913

Метод парабол 0,005 0,00005024 0,00001565 0,00006398

Метод парабол 0,001 0,00002230 0,00001487 0,00003335

Метод парабол 0,0005 0,00002179 0,00000747 0,00003266

Метод парабол 0,0001 0,00002129 0,00000729 0,00003199

Метод оск. элемент. 0,01 0,00027231 0,00006325 0,00030166

Метод оск. элемент. 0,005 0,00006964 0,00002702 0,00009983

Метод оск. элемент. 0,001 0,00002310 0,00000833 0,00003417

В табл. 4 приведены значения погрешностей в координатах для астероида Apollo 2004 TD18 на момент времени 2514100,5 J.D. Внутри интервала интегрирования 6 марта 2171 г. произойдёт сближение этого астероида с Землёй на расстоянии 0,000297 а. е. Однако это сближение существенно не отразится на результатах прогноза движения с помощью рассматриваемых методов в пределах 100 дней при шаге интегрирования 0,005 дня.

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

Таблица 4

Абсолютные погрешности в координатах астероида Apollo 2004 TD18

Метод расчёта Шаг (д.) Ах (а.е.) е. (а. А Az (а.е.)

Метод трапеций Метод трапеций Метод трапеций Метод парабол Метод парабол Метод парабол Метод оск. элемент. Метод оск. элемент. Метод оск. элемент. 0,05 0,01 0,005 0,05 0,01 0,005 0,05 0,01 0,005 0,00035603 0,00000027 0,00000004 0,00009784 0,00000021 0,00000017 0,00029640 0,00000143 0,00000015 0,00052909 0,00000786 0,00000282 0,00009602 0,00000436 0,00000192 0,00053217 0,00000810 0,00000297 0,00001556 0,00000084 0,00000046 0,00000060 0,00000062 0,00000039 0,00044770 0,00000090 0,00000045

сближений астероида с Землёй). Однако при сближении астероида с Землёй на расстояние менее чем 0,0002 а. е., более целесообразно использовать метод Эверхарта 27 порядка точности, чем рассматриваемые методы.

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

Работа выполнена при финансовой поддержке Федерального агентства по образованию (проект № РНП.2.1.1.1689).

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Заусаев, А. Ф. Каталог орбитальной эволюции короткопериодических комет с 1800 по 2204 гг. [Текст] / А. Ф. Заусаев, А. А. Заусаев.—М.: Машиностроение-1, 2007.—410 с.

2. Заусаев, А. Ф. Каталог орбитальной эволюции астероидов, сближающихся с Землёй с 1800 по 2204 гг. [Текст] / А. Ф. Заусаев, В. В. Абрамов, С. С. Денисов. —М.: Машиностроение-1, 2007. — 607 с.

3. Заусаев, А. Ф. Теория движения n материальных тел, основанная на новом принципе взаимодействия [Текст] / А. Ф. Заусаев // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. — 2006. — № 43. — С. 132-139.

4. Standish E. M. JPL Planetary and Lunar Ephemerides, DE405/LE405 [Text] / E. M. Standish // Jet Prop Lab Technical Report. IOM 312. F-048. —1998. — P. 1-7.

5. Абалакин, В. К. Справочное руководство по небесной механике и астродинамике [Текст] / В. К. Абалакин, Е. П. Аксенов, Е. А. Гребенников и др. —М.: Наука, 1976. — 862 с.

Поступила в редакцию 24/IX/2008; в окончательном варианте — 30/IX/2008.

MSC: 85-08

INTERPOLATION METHODS USED TO OBTAIN ORBIT COORDINATES AND ELEMENTS OF THE SOLAR SYSTEM LARGE PLANETS AND SMALL BODIES

A. Ph. Zausaev, D. A. Zausaev

Samara State Technical University,

443100, Samara, Molodogvardeyskaya str., 244

E-mail: Zausaev_AFamail.ru

Various methods of interpolation for receiving coordinates, velocities and elements of large planets and small bodies orbits of the Solar system are studied. New osculating elements method is introduced. This method has proved to be the most effective in comparison with trapezium and parabola methods.

Key words: numerical integration, differential equation of motion, trapezium method, parabolas method, osculation elements method.

Original article submitted 24/IX/2008; revision submitted 30/IX/2008.

Zausaev Anatoliy Phedorovich, Dr. Sci. (Phis. & Math.) Prof., Dept. of Applied Mathematics and Computer Science of Samara state Technical University.

Zausaev Dmitriy Anatolievich, Student of Samara state Technical University.