Научная работа: Эволюция орбиты Марса на интервале времени в сто миллионов лет

РОССИЙСКАЯ АКАДЕМИЯ НАУК

ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР им. А.А.ДОРОДНИЦЫНА

ИНСТИТУТ КРИОСФЕРЫ ЗЕМЛИ СО РАН

СООБЩЕНИЯ ПО ПРИКЛАДНОЙ МАТЕМАТИКЕ

Гребеников Е.А., Смульский И.И.

ЭВОЛЮЦИЯ ОРБИТЫ МАРСА НА ИНТЕРВАЛЕ ВРЕМЕНИ В СТО МИЛЛИОНОВ ЛЕТ

Москва-Тюмень, 2007 г.

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

1. Введение

1.1 Формулировка проблемы

орбита марс эволюция амплитуда

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

Сошлемся здесь, например, на исследования эволюции решений уравнений движения больших планет Солнечной Системы на весьма больших интервалах времени, выполненные Ж.Ляскаром, Т.Куинном и другими известными специалистами [1-4]. В этих работах показано, что на промежутках времени порядка 100 млн. лет, 250 млн. лет и так далее, до 15 млрд. лет, элементы орбит отдельных планет (Меркурия, Марса, Венеры, Плутона и др.) достигали значений, существенно больших чем их «начальные» значения. Это позволило упомянутым исследователям сделать вывод о возможной неустойчивости Солнечной Системы и о хаотичном характере движений в ней (некоторые выдержки из работ Ж Ляскара мы приводим в Приложении 1). Эта проблема исследовалась ими на основе применения аналитических методов теории вековых возмущений и методов численного интегрирования осредненных по угловым переменным уравнений движения планет. Кроме того, Т.Куинн и его коллеги [3] решали эту задачу с помощью прямых методов численного интегрирования первоначальных уравнений движения на мощных суперкомпьютерах.

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

Астрономическая теория ледниковых периодов позволяет определить составляющую эволюции климата планеты, которая обусловлена изменением ее орбитального и вращательного движений. Для ее исследования мы численно решали задачу о взаимодействии одиннадцати тел Солнечной Системы (девять больших планет, Солнце и Луна ) на интервале времени в 100 млн. лет [5-7] , т.е. исследовали классическую ньютонову проблему 11-ти тел.

Согласно закону всемирного тяготения, тело с номером i притягивается телом с номером k силой

, (1)

где G - гравитационная постоянная,

- радиус-вектор до тела с массой m_i от тела с массой m_k.

Если количество тел равно n, то их воздействие на i-оe тело выражается суммарной силой

(2)

и, в соответствии с классическими законами динамики [8] , движение тел в инерциальной системе отсчета определяется системой дифференциальных уравнений порядка 6n

i =1,2,…,n, (3)

где - радиус-вектор тела относительно центра масс (например, относительно барицентра Солнечной Системы).

Для n=11 (девять больших планет, Солнце и Луна) уравнения (3) представляют собой систему дифференциальных уравнений, описывающих всевозможные движения в выбранной модели Солнечной Системы. Для нахождения какого-либо ее частного решения мы, очевидно, должны задать 3n значений координат и 3n значений компонент скоростей на определенную дату, которую в дальнейшем будем называть начальной эпохой с T₀ = 0. Мы использовали различные варианты методов численного интегрирования, при этом в качестве системы координат мы выбрали барицентрическую экваториальную систему координат, ось x которой направлена на точку весеннего равноденствия _о на эпоху 1950.0 г. с номером юлианского дня JD_S = 2433282.4234.

При рассмотрении взаимодействий тел Солнечной Системы часто, кроме гравитационных сил для материальных точек, учитывается множество «слабых» дополнительных физических и динамических факторов (например, воздействие астероидов и спутников планет, отличие формы планет от однородных шаров, очевидно равнозначных материальным точкам, приливные силы, релятивистские эффекты, изменение массы тел со временем и т.п.). К примеру, в работах [3, 4] получена оценка влияния некоторых из указанных факторов и показано, что ими в принципе можно пренебречь. Константы взаимодействия многих дополнительных факторов определяются, как правило, по величине невязки между результатами, обусловленными ньютоновым взаимодействием и наблюдениями. Из этого, очевидно, вытекает, что в начале мы должны изучить, как можно точнее, эволюцию орбит планет в результате их гравитационного взаимодействия (1) и только после этого следует учесть воздействие любого из перечисленных выше более слабых факторов.

1.3 Начальные условия

Чтобы проще реализовать сравнительный анализ наших результатов с результатами других исследователей, мы выбрали в качестве начального момента дату 30.12.1949 г. с номером юлианского дня JD₀ = 2433280.5.

При этом мы сделали вычислительный эксперимент с использованием двух вариантов начальных условий. В основу первого варианта начальных условий были положены эфемериды DE 19 Лаборатории реактивного движения США (сокращенно - ЛРД, в английской транскрипции -JPL), взятые из Справочного руководства [8]. Для определения масс тел мы выполнили отдельное исследование имеющихся современных сведений о массах планет. К сожалению, численные данные о массах планет в специальной литературе сильно отличаются между собой. Например, для Плутона эта разница достигает порядка величины и более. В результате такого анализа мы определили значения масс планет и начальные значения их координат (они приведены в Таблице 1 Приложения 2). В статье Т.Куинна и др. [3] были использованы значения масс и начальные условия на ту же эпоху, но определенные по эфемериде DE102. Сравнивая наши и эти начальные условия, мы установили, что наибольшее отличие в массах имеется для массы Плутона, которое составляет почти 40%, в то время как массы Меркурия, Урана и Нептуна отличаются на десятые доли процента, а масса Солнца еще меньше - на 0.01%. Относительная разница в координатах положений планет и их скоростей составляет величину, меньшую 5·10^-5, что дает нам право считать, что сравнение полученных нами результатов с результатами других авторов более или менее корректно.

Для получения численных результатов, характеризующих динамическую эволюцию Солнечной Системы на промежутке времени в 100 млн. лет, сначала мы использовали первый вариант начальных условий и для проверки метода, сравнили полученные численные результаты с современными данными наблюдений. Для того, чтобы убедиться в достоверности наших результатов, мы сочли необходимым использовать также второй вариант начальных условий, при котором координаты и скорости тел были определены на ту же начальную эпоху -30.12.1949.0 г., но в другой системе отсчета (а именно, 2000.0 г. c JD_S = 2451544 и начальные условия для второго варианта были определены по известной JPL-теории DE 406/LE406, сведения о которой и программы для расчета представлены на сайте JPL Solar System Dynamics: http://ssd.jpl.nasa.gov/). Относительные значения масс были взяты из системы DE405, приведенные в работе Cтэндиша Е.М. [9]. После этого мы модифицировали их значения с учетом величины GM_E и с использованием констант известной системы IERS (см. стр. 426 Труды ИПА РАН [10]). Эти исходные данные и начальные условия представлены в табл. 2 Приложения 2.

1.4 Метод решения

Несмотря на то, что в математической литературе существует большое разнообразие численных методов интегрирования дифференциальных уравнений (см., например, монографию Т.В.Бордовицыной [11]), мы пришли к выводу, что конечно-разностные методы интегрирования для нашего исследования не обеспечивают необходимую точность. Это обусловило необходимость разработки специального алгоритма и специальной программы, названной нами «Galactica», суть которого состоит в том, что значение вычисляемой функции в «следующий» момент времени t=t₀ + t вычисляется более точно чем в «обычных» схемах численного интегрирования, а именно, мы определяли значение координаты x с помощью ряда Тейлора

, (4)

где x₀^(k) - производная порядка k в момент t₀, а целое число

K является «порядком» наивысшей производной.

Значение скорости (т.е. первая производная координаты x) определяется по аналогичной формуле, а ускорение - по формуле (3). Производные более высокого порядка x₀^(k) мы определяли путем непосредственного дифференцирования правых частей уравнений (3) в аналитическом виде. Отметим, что в небесной механике методы, основанные на разложении Тейлора, применялись и ранее (см., например, [12, 13]).

В расчетах мы использовали числа с двойной точностью, при которой они выражаются 17-ю десятичными знаками. Вычисления показали, что в выражении (4) можно ограничиться производными шестого порядка, т.е. мы брали K=6. Вычисления на большие промежутки времени были реализованы на суперкомпьютерах RM-600 и МВС-1000 в Новосибирском ВЦ СО РАН. В программе предусмотрена возможность записи результатов «в файл» через определенное количество шагов, например, через каждые 10 тысяч лет. Для анализа параметров орбиты, по этим данным просчитывался один оборот планеты вокруг Солнца (см. [14]) и по результатам расчетов определялись угол наклона орбиты i, долгота восходящего узла ц_Щ,, эксцентриситет e и долгота перигелия ц_р. (см. рис. 1). Кроме того, вычислялись и другие параметры, в том числе: радиусы, долготы и моменты прохождения через перигелий и афелий, средний за один оборот момент количества движения и наибольшее отклонение точек орбиты от её средней плоскости.

Время счета положений небесных тел на интервале времени в 10000 лет с шагом dt = 1·10^-4 года на суперкомпьютере МВС-1000 с процессорами DEC Alpha и частотой 833 МГц составляет 65 минут, т.е. время счета одного шага интегрирования равно t_1st = 3.9·10^-5 сек, а решение задачи на интервале времени, равном 100 млн. лет, заняло 2 года. Для сравнения приведем время счета одного шага интегрирования в работе Ж.Ляскара [2] на суперкомпьютере Сompaq alpha workstation с характеристиками (частота 833 МГц и t_1st = 8.64·10^-5 сек) и в работе Т. Куинна с соавторами [3] на суперкомпьютере Silicon Graphics 4D-25 workstation (t_1st = 3.74·10^-3 сек).

1.5 Некоторые соображения о достоверности результатов

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

Эффективным методом тестирования можно считать численное интегрирование дифференциальных уравнений осесимметричной задачи n-тел, точное аналитическое решение которой изложено в работах [5, 15, 16]. Некоторые результаты и доказательства достоверности примененного нами метода описаны в наших работах [7, 14]. При шаге интегрирования dt = 1?10^-4 года относительная погрешность момента количества движения Солнечной Системы за 100 млн. лет составила M = 8·10^-11, что на много меньше погрешности момента количества движения Солнечной Системы, приведенной, приведенной в работе [3] и составившей в которой на интервале времени интегрирования в 3.056 млн. лет величину дM = 1.55·10^-7. За этот же интервал времени погрешность наших расчетов для момента количества движения составила еще меньшую величину, равную дM = 3.93·10^-12. Ниже мы покажем, что величину погрешности момента количества движения можно уменьшить даже до величины дM = 1.47·10^-15 за 100 млн. лет.

Рис. 1. Параметры орбиты Марса в неподвижной экваториальной x_,yz и подвижной эклиптической x_ey_ez_e гелиоцентрических системах координат.

1 - небесная сфера; 2 - плоскость экватора Земли на 1950 г.; 3 - плоскость орбиты Земли на 1950 г. (плоскость эклиптики); 4 - плоскость орбиты Марса в в будущую эпоху Т; 5 - плоскость экватора Земли в эпоху Т; 6 - плоскость орбиты Земли в эпоху Т (наклон для наглядности увеличен) ; N - северный полюс мира; - северный полюс подвижной эклиптики; - точка весеннего равноденствия 1950 г.; - точка на линии пересечения подвижного экватора в эпоху Т с подвижной эклиптикой (точка весеннего равноденствия в эпоху Т); G - дуга большого круга, перпендикулярного плоскости орбиты Марса; B - гелиоцентрическая проекция перигелия Марса на небесной сфере; А - восходящий узел орбиты Марса на подвижной эклиптике; D - восходящий узел орбиты Марса на неподвижном экваторе 1950 г.; параметры орбиты Марса в инерциальной системе: = ₀D; _р.= DB; i = ₀DG; и в подвижной эклиптической системе: Щ_a = A; щ_a = AB; р_a = A + AB = Щ_a + щ_a; i_еа = i_е = AG; индекс «_a» - по данным наблюдений.