CC BY
71
А.Ф. Заусаев, А. С. Исуткин
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ БОЛЬШИХ ПЛАНЕТ (МЕРКУРИЙ-НЕПТУН) И ЛУНЫ МЕТОДОМ ТЕЙЛОРА
Получены дифференциальные уравнения движения больших планет (Меркурий-Нептун) и Луны с учетом релятивистских эффектов и вращения системы координат. Предлагается метод разложения в ряд Тейлора для интегрирования уравнений движения. На интервале времени 17502050 гг. создан банк данных координат и скоростей больших планет и Луны с использованием высокоточных радиолокационных измерений внутренних планет.
Разработка высокоточной численной теории движения больших планет (Меркурий-Нептун) включает в себя решение следующих основных задач: а) получение
дифференциальных уравнений движения; б) выбор метода решения; в) сопоставление результатов вычислений с наблюдениями. Каждый из этих этапов не является в настоящее время окончательно изученным, поэтому актуальность данной работы не вызывает сомнений.
Целью данной работы является построение численной теории движения больших планет (Меркурий-Нептун) на интервале времени 1750-2050 гг., согласующейся как с радиолокационными, так и с телескопическими наблюдениями. Новизна работы заключается в построении метода численного интегрирования уравнений движения, позволяющего строить решение путем разложения правых частей дифференциальных уравнений в ряды Тейлора.
Практическая ценность работы определяется возможностью использования полученного банка данных координат и скоростей больших планет в задачах по исследованию эволюции орбит малых тел Солнечной системы.
Дифференциальные уравнения движения в планетной задаче
Дифференциальные уравнения движения в ньютоновской задаче п тел в прямоугольных координатах с началом в центре Солнца могут быть приведены к одному уравнению матричного вида [12]:
^ = -к2 • (1 + m) • ^ + £k2 • m1 • (^^ _X.), (1)
г і Дз г
где Х - матрица-столбец с элементами х, у, і; X, - матрица-столбец с элементами х, ,у, т, х, у,
2 - масса и гелиоцентрические координаты возмущаемой планеты; т,, х,, у, 2 - массы и гелиоцентрические координаты больших планет; г, Д і, Гі - расстояния, вычисляемые по формулам: г2=х2+у2+і2, А2=(хгх)2+(угу)2+(2Г2)2, г2=х2+у2+22.
Уравнения (1) является нерелятивистскими уравнениями движения. Однако в середине XIX века Леверье впервые установил, что теоретическое вековое смещение перигелия Меркурия, полученное на основании решения уравнения (1) отличается от наблюдаемого. Полученный результат дал основание усомнится в достоверности нерелятивистских уравнений (1) и побудил исследователей к построению новых теорий движения больших планет.
Разработанная Эйнштейном общая теория относительности позволила в дальнейшем создать релятивистскую небесную механику, свободную от вышеуказанного противоречия. Гелиоцентрические релятивистские уравнения движения планет с учетом лишь ньютоновских и шварцшильдовских членов можно записать в виде [3]
^=_к2 • (1+т) • X+Х*2 • т • (хдх -х)+
^ ґ і Д г
к2 к2 X2 XX XX
+—2((4 - 2а)—X - (1+a)-^X+3а— X+(4 - 2а)—X), (2)
с г гг г
где а - параметр, характеризующий выбор системы координат (для стандартной прямоугольной
системы координат а=1); с - скорость света в пустоте; X2 = х2 + у2 + 22; XX = хх + уу + 22 ;
х, у, 2 - проекции скорости. Уравнения (1) и (2) записаны в гелиоцентрической системе
координат, т. е. с началом в центре Солнца. Как известно, система координат, связанная с центром Солнца не является инерциальной и имеет вращательное движение в Галактике, совершая полный оборот примерно за 200 млн лет. Уравнения движения во вращающихся осях
можно вывести для фиксированных осей путем добавления центробежных и кориолисовых членов. Уравнения (2) во вращающейся системе примут следующий вид:
х , 2 х - X х. к2 „ „к2
x=-k2 • (1+m) • — +^k2 • mt • (-—-3) + -((4-2a)—x-
(3)
X2 XX xx 2
-(1+a)—r x+3a^— x + (4-2a)—x)-w x-2wy ; r r r
y=-k2 • (1+m)^ +^k2 • mt ■ (^-l -l) + ^_((4-2a)kJy-r i — r c2 r4
X2 XX XX
- (1 + a)—- y+3ay + (4 - 2a)— y)-w2 y - 2w&; r r r
z = -k • (1+m)^ +£k2• m ■(Z’-3Z-4) + kJ((4-2a)^z-
r i - f c r
X2 XX XX
-(1+a)—3- z+3a—^ z+(4-2a)—5- z), r r r
где w - угловая скорость вращения системы относительно фиксированных осей; X2 = x2 + у2 + z2; XX = xx + yy + zz .
При исследовании движения Луны учитывалось возмущение от несферичности центрального тела. Потенциал для тел произвольной формы вычисляется по следующей формуле [2]:
и = — I1 -jr Jn fr0 Рп (sinj) + ХХ Pr(k)(sin^)[c„k cos k1 + S„k sin ы]|. (4)
r [ n= 2 V r 0 n=2 k=1 V r 0 J
В случае Земли: m — ее масса; r0 — средний экваториальный радиус; j и 1 — широта и долгота соответственно; Рп иР(г к) — полином и присоединенная функция Лежандра, причем
Р z) = JL dn(z2-1)” ; p(k)(z) = (1 - z2)k/2 dp# (5)
2n! dzn dzk
Коэффициенты Jn, Cnk и Snk зональных и тессеральных гармоник зависят от формы и распределения масс. Численные значения этих коэффициентов до шестого порядка включительно взяты из работы [2].
Системы дифференциальных уравнений (1), (2) и (3) с учетом выражения (4) являются основными уравнениями при изучении движения больших планет и луны.
Одним из недостатков нерелятивистских уравнений Ньютона является принцип дальнодействия, допускающий возможность мгновенного воздействия данного тела на сколь угодно большом расстоянии без посредства промежуточной среды.
Первоначально авторами была предпринята попытка в уравнениях (1) учесть конечность распространения воздействия путем введения запаздывания в силовых полях. Уравнения движения в центральном поле с учетом запаздывания взаимодействия были рассмотрены В.И. Зубовым [5]. Им же было отмечено, на наш взгляд, очень важное обстоятельство, что конечность скорости распространения возмущений и воздействий приводит к тому, что орбиты отдельных материальных точек и целые конфигурации оказываются квантовыми.
Дифференциальные уравнения движения с учетом запаздывания не меняют своего вида по сравнению с уравнениями (1), однако координаты планет, входящие в уравнения, должны быть смещены по времени на ri /с. Это напоминает процесс синхронизации часов с помощью световых сигналов, который был впервые предложен Пуанкаре, а в последствии использовался Эйнштейном.
Численное интегрирование уравнений движения девяти больших планет (Меркурий-Плутон) с учетом запаздывания взаимодействия на интервале времени 100 лет (1900-2000) [4] показало, что для планет от Венеры до Плутона при сопоставлении с современными численными теориями движения больших планет [13,14] расхождения в гелиоцентрической долготе и широте на всем исследуемом интервале превосходит 1 угловой секунды. Однако смещения перигелия Меркурия, которое имеет место в общей теории относительности, не обнаружено. Кроме того, учет запаздывания взаимодействия в силовых полях приводит к
вековому изменению больших полуосей планетных орбит. Этот результат побудил отказаться от учета принципа конечной скорости распространения гравитации в уравнениях (1), и в дальнейшем при построении численной теории движения больших планет (Меркурий-Нептун) использовались релятивистские уравнения (3), а для Луны также учитывались соотношения (4),
(5).
Интегрирование уравнений движения
Дифференциальные уравнения (2), т.е. их правые части состоят из двух частей -ньютоновской и шварцшильдовской. Как видно из записи уравнений (2), ньютоновская часть полностью совпадает с правой частью уравнения (1). Нетрудно оценить, что вклад от шварцшильдовских членов, а также сил Кориолиса незначителен по сравнению с ньютоновскими силами, поэтому интегрирование уравнений (1) проводилось с максимально возможной точностью, а члены Шварцшильда и силы Кориолиса учитывались как возмущения.
Решение уравнений (1) в точке х^+Ъ), если оно известно в точке х^о), методом Тейлора осуществляется по формулам
Сх('о) ^ 1 С2х('о) ^2
2 "
х('о + 'і) - х('о) + ~
&
+-
&
+...
(6)
2 3
Сх('о + Н) _ Сх('о) С х(*о) 1 С х(Го).2 ,
+--------— '1 + ------------------о— '1 +...
& & Ж2 1 2 йг1' (7)
В методах Адамса и Коуэлла значения производных вычисляются с помощью интерполяционных формул Ньютона или Стирлинга через разности от правых частей уравнений (1).
Описывая источники ошибок в численных методах, В.Ф. Мячин [10] отмечает, что замена производных разностями значительно снижает точность формул (6) и (7). Сам В.Ф. Мячин предлагает изящный по математической форме способ (впервые предложенный Стефенсоном) для нахождения производных. Однако метод Стефенсона-Тейлора, изложенный в работе [10], применим только для решения уравнения (1).
В предлагаемом алгоритме получены общие формулы для нахождения производных от степенной функции. Это позволяет строить решения численным методом путем разложения в
йх
ряды Тейлора для довольно обширного класса уравнений вида — = /(х), где /(х) - матрица-
йг
столбец из рациональных функций.
Полагая в уравнении (1)
1
= Я:
1
А3
= &
запишем его в виде
С2 х
= АХЯ + £ В, ((х - х)Бг - х,Я,).
(8)
(9)
& 2
ш ,=о
Дифференцируя соотношения (9), выражение для вычисления производных можно записать в следующем виде:
Сп+2 х
Жп+2 ё'п
или, применяя формулу Лейбница, имеем
ёп (х, - х)Б, ёп (х,Я,)
А
С'п
сіі”
Сп+2 х Сіп+2
= аЕс
сп-}х СзЯ
(
+
(Хі - х) ё:1Б1 _1 ё^х, ё]Я
\
&
п - з
&
п - з
(1о)
(11)
у=0 ьи ы .=о ^ у =о ш ш у=о
В силу соотношений (8) Я, и Я являются степенными функциями вида Рщ, где Рщ ■ матрица-столбец с элементами (х2+у2+/2)-3/2, ((х1-х)2+(у1-у)2+(71-/)2)-32, (х12+у12+/12)-32, Ш1= - 3/2.
3
3
г
г
,=о
с
Для нахождения производных по формулам (11) необходимо уметь вычислять
С" Рщ
производные ---------. Для вычисления этих производных были получены следующие
Сг"
рекуррентные соотношения:
С'Рт1 С'Р
—— = т1 • Рт1 -1 • — + т1 • (т1 -1) • Рт1 -2 • ^._п + т, • (т1 -1) • (т1 -2) • Рт1 -3 • ^ +... +
Ж Ж
+ т1 • (т1 -1) • (т1 - 2) •... • (т1 -1 +1) • Рт-' •[Ср\> > (12)
СР я С2 Р к С"Р к к к к к к к
где сТ^о; =^ - ; сt^ = F",o; ^ = ^ ' ^'=^Г+^
где /=2,3,4,- ; у=1,2,3,4,-С"Р
Компоненты с п находятся по формулам
С"(х2 + у2 + 72) С"-'х й1х сп-/у С'у С"-^ С'г
— ------ ------- = £ с "------:---- + £ с "---------- + £ с "--- —-;
Л" £0 Сг"-' Л' £0 аг"-1 Л' £0 аг"-1 Л'
л" ((X/ - х)2 + (у. - у)2 + (г. - г)2) = £с< ^ С"-' (х, - х) ^ С (х, - х) +
Сг" £ " Сг"-' Сг'
+ ус> • С"-' (у1 - у) • С (у1 - у) +ус> • С"-' (7' - г) • С' (г^ - г) _ (13)
£0 " Л"-' Сг' £ " Сг"-' Л' ;
С"(х'2 + у'2 + г'2) = £с, • С% + £с, ^ ^ + £с, ^ С"% С%
Сг" £ " Сг"-' Л' £ " Сг"-' Л' £ " Сг"-' Сг' '
Радиолокационные наблюдения внутренних планет и Марса
Радиолокационные наблюдения планет, проводимые как в нашей стране, так и в США, позволили значительно уточнить величину астрономической единицы, а также построить численные теории движения внутренних планет и Луны, на два-три порядка превосходящие по точности теории, основанные на оптических наблюдениях [6,9]. Радиолокационный метод позволяет измерить топоцентрическое расстояние центра видимого диска планеты и радиальную скорость изменения этого расстояния.
Полагая, что электромагнитное излучение покинуло передающую антенну в момент її , достигло центра видимого диска планеты в момент 1 и, отразившись, вернулось на Землю к приемной антенне в момент 1 , обозначим векторами г1(11 ), г2(12) и г3(13), соответственно, положение передающей антенны, центра масс планеты и приемной антенны в соответствующие моменты времени в гелиоцентрической системе координат и получим, что полное время распространения волн А1 определяется по формуле [1]
А = 1 ( & ) - ((()| + \гз (/3) - Г2 (/2) - 2Д|), (14)
где с - скорость света в вакууме; Я - радиус планеты. Полагая скорость радиального смещения планеты в направлении от передающей антенны на планету равной
т I „ _ „ О
Г12 = ~т7Г2 (/2) — Г1 (/1^ ,
т
а скорость радиального смещения планеты в направлении от планеты на приёмную антенну -
т і
Г23 _ \Г3 (/ 3 ) — Г2(/ 2) ,
т
получаем, что полное доплеровское смещение частоты сигнала определяется по формуле [1]
Д/ = /3 - /, @ + Г2У ). (15)
где / - частота начального излучения, а/3 - частота принятого сигнала.
Результатами измерения являются: время излучения, время запаздывания и доплеровское смещение. На основе наблюдательных данных как в нашей стране, так и в США, сотрудники ИТА РАН составили фонд радионаблюдений для Меркурия, Венеры и Марса [8,11]. Этот фонд в виде распечаток нам любезно был предоставлен Г. А. Красинским.
На основании имеющегося фонда радиолокационных данных для каждой из планет (Меркурий, Венера, Марс) был составлен банк радиолокационных данных, куда помимо вышеуказанных величин входили значения топоцентрических координат передающей и приемной антенны. Для удобства сопоставления наблюдательных и теоретических данных время излучения из всемирного переведено в эфемероидное и представлено в юлианских днях. В качестве единицы измерения временной задержки использовалась 1 мкс, что соответствует 6 мкс ~ 1 км.
Краткое описание программы численного интегрирования уравнений движения больших
планет (Меркурий-Нептун) и Луны на интервале времени с 1750 по 2050 гг.
Построение численной теории движения больших планет можно разбить на ряд самостоятельных этапов: а) определение начальных орбит или нахождение начальных данных координат и скоростей планет; б) численное интегрирование уравнений движения; в) сопоставление координат и скоростей, полученных из вычислений с наблюдениями. Вследствие того, что к моменту начала исследований имелся ряд численных и аналитических теорий движения больших планет, был выбор - определить начальные данные элементов орбит планет непосредственно из наблюдений или, используя начальные данные предшественников, путем коррекции согласовать их с радиолокационными наблюдениями. Авторами был выбран второй путь и за основу взята теория движения больших планет (Меркурий-Нептун), разработанная французским бюро долгот [13], где координаты и скорости вычислялись с помощью полиномов Чебышева, коэффициенты которых хранились на магнитных лентах и дисках ЭВМ БЭСМ-6 НИВЦ МГУ. Из банка данных (ББЬ) были взяты начальные координаты и скорости больших планет.
Выбор начального момента, на который из ББЬ были взяты координаты и скорости больших планет, обусловлен тем, что вблизи него имеются точные радиолокационные наблюдения этих планет. В данном исследовании за начальный момент интегрирования выбрана юлианская дата 2438985,20524 по эфемеридному времени. Эта дата совпадает с радиолокационными наблюдениями Меркурия, и вблизи нее имеются наблюдения Венеры и Марса.
Значения обратных масс планет приведены в табл. 1. Постоянные величины: скорость света, астрономическая единица имели следующие значения: с=299792,456 км/с,
1а.е .=149597877,0 км.
Таблица 1
Обратные массы планет
Меркурий 6023600,0 Марс 3098710,0 Уран 22869,0
Венера 408523,5 Юпитер 1047,355 Нептун 19314,0
Земля 332945,662 Сатурн 3498,5 Луна 27069007,130
В табл. 2 приведены начальные данные экваториальных координат и скоростей больших планет (Меркурий-Нептун) и Луны на начальный момент интегрирования, приведенные к эпохе 1950,0. Номерами с 1 по 8 обозначены гелиоцентрические координаты и скорости планет с Меркурия по Нептун, под номером 9 приведены геоцентрические координаты и скорости Луны.
Имея набор начальных данных и проводя численное интегрирование уравнений (3) с учетом формулы (4), в нашу задачу входило согласовать численную теорию движения больших планет как с радиолокационными наблюдениями, так и с данными ББЬ. Последние, как известно, периодически корректировались на различные моменты времени. Для получения
согласования численной теории как с радиолокационными наблюдениями, так и с ББЬ, был положен параметр ш (угловая скорость вращения системы отсчета), равный - 0,1 • 10-9.
Таблица 2
Исходные данные экваториальных координат и скоростей больших планет (Меркурий-Нептун) и Луны на момент .1.В.=2438995.20524 по эфемеридному времени (эпоха равноденствия 1950.0)
У У г
1 0.2767824843240 -0.2663109230427 -0.1710542991008
2 -0.5562338371513 -0.4334777516141 -0.2624071588727
3 0.7708799972440 -0.6031851105690 -0.2615706859069
4 -0.70541498193 -1.228088178385 -0.5447665024727
5 1.315251047812 4.515557687907 1.905038796205
6 9.245711130382 -2.4992740473 -1.433026177260
7 -17.63514760806 4.333491971666 2.147945785814
8 -19.88966976784 -21.35132182129 -8.245616118940
9 0.002167101739566 -0.001384211412704 -0.000846731754166
X У 2
1 0.01548861824641 0.018320790135613 0.008219892692185
2 0.01277015692272 -0.01398792915761 -0.00710976591121
3 0.01088318279322 0.01195425995239 0.005183893712663
4 0.012937477745 -0.004706404681947 -0.0025076995549196
5 -0.007377463032567 0.002060465996562 0.001064264639332
6 0.001336765652834 0.004941971146370 0.001985552817801
7 -0.001076509155149 -0.003645055901199 -0.001581921432578
8 0.002342947947438 -0.001864424749196 -0.0008226848154466
9 0.000327489552997 0.000425707818095 0.000169080796176
Программа численного интегрирования уравнений движения состоит из основной программы и двух подпрограмм. С помощью первой подпрограммы вычисляются значения пой производной от произведения двух функций, вторая подпрограмма вычисляет п-ую производную от степенной функции по формулам (12). В основной или управляющей программе проводятся необходимые вычисления по формулам (6)-(13), а также по формулам (14) и (15) проводятся сопоставления с радиолокационными наблюдениями. Основной вопрос при суммировании рядов Тейлора - это вопрос о величине радиуса сходимости данного ряда. Если для планет от Венеры до Нептуна коэффициенты ряда Тейлора довольно быстро убывают и при достаточно большом отрезке степенного ряда шаг можно выбрать значительно больше одного дня, то для Меркурия этого сказать нельзя. Лишь коэффициенты начальных членов ряда убывают довольно быстро, затем этот процесс замедляется, так, например, 20-я производная от координат Меркурия составляет около 10-12. Несмотря на то, что ряд Тейлора знакочередующийся, увеличение числа членов, удерживаемых в ряде, не может существенно увеличить шаг интегрирования. В данном исследовании шаг интегрирования равнялся одному дню. Количество членов ряда для вычисления координат и скоростей равнялось 12, что обеспечивало точность 10-13 для Меркурия на каждом шаге интегрирования. В программе проводилось сравнение результатов вычислений с радиолокационными наблюдениями. Это позволило корректировать начальные данные в уравнениях (1) и (2).
В процессе интегрирования координаты и скорости планет на стандартные моменты времени через каждые 10 дней записывались на магнитный диск ПК.
Следует отметить, что для Луны мы не располагали точными наблюдениями, поэтому координаты и скорости стремились согласовать с координатами и скоростями из ББЬ [13].
Таким образом, на интервале времени 1750-2050 гг. с шагом 10 дней был создан банк данных координат и скоростей больших планет (Меркурий-Плутон) и Луны, который предназначен для эффективного использования работы программ численного интегрирования уравнений движения малых тел Солнечной системы.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Александров Ю.Н., Кузнецов Б.И., Петров Г.М., Ржига О.И. Методика радиолокационных астрономических наблюдений // Астрон. журн.. 1972. Т. 49. С. 175-185.
2. Аксенов Е.П. Теория движения искусственных спутников Земли. М.: Наука, 1977. 360 с.
3. БрумбергВ.А. Релятивистская небесная механика. М.: Наука, 1972. 382 с.
4. Заусаев А.Ф. Интегрирование уравнений движения больших планет с учетом запаздывания взаимодействия в силовых полях // Всесоюз. конф.: Методы исследования движения, физика и динамика малых тел Солнечной системы: Всесоюз. конф.: Тез. докл. Душанбе, 1989. 31с.
5. ЗубовВ.И. Аналитическая динамика системы тел. Л.: Изд-во ЛГУ, 1983. 385 с.
6. Кислик М.Д., Комока Ю. Ф., Котельников В.А., Петров Г.М., Тихонов В. Ф. Определение орбит Земли и Венеры, астрономической единицы и радиуса Венеры на основе радиолокационных наблюдений Венеры в 1962-1977 гг // ДАН, 1978. Т.241. № 5. С.1046-1049.
7. Киттель Ч., Наут У., Рудерман К. Механика. М.: Наука, 1975. 282 с.
8. Красинский Г.А. Исследование движения больших планет и Луны на интервале 1715-2000 гг. и уточнение некоторых астрономических постоянных: Автореф. дисс. д-ра наук. Киев, 1988.
9. Красинский Г.А., Питьева Е.В., Свешников М.Л., Свешникова Е.С. Уточнение эфемерид внутренних планет и Луны по радиолокационным, лазерным и меридианным наблюдениям 1961-1980 гг. // Бюл. ИТА. 1982. Т.15. № 3(166). С.145-175.
10. Мячин В.Ф., Сизова О.А. Современное интегрирование уравнений небесной механики численным методом Тейлора-Стефенсона // Бюл. ИТА. 1970. Т.12. № 5(138). С.380-400.
11. Питьева Е. В. Учет топографии Марса и Венеры при обработке радиолокационных наблюдений // Бюл. ИТА. 1982. Т.15. № 3(166). С.169-175.
12. Чеботарев Г.А. Аналитические и численные методы небесной механики. М.:Наука, 1965. 367 с.
13. Bretagnon P. Theorie du mouvement de l’ensemble des planetes, solution VS0P82 // Astron. Astroph. 1982. V.144. P.278-288.
14. Oestervinter C., Cohen C. New orbital elements of Moon and planets // Celest.Mech. 1972. V.5. P.317-395.
