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

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

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

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

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

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

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

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

Москва-Тюмень, 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» - по данным наблюдений.

2. Сравнение результатов вычислений с аппроксимацией данных наблюдений по Ньюкомбу

2.1 Приведение вычисленной долготы перигелия к неподвижной точке на орбите

При сравнении результатов интегрирования уравнений (3) с наблюдениями возникает ряд трудностей, особенно если мы увеличиваем точность сравнения. Можно, например, сравнивать координаты тел, периоды их обращения вокруг Солнца, величины вековых возмущений. Сравнение положений тел, полученных в результате интегрирования в интервале времени, равном 50-ти годам, с эфемеридами DE406/LE406, указывает на их хорошее совпадение (см. табл. 4 из Приложения 3В). Под вековыми возмущениями мы понимаем законы изменения, прежде всего, угловых параметров орбит на интервале времени порядка одного тысячелетия. Так как вековые возмущения С. Ньюкомба, представленные им в виде «законов изменения» средних элементов эллиптической орбиты, определены на основании аппроксимации известных в астрономии наблюдений, то сравнение решений уравнений (3) с последними, позволяет убедиться в том, что они, вообще говоря, пригодны и на большем интервале времени.

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

На рис. 1 показаны основные параметры орбиты Марса (i, _р и ), которые определяются в результате интегрирования уравнений (3). Так как долгота перигелия _р = DB отсчитывается от подвижного узла D, то при изменении углов, определяющих положение плоскости орбиты (то есть величин i и ), точка D будет перемещаться вдоль линии DAB и тем самым будет «вносить вклад» в величину _р. К сожалению, невозможно зафиксировать некоторую неподвижную точку на орбите, отсчет от которой позволил бы определить чистое перемещение перигелия, не искаженное колебаниями самой орбиты. Только при совершении перигелием нескольких оборотов можно определить, более или менее точно, среднюю скорость его перемещения. Чтобы точнее «зафиксировать» мгновенное перемещение перигелия, можно, например, выбрать «условно неподвижную» точку на орбите и пересчитать положение перигелия по отношению к ней. С этой целью, через восходящий узел орбиты Земли ₀ на неподвижном экваторе (см. рис. 1), проведем большой круг, перпендикулярный к орбите Марса. Обозначим точку их пересечения буквой G. Если положение плоскости орбиты планеты относительно вектора (момента количества движения всей Солнечной Системы) определять с помощью углов Эйлера (угол прецессии -? и угол нутации- ?), тогда эволюция плоскости орбиты планеты будет определяться прецессионным движением со скоростью вокруг вектора и колебаниями угла нутации ? [17, 18]. Так как вектор близок к перпендикуляру, проведенному к плоскости орбиты Марса, то, при колебаниях орбиты GDAB, точка G будет иметь пренебрежимо малые смещения по орбите, что дает право отсчитывать от нее угловое расстояние перигелия _р0. Величина _р0 состоит из двух дуг

_p0 = GD + DB = GD + _р. (7)

Найдем из прямоугольного треугольника ₀GD дугу GD, если учесть, что дуга ₀D =,, угол GD = /2, а угол D = i. Из сферической астрономии (см., например, Справочное руководство, [8]) имеем, что

sin₀G = sin₀D·sini, (8)

cosGD = cos₀D/cos₀G. (9)

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

. (10)

2.2 Приведение данных наблюдений к неподвижной системе координат

2.2.1 Параметры орбиты Марса в подвижной системе координат по данным наблюдений

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

В работах Ньюкомба [21], опубликованных в 1895-1898 годах, приведены приближенные формулы для осредненных параметров орбиты Марса (эксцентриситета, долготы перигелия, долготы узла и наклона) в подвижной эклиптической системе координат, которые он получил на основе анализа и аппроксимации, имеющихся в его распоряжении, наблюдений:

e_a = 0.09331290 + 0.000092064·T_j - 0.000000077·T_j ²; (11)

р_a = 334^_1305.53 + 6626.73·T_j + 0.4675·T_j ² - 0.0043·T_j ³; (12)

Щ_a = 48^_4711.19 + 2775.57·T_j - 0.005·T_j ² - 00192·T_j ³; (13)

i_еа = 1^_5101.20 - 2.430·T_j + 0.0454·T_j ²; (14)

где T_j -время, отсчитываемое в юлианских столетиях по 36525 «эфемеридных» суток от фундаментальной эпохи 1900, январь, 12^h UT с номером юлианского дня JD_b = 2415020.3134. Индекс a означает, что коэффициенты в формулах получены в результате обработки наблюдений.

Рис. 2. Геометрические схемы: а - изображение в плоскости подвижной эклиптики; б - на небесной сфере. Другие обозначения даны на рис. 1. Символ A_MrA_Mr' - плоскость экватора Марса в текущую эпоху.

Формулы (11) - (14) заимствованы из Справочного руководства [8], в которых учтены некоторых незначительных уточнений Росса. Долгота Щ_a (см. рис. 1) также отсчитывается в плоскости подвижной орбиты Земли от ее подвижного восходящего узла г. Отметим, что р_a = Щ_a + щ_a является некоторой «условной» долготой перигелия, так как величины , определяются в разных плоскостях (щ_a - угловое расстояние точки перигелия B в плоскости орбиты Марса от его восходящего узла A, а величина Щ_a вычисляется в другой плоскости.

2.2.2 Параметры орбиты Земли в неподвижной системе координат

Так как элементы орбиты Марса отнесены к подвижной орбите Земли, сначала найдем их значения в неподвижной системе координат. Известно [8], что плоскость эклиптики «поворачивается» в своей плоскости (см. рис. 2 a) относительно мгновенной оси вращения, расположенной под углом П к оси x_e,,

П = 1735703 + 3287·T_t + 0,6·T_t², (15)

с годичной скоростью

= 0,47107 - 0,00068·T_t, (16)

где T_t - время в тропических столетиях по 36524.22 «эфемеридных» суток от фундаментальной эпохи с номером юлианского дня JD_b.

Если мы фиксируем ось x_e неподвижной системы координат в эпоху с номером юлианского дня JD_S, то проекция угловой скорости ·cosП на эту ось будет означать угловую скорость поворота эклиптики в неподвижной системе координат. Так как в неподвижной системе координат положение экватора А₀А₀ (см. рис. 2б), к которому подвижная эклиптика наклонена под углом i_Ea, известно, отсюда вытекает, что скорость изменения этого угла будет равна

·cos . (17)

Обозначим угол между эклиптикой и неподвижным экватором в неподвижной системе координат в эпоху JD через = е_a(T_S). Эта величина является начальным условием при интегрировании уравнения (17). Учтем при этом, что угол наклона подвижной эклиптики к подвижному экватору Земли определяется выражением

е_a(T_j) = 23^_2708,26 - 46,845·T_j - 0,0059·T_j² + 0,00181·T_j³,

а величина T_S = (JD_S - JD_b)/36525 представляет собой промежуток времени в юлианских столетиях от фундаментальной эпохи JD_b до эпохи JD_S.

Так как равенства (15) и (16) содержат время T, то оно входит в оба сомножителя соотношения (17). Используя среднее значение угла П для среднего момента (T_S + T_j)/2, проинтегрируем уравнение (17), в результате чего получим:

i_Ea = ?_S + 100[0,47107·(T_t - T_St) - 0,00034·(T_t - T_St)²]*

*cos1735703 + 3287·0.5·(T_St + T_t) + 0,6·0.25·(T_St + T_t)², (18)

где ?_S - наклон эклиптики к экватору в эпоху T_S, а величины

T_St = T_S·k_j, k_j = 36525/36524.22.

Формула (18) определяет «закон» изменения наклона эклиптики по отношению к неподвижной, «фиксированной» на эпоху T_S, плоскости экватора плоскости экватора, согласно данным наблюдения.

В результате решения уравнений (3) мы получаем переменные значения параметров орбиты Земли и, в частности, закон изменения угла наклона i_E плоскости её орбиты к неподвижному экватору и движение восходящего узла _E по неподвижному экватору. Законы изменения наклонов i_E и i_Ea с высокой точностью совпадают, так как их относительная разность на интервале времени 1000 лет не превышает величины .

Вывод зависимости восходящего узла эклиптики от вековых возмущений Ньюкомба относительно неподвижных координат затруднительно, поэтому целесообразно использовать полученное ранее в результате решения уравнений (3) выражение для перемещения восходящего узла орбиты Земли по неподвижному экватору _о2 (см. рис. 2 б) в виде среднеквадратичной аппроксимации по параболе:

_о2 = aT_n² + b T_n + c при -25 T_n 25. (19)

Здесь T_n - используемое при численном интегрировании время в юлианских столетиях от начальной эпохи JD₀,

a = 2.397234912503438·10^-6; b = 5.535655301401231·10^-5; c = -2.089118512971404·10^-6.

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

a = 2.376361216684775·10^-6; b = 4.788549396485488·10^-5; c = -2.325022430798324·10^-5.

Следует отметить, что перемещение восходящего узла орбиты Земли по неподвижному экватору _о2a мы также получили, используя средние элементы орбиты Земли, приведенные к системе 2000.0 г., в теории движения планет, разработанной во Франции

Симоном Дж.Л. и др. [20]. Отметим при этом, что величины дуг _о2 и _о2a совпадают с высокой точностью в любой из моментов рассмотренного диапазона времени -25 T_n 25.

В результате прецессии земной оси восходящий узел (см. рис. 2 б) перемещается против движения Земли, по подвижной эклиптике EE', с угловой скоростью, отнесенной к одному году (см. например, формулу из [ 8 ]),

p = 50,25641 + 0,02223 T_t. (20)

Тогда в любой момент времени Т смещение точки весеннего равноденствия по подвижной эклиптике относительно момента Т_S будет определяться дугой ₂ (см. рис. 2 б), полученной в результате интегрирования равенства (20), в следующем виде:

₂ = 5025.641·(T_t - T_St) + 2.223·(T_t - T_St)²/2. (21)

Для сравнения с результатами интегрирования уравнений (3) необходимо «время» в тропических столетиях T_t в зависимостях (18) и (21) заменить формулой: T_t = k_j·[T_n + (JD₀ - JD_b)/36525]. Таким образом, наклон подвижной эклиптики к неподвижному экватору определяется выражением (18), перемещение ее восходящего узла - формулой (19), а положение подвижного экватора - равенством (21).

2.2.3 Параметры орбиты Марса в неподвижной системе координат

Долготы _а и _а восходящего узла (точки D на рис. 2 б) и перигелия (точка B) Марса, соответственно, отсчитываются от точки , поэтому на их величины влияет прецессия ₂ восходящего узла Земли. Аналогично, угол наклона плоскости орбиты Марса i_еа зависит от изменения угла наклона орбиты Земли i_Ea, поэтому целесообразно вычислить вековые возмущения орбиты Марса относительно неподвижной плоскости экватора 1950.0 г. A₀A₀' (см. рис. 2 б) как функции от следующих параметров: 1) углового расстояния восходящего узла ^a=г₀D, отсчитываемого в плоскости экватора 1950 г. от восходящего узла Земли ₀; 2) _p0^a = GB - углового положения перигелия, отсчитываемого в плоскости орбиты Марса от неподвижной точки G; 3) угла наклона i^a = DG плоскости орбиты к неподвижной плоскости экватора.

Из рис. 2 б видно, что в неподвижных координатах угловое положение перигелия в плоскости его орбиты определяется дугами GD, DA и AB:

_р0^a = GD +DA+ AB. (22)

Для определения этих дуг рассмотрим сферический треугольник ₂DА, в котором известны два угла: ₂ = i_Ea, А = i_еа, и сторона ₂А = _а - ₂, лежащая против неизвестного угла D. Его можно определить по «теореме косинусов» (формула (1.1.010) из [8]):

cosD = -cos i_Ea·cosi_eа + sin i_Ea·sini_eа·cos(_а - ₂). (23)

Так как угол наклона орбиты к плоскости неподвижного экватора i^а = - D, то тогда из (23) получаем его значение в виде:

i^a = - arccos(-cos i_Ea·cosi_еа + sin i_Ea·sini_eа·cos(_а - ₂)). (24)

Далее, из треугольника ₂DА, по теореме синусов (см. формулу (1.1.007),

[8]) находим величину дуги ₂D и положение восходящего узла в виде формул:

= ₀D =_{02 + 2}D = ₀₂ + arcsin[sin(_а - ₂) sin i_Ea /sini^a]. (25)

Для определения положения перигелия на орбите относительно неподвижной точки G (см. рис. 2 б), согласно (22), найдём аналогично дугу DA из ₂DA :

sin i_Ea / sin DA = sin ( - i^а) / sin ₂A.

Из этого вытекает, что дуга

DA = arcsin[sin i_Ea·sin(_а - ₂) / sini^а]. (26)

Дуга AB, фигурирующая в выражении (22), определяется по формуле

AB = щ = _а - _а. (27)

Из прямоугольного сферического треугольника ₀GD, в котором известны дуга ₀D = и угол D = i^а, определяем дугу ₀G согласно формуле

sin _оG = sin _оD sin i^a; cos GD = cos _оD/ cos _оG,

и далее, определяем дугу

GD = arcos{cos ^a/[1 - (sin ·sin i^a)2]^0.5}. (28)

После подстановки значений дуг DA, AB, и GD в соотношение (22), получаем угловое расстояние перигелия Марса в неподвижных координатах, которое отсчитывается от неподвижной точки G на его орбите, в следующем виде:

_p0 = _a - _а + arcsin[sini_Ea sin(_а - ₂)/sin i^а] + arccos[cos /(1 - (sin sin i^а)²)^0.5]. (29)

Таким образом, используя наблюдения, можно вычислить наклон орбиты Марса к неподвижной плоскости экватора i^а, долготу восходящего узла ^a с помощью соотношений (24) и (25), соответственно, и угловое расстояние перигелия _p0 от неподвижной точки G на орбите Марса с помощью соотношения (29).

2.3 Сравнение вычисленных вековых возмущений элементов орбиты Марса со значениями, полученными из наблюдений

Для этого мы проинтегрировали дифференциальные уравнения (3) с использованием, приведенных выше, двух вариантов начальных условий на интервалах времени: -2.6 тыс. лет T 2 тыс. лет для первого варианта и -3.4 тыс. лет T 3.6 тыс. лет - для второго. Мы решали эти задачи с шагом интегрирования dt = 1?10^-5 года и с «расширенной длиной» чисел ( с четырех кратной точностью), при которой вычисляемые координаты представляются числами с 35-ю десятичными знаками. Вычисления показали, что погрешность в этом случае пропорциональна [8] интервалу времени интегрирования и изменяется со скоростью, примерно равной dM/dt = 1.47·10^-21 в столетие. Вычисления с заданными начальными условиями и параметрами на интервале времени в 100 млн. лет, привели бы к относительной погрешности, примерно равной, M = 1.47·10^-15, что значительно меньше, приведенной нами ранее, величины M = 8 10^-11. Не смотря на то, что погрешность уменьшилась в 54000 раз, результаты вычислений, по сравнению с результатами интегрирования «с двойной точностью» и с шагом интегрирования dt = 1?10^-4 года, практически совпали. Это соображение говорит в пользу достоверности проведенных вычислений.

На рис. 3 точками представлена динамика параметров орбиты Марса: эксцентриситета e, долготы восходящего узла , угла наклона плоскости орбиты i и угла положения перигелия _р0 в системе неподвижного экватора, вычисленная методам численного интегрирования системы (3), описывающей взаимодействия одиннадцати тел Солнечной Системы. Непрерывные линии представляют те же величины (e_a, i^а, ^а и _р0^а), полученные из наблюдений в промежутке времени, порядка 2500 лет, после некоторой математической обработки. По этой причине их достоверность ограничена этим периодом времени. Из графиков видно, что результаты наблюдений хорошо согласуются с вычисленными из дифференциальных уравнений значениями элементов e, i, и _р в пределах нескольких тысяч лет «в прошлое» и «в будущее» от исходной эпохи. Наибольшее отклонение имеет наклон i при Т = -2.6 тыс. лет с относительной погрешностью (i -i^а)/ i^а = 7.4·10^-4.

Рис. 3. Сопоставление рассчитанной (1) эволюции орбиты Марса за -2.6 ч +2 тыс. лет с аппроксимацией С. Ньюкомба данных наблюдений (2): e - эксцентриситет; i -наклон плоскости орбиты к плоскости экватора 1950.0 г. в радианах; - угловое положение восходящего узла орбиты от оси x на эпоху 1950.0 г. в радианах; угловое положение перигелия в плоскости орбиты: _p0 - от неподвижной точки G на орбите (1, 2) и _p - от восходящего узла в радианах (3). Т - время в тысячах юлианских лет от 30.12.1949 г.; интервал между точками - 200 лет.

Близость результатов «теоретических» вычислений с данными наблюдений, представленных «эмпирическими» соотношениями, свидетельствует в пользу достоверности наших расчетов. Анализ также показывает, что при изучении эволюции орбиты Марса в течение 300 непрерывных оборотов, появляются небольшие короткопериодические колебания в элементах орбиты и, прежде всего, в эксцентриситете орбиты Марса, период которых равен двенадцати годам. Эти колебания показаны на рис. 3.

На рис. 3 представлены результаты вычислений с использованием первого варианта начальных условий (согласно табл.1), отнесенных к эпохе 1950.0 г. Такие же вычисления были нами выполнены и со вторым вариантом начальных условий (согласно табл.2), но отнесенных к эпохе 2000.0 г. Их результаты еще ближе к «эмпирическим» соотношениям. Эти сравнения убеждают нас в справедливости, приведенных выше, количественных и качественных выводов.

3. Эволюция орбиты Марса в интервале времени 3 млн. лет

Эволюция вычисленных элементов орбиты Марса за 3 млн. лет «в прошлое» представлена на рис. 4. Она показывает, что эксцентриситет орбиты е испытывает короткопериодические изменения с амплитудой 0.019, с «главным» периодом, равным Т_е1 = 95.2 тыс. лет при среднем (за 50 млн. лет) значении е_m = 0.066. Долгота восходящего узла ц_Щ меняется с «средним» периодом Т_Щ = 73.1 тыс. лет вокруг среднего значения ц_Щт = 0.068 радиан. Угол наклона плоскости орбиты і к плоскости экватора (эпоха 1950 г.) испытывает колебания с таким же периодом Т_і = 73.1 тыс. лет вокруг среднего значения і_m = 0.405 радиан. Долгота перигелия ц_р практически линейно увеличивается со временем, т.е. перигелий перемещается в направлении обращения Марса вокруг Солнца, совершая в среднем один оборот за T_p = 76.8 тыс. лет. Угловая скорость ?_p вращения перигелия колеблется вокруг среднего значения ?_pm = 1687 " за столетие, при этом в момент времени t = -1.35 млн. лет она принимает отрицательное значение.

Рис. 4. Эволюция орбиты Марса за 3 млн. лет: обозначения те же, что и на рис. 3; Т - время в миллионах лет; интервал между соседними точками равен 10 тыс. лет; ?_p - угловая скорость вращения перигелия в "/столетие на интервале времени в 20 тыс. лет; средняя за 50 млн. лет угловая скорость вращения перигелия ?_pm = 1687 "/столетие - средняя за 50 млн. лет.

4. Эволюция орбиты Марса за 100 млн. лет.

На рис. 5а представлено изменение параметров орбиты Марса за 50 млн. лет. На нем отчетливо виден второй период изменения эксцентриситета, равный Т_е2 = 2.31 млн. лет.

Рис. 5. Эволюция орбиты Марса за 100 млн. лет (обозначения те же, что и на рис. 3 - 4).

Размещено на http://www.allbest.ru/

При этих колебаниях достигаются «крайние» значения эксцентриситета е = 0.0014 и е = 0.126. Колебания угла наклона i происходят в пределах от 0.288 до 0.521 радиан, их диапазон составляет 13.37є. Долгопериодические колебания восходящего узла ц_Щ и наклона орбиты i обусловлены колебаниями ее угла нутации ? со средним периодом, равным Т_? = 1.15 млн. лет. Следует отметить, что период Т_? в два раза меньше Т_е2. Перигелий со средним за 50 млн. лет периодом T_p перемещается в направлении обращения Марса вокруг Солнца. Угловая скорость ?_p вращения перигелия колеблется вокруг среднего значения ?_pm. При сопоставлении с графиком для e видно, что возвратные движения перигелия происходят, когда эксцентриситет орбиты приближается к нулю. Из представленных графиков следует, что амплитуды колебаний и периоды стабильны, т.е. орбита Марса совершает «стационарные» движения.

На рис. 5б показано изменение параметров орбиты Марса в диапазоне времени -50 млн. лет T -100 млн. лет. Из сравнения графиков, на рисунках 5а и 5б, видно, что характер эволюции параметров орбиты в течение двух периодов времени по 50 миллионов лет каждый сохраняется. Периоды колебаний эксцентриситета е, восходящего узла орбиты и угла наклона её плоскости i к плоскости неподвижного экватора, а также их амплитуды и период Т_p обращения перигелия остаются теми же. Аналогичные результаты получены для орбит других планет на том же промежутке времени в 100 млн. лет. Это свидетельствует об устойчивом характере планетных движений в Солнечной Системе.

5. Сравнение наших результатов с результатами других исследователей

В 1950 г. Д. Брауэр и А. Вурком опубликовали новую аналитическую теорию вековых возмущений элементов орбит планет Солнечной Системы [21]. На ее основе К. Остервинтер и др. [22] вычислили изменения элементов орбит больших планет, от Меркурия до Нептуна, за период времени -5 млн. лет T 5 млн. лет. В этой работе приведены, в частности, углы наклона плоскости орбиты к неподвижной эклиптике _1оА₁ (см. рис. 2б). Для сравнения с ними мы определили значения наклонов планетных орбит к плоскости неподвижной эклиптики по значениям рассчитанных нами наклонов к плоскости неподвижного экватора. С этой целью мы рассмотрели треугольник ₀DA₁, в котором известны два угла = ₀ и DA₁ = - і, сторона между ними ₀D = _. И на основе сферической тригонометрии вычислили другие его параметры Угол наклона орбиты і_е0 = A₁D к неподвижной эклиптике может быть вычислен по «теореме косинусов» [8] и будет равен

і_е0 = arccos[-cos₀ cos( - і) + sin₀ sin( - і) cos_.]. (30)

Из этого сравнения видно, что значения эксцентриситетов е в начале (при T = 0) совпадают количественно, затем («после» -1 млн. лет) начинают различаться качественно из-за того, что второй период изменения эксцентриситета Т_е2 = 2 млн. лет меньше найденного нами периода (Т_е2 = 2.31 млн. лет), который, как будет показано на рис. 7, совпадает с периодом колебаний Ж.Ляскара и др. [4]. Это говорит о том, что в работе [22]

Рис. 6. Сравнение численных решений (1) эксцентриситета и наклона орбиты Марса в промежутке времени, равном 5-ти миллионам лет (от -5 млн. лет до «сегодняшнего дня», Т=0) с расчетами (2) К. Остервинтера и др. [22], полученных ими на основании теории Брауэра и Вуркома [21].

долгопериодические изменения эксцентриситета вычислены, по-видимому, с большой погрешностью. Что касается углов наклона, наши результаты и результаты, приведенные в работе [22], фактически совпадают. В статье [25] период больших колебаний наклона i_e0 равен, Т_ie2 = 1.166 млн. лет, т.е. близок к «нашему периоду колебаний», Т_? = 1.15 млн. лет, для угла нутации ?. Однако периоды коротких колебаний, приведенные в статье [22], в два с лишним раза больше вычисленных нами периодов колебаний и периодов колебаний, приведенных в работе [4] (см. рис. 7). Короткопериодические колебания углов наклона, приведенные в работе [22], по-видимому, вычислены с большей погрешностью. Однако, в отличие от работы [4], эволюция колебаний плоскости орбиты Марса, описанная в статье [22] так же, как и у нас, имеет устойчивый характер.

В работе [4] долгота перигелия _L отсчитывается от подвижного восходящего узла (moving equinox), а угол наклона I_Mr плоскости орбиты (obliquity) отнесен к подвижной плоскости экватора Марса. Для дальнейшего необходимо привести эти результаты к неподвижному земному экватору и с этой целью выполним некоторые геометрические построения.

Восходящий узел находится на пересечении кругов обиты Марса DАВ и его экватора А_Mr А?_Mr (см. рис. 2б). Восходящий узел _Mra расположен на другом конце диаметра от нисходящего узла _Mrd, показанного на этом рисунке. Поэтому долгота перигелия _L определяется дугой _MraD_MrdB, которую можно записать в виде суммы дуг:

...