Математическое моделирование тесных сближений малых тел Солнечной системы с большими планетами и Луной Текст научной статьи по специальности «Математика»

УДК 521.1; 523.642 В. В. Абрамов

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕСНЫХ СБЛИЖЕНИЙ МАЛЫХ ТЕЛ СОЛНЕЧНОЙ СИСТЕМЫ С БОЛЬШИМИ ПЛАНЕТАМИ И ЛУНОЙ

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

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

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

Поскольку моменты тесных сближений являются кратковременными по сравнению с общей продолжительностью вычислительного процесса, применение методов Адамса является вполне оправданным для численного интегрирования уравнений движения малых тел при отсутствии тесных сближений. Это связано с высокой эффективностью многошаговых методов. Например, с помощью неявного метода Адамса—Мултона 11-го порядка можно добиться практически такой же точности, что и в модифицированном методе Эверхарта 27-го порядка. При отсутствии тесных сближений небесных тел друг с другом процесс численного интегрирования методом Адамса—Мултона с шагом 0,2 дня происходит в 4-5 раз быстрее, чем методом Эверхарта с шагом в 1 день [2]. Поэтому метод Адамса—Мултона был выбран для решения уравнений движения малых тел в отсутствие сближений.

Данный многошаговый метод является неявным, следовательно, при его использовании необходимо решать нелинейное уравнение, например, с помощью итерационного процесса [1]:

Р : уП+1 = Уп + Ь1-Вг/(Хп-г, Уп-г), г =0

Е: П1 = / Х+ь у™),

к

С '■ У(п11 = Уп + шо/п+1 + й £ М{/(Хп+1-г, Уп+1-г),

г=1

. (1)

Е : /п++1 = / (Хп+Ь У|г+11))>

к

с: У(+)1 = Уп + йМо /(++1 + й £ М{/(Хп+1-г, Уп+1-г) = Уп+1,

г=1

Е : /п+1 = /{Хп+1, Уп+1)•

Процесс вычислений (1) называется методом прогноза и коррекции Р(ЕС)УЕ. Запись Р означает применение предсказывающей формулы метода Адамса—Бэшфорта, Е — вычисление правой части дифференциального уравнения, с — применение исправляющей формулы метода Адамса—Мултона; V —количество итераций, Б{ и Мг- — коэффициенты метода Адамса—Бэшфорта

и метода Адамса—Мултона соответственно. Таким образом, сначала производится вычисление начального приближения У(г°+1 по формуле метода Адамса—Бэшфорта, которое затем используется в вычислении правой части дифференциального уравнения — функции /, входящей в формулу метода Адамса—Мултона. Итерации продолжаются до тех пор, пока не будет достигнута сходимость, либо фиксированное число раз.

Если продолжать численное интегрирование с помощью метода Адамса—Мултона в моменты тесных сближений, то необходимо увеличить порядок аппроксимации с 11 до 16, а также уменьшить шаг интегрирования [3].

Проблему сохранения устойчивости также можно решить путём применения в моменты тесных сближений метода Эверхарта [4], который обладает большей устойчивостью по сравнению с многошаговыми методами. Эффективность данного неявного одношагового метода была подтверждена путём сопоставления результатов вычислений координат и скоростей небесных тел с базой данных БЕ405, которая в свою очередь согласована с радиолокационными наблюдениями [5]. Другим преимуществом метода Эверхарта является простота смены шага интегрирования. То есть данный метод, как и любой другой одношаговый, позволяет плавно изменять шаг, например, уменьшать его при увеличении модуля ускорения малого тела.

Для сравнения эффективности методов Адамса и Эверхарта, они применялись к одним и тем же дифференциальным уравнениям движения, описывающим математическую модель движения небесных тел. При вычислениях, помимо взаимного влияния небесных тел, учитывались релятивистские эффекты. Для этого необходимо было использовать дифференциальные уравнения движения в барицентрической системе координат с учётом ньютоновских и шварцшиль-довских членов [6]:

,2^ г ] - г г (п 2к2(в + у) „ Шк к2 (2 в - 1) V- Шк

г г = к X щ —— 1--------2--X —--------2---Е —+

I^I2 П лГ^ I2 2(1 + г). . 3 ((г г- г 1)г г 12 1 , ,.-1

ГЫ +(1+г)Ы --^"гг-+2С2(Г1-гг)г1)+

+ ^2 £ Ш ((гг - г]) х ((2 + 2Г)гг - (1 + 2Т)г]^(гг - г]) + ^С2” ^^7• (2)

1 Иг] 2С 1 'г1

Здесь г г, Г г, Г { — координаты, скорость, ускорение в барицентрической системе координат г-го возмущаемого тела; г], Г], Г] — координаты, скорость, ускорение в барицентрической системе координат 1-го возмущающего тела; к2 — гравитационная постоянная, Ш1 — масса 1-го тела; гч = 1Г1 - Г11, = I г г I; в и у — релятивистские параметры (в = у =1); с — скорость света.

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

Для первого случая был выбран следующий алгоритм определения шага интегрирования:

к°, а ^ а° (метод Адамса),

Н°\/а°, а > а° (метод Эверхарта),

к(а) = 17_ _ ,_____________________________^_______________ч (3)

где

а° = к2 Ш°. (4)

г°

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

а = тах а[, (5)

г=1,2,..., п

где а{ — модуль ускорения, создаваемого г-м телом, влияющим на ]-е тело, движение которого исследуется (г = ]). При этом Солнце в множество объектов а{ не входит.

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

расстояние между исследуемым телом и некоторым телом массой ш°, при котором не нарушается устойчивость, равно г°. Таким образом, величина а° является «критическим» ускорением, при возникновении которого следует начать уменьшать шаг. Например, численное интегрирование методом Адамса—Мултона с шагом /г° = °,2 дня можно производить до тех пор, пока ис-

Тогда а° = 2,22192• 1°-6 — согласно (4), критическое ускорение, при котором необходимо начать уменьшать шаг.

При возникновении максимального ускорения а, равного критическому а°, сообщаемого одним из тел, за исключением Солнца, необходимо переключиться на метод Эверхарта и в дальнейшем считать с переменным шагом Н, согласно второму условию формулы (3), до тех пор, пока ускорение а снова не станет меньше критического а°. В этот момент необходимо вновь продолжить интегрирование с помощью метода Адамса—Мултона.

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

Таким образом, в моменты тесных сближений величину шага интегрирования Н° необходимо уменьшить в / раз. Число / также определяется опытным путём и является достаточно большим. Для очень тесных сближений шаг интегрирования следует уменьшать в пять и более тысяч раз.

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

Для каждого способа численное интегрирование дифференциальных уравнений движения (2) осуществлялось с различным шагом. Закладка начальных данных для старта многошагового метода производилась с помощью метода Эверхарта. Для первого способа значение критического ускорения в формуле (3) было выбрано равным а° = 2,22192 • 1°-6. Для второго способа в формуле (6) величина / = 5°°°.

Сначала численное интегрирование производилось с помощью первого способа с начальным шагом /г°1 = °,2. В результате вычислялись координаты х^), у^), z1(t) астероида АрорЫз для некоторых моментов времени t. Затем интегрирование проводилось заново с другим начальным шагом Н°2 = °,125. В результате определялись координаты Х2(t), У2Ш, Z2(t) несколько отличающиеся от координат х^), у^), Zl(t) для тех же моментов времени. Для сравнения полученных координат вычислялись расстояния между положениями астероида, соответствующих координатам Х1(0, у^), Zl(t) и координатам Х2^), У2Ш, Z2(t), то есть определялись величины

по которым можно судить о точности вычислений. Аналогично определялись величины г (О по формуле (7) и для второго способа.

В таблице для различных моментов времени представлены отклонения г(0 в километрах, вычисленные по формуле (7), для астероида 99942 АрорЫз на интервале интегрирования с 31 декабря 2006 года по 1 июня 2207 года при использовании в моменты тесных сближений либо метода Эверхарта, либо Адамса—Мултона. Как видно из таблицы, при использовании в моменты сближений метода Эверхарта сходимость решения после 2013 года становится хуже, чем при использовании метода Адамса—Мултона. Расхождения г(^ во втором столбце становятся примерно на порядок больше, чем в третьем. Связано это с тем, что в 2013 году произойдет одно из сближений астероида 99942 АрорЫз с Землей, не являющееся тесным [7].

Однако после тесного сближения данного астероида с нашей планетой 13 апреля 2029 года [8] дальнейшего ухудшения сходимости решения методом Эверхарта не происходит. Расхождения г(0 для первого способа продолжают превышать расхождения г(0 для второго способа до конца исследуемого временного интервала примерно на порядок.

На рисунке для астероида 99942 АрорЫз изображены графики логарифмов отклонений г при шаге интегрирования 0,2 и 0,125 дня в зависимости от интервала интегрирования t, представленного в днях, при использовании в моменты тесных сближений либо метода Эверхарта,

следуемый объект не сблизится с Землёй (масса ш° = 3,°°349 • 1° 6 кг) на расстояние г° = °,°2 а. е.

(6)

г и) = г2й))2,

(7)

Отклонения г (км) при использовании различных методов в моменты тесных сближений при шаге 0,2 и 0,125 дня

Дата Метод Эверхарта Метод Адамса

22.06.2012 0,0000017 0,0000015

27.07.2013 0,0000080 0,0000028

31.05.2020 0,0000791 0,0000107

17.09.2026 0,0003519 0,0000369

24.02.2046 76,261015 6,5681913

03.08.2065 2,8712877 0,2488104

10.01.2085 102,50660 8,6936867

20.06.2104 24,886484 2,1058441

28.11.2123 100,97956 8,5690589

07.05.2143 851,69127 72,216922

27.04.2160 2056,1727 174,80432

05.10.2179 4899,2355 415,65118

14.03.2199 78420,181 6705,8011

01.06.2207 127471,42 12048,359

Igr,

км

4

-2

-4

к

\ / \ / ^

v у

г

Р'

О 2000 4000 6000 Г, дней

Логарифм отклонения г (км) при шаге 0,2 и 0,125 дня в зависимости от интервала интегрирования t (дней) при использовании в моменты тесных сближений: 1 — метода Эверхарта; 2—метода Адамса

либо метода Адамса—Мултона.

Из рисунка также видно, что дальнейшего ухудшения сходимости решения методом Эверхарта не происходит до конца исследуемого промежутка, причем кривые изменения расхождения гпосле 2013 года являются практически параллельными. По-видимому, это связано с тем, что критерий изменения шага в методе Эверхарта не достаточно приемлемо работает для умеренных сближений, то есть таких, как сближение астероида АрорЫз с Землей в 2013 году.

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

Несмотря на то, что моменты сближений составляют незначительную по времени часть от суммарного исследуемого интервала, метод Эверхарта сильно замедляет общее время вычислительного процесса. В результате, применение первого из рассмотренных выше способов оказывается менее предпочтительным. Вычисление параметров орбиты астероида 99942 АрорЫз на интервале времени 200 лет с помощью первого способа продолжалось 453 секунды, а с помощью второго — только 283 секунды.

Таким образом, если при исследовании движения малых тел, сближающихся с большими планетами или Луной, решающим фактором является быстродействие, то оправданным является применение метода Адамса—Мултона и при их отсутствии. А метод Эверхарта следует использовать

и в моменты тесных сближений только для закладки начальных значений многошагового метода.

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

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

1. Холл, Дж. Современные численные методы решения обыкновенных дифференциальных уравнений [Текст] / Дж. Холл, Дж. Уатт. — М.: Мир, 1979. — 312 с.

2. Абрамов, В. В. Применение методов Адамса к решению уравнений движения больших планет, Луны и Солнца [Текст] / В. В. Абрамов / Мат. моделирование и краевые задачи: Тр. Третьей Всерос. научн. конф. — Самара: СамГТУ, 2006. — Ч. 3. — С. 13-19. — ISBN 5--7964--0802--X.

3. Абрамов В. В. Эффективность метода Адамса—Мултона при математическом моделировании движения малых тел Солнечной системы [Текст] / В. В. Абрамов / Нелинейный динамический анализ-2007: Тез. докл. междунар. конгр. — СПб.: СПбГУ, 2007.—С. 184.

4. Everhart, E. Implicit single methods for integrating orbits [Text] / E.Everhart // Celestial mechanics. — 1974. — No. 10.— P. 35-55.

5. Заусаев, А. Ф. Применение метода Эверхарта 31-го порядка для решения уравнений движения больших планет [Текст] / А.Ф. Заусаев, А. А. Заусаев, А. Г. Ольхин // Тр. ГАИШ. — 2004. — Т. LXXV: ВАК-2004. —С. 209-210.

6. Newhall, X.X. DE 102: a numerically integrated ephemeris of the Moon and planets spanning forty-four centuries [Text] / X.X. Newhall, E.M. Standish, Jr. and J.G. Williams // Astron. Astrophys. — 1983. — No. 125. — P. 150-167.

7. Абрамов, В. В. Моделирование сближений астероида 99942 Apophis с внутренними планетами и Луной [Текст] / В. В. Абрамов, С. С. Денисов, Л. А. Соловьев / Тез. докл. XXXII Самарской обл. студ. научн. конф. — Самара, 2006. —Ч. 1. —С. 101.

8. Абрамов, В. В. Математическое моделирование движения малых тел Солнечной системы на основе методов Адамса [Текст] / В. В. Абрамов // Вестн. Сам. госуд. техн. ун-та. Сер.: Физ.-мат. науки. — 2006. — № 43. — С. 192194. — ISBN 5-7964-0877-1.

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

Поступила 21.04.2007