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

_L = _MraD_Mrd+_MrdB = +_MrdB = +_pL -D_Mrd, (31)

где ц_pL = DB - долгота перигелия от «опорного» земного экватора А₀А?₀, поэтому долгота перигелия от неподвижного земного экватора запишется в виде [4]:

ц_pL = _L + D_Mrd - + 2, (32)

где величина 2 добавлена для преобразования отрицательного угла долготы перигелия в начальную эпоху (так как _L(0) = 1.24, D_Mrd(0) = 1.42) в положительный угол.

Дуга D_Mrd изменяется за счет прецессии марсианского экватора, принятой в статье [4] равной p_Mr = 7.597?/год. Прецессия происходит «за стрелкой часов», т.е. точка _Mrd приближается к точке D со скоростью p_Mr. Если обозначить положение нисходящего узла в начальную эпоху через _0Mrd, тогда смещение _Mrd по орбите Марса в радианах запишется в следующем виде:

_0MrdMrd = -100·2··p_Mr·T_j/(3600·360), (33)

где T_j - время в юлианских столетиях от начальной эпохи. Тогда дуга D_Mrd в любой момент времени выражается в виде суммы двух дуг:

D_Mrd = D_0Mrd + _0MrdMrd, (34)

где D_0Mrd - величина этой дуги в начальную эпоху.

Для определения дуги D_0Mrd воспользуемся средними экваториальными (земными) координатами северного полюса N_Mr Марса, взятыми из Справочника [8], стр. 65]:

б_NMr = 317°.32 - 0°.1011T_t; ?_NMr = 52°.68 - 0°.0570T_t (35)

и угол наклона экватора Марса к плоскости его орбиты

I_Mr^a = 25°.19969 + 0°.01219T_j + 0°.00006T_j², (36)

где T_t и T_j время в тропических и юлианских столетиях от начальной эпохи 1950.0

В треугольнике DK_Mrd с помощью параметров (35) - (36) мы определяем два угла K = д_NMr + /2 и _Mrd =I_Mr^a, а также величину

DK = б_NMr + /2 - ц_?^а, (37)

где величина ц_?^а определяется из выражения (25).

Используя основные теоремы сферической астрономии (например, соотношение sinD_Mrd^a / sinK = sinDK /sin_Mrd), находим величину дуги

D_Mrd ^a = arcsin[sin(/2 + д_NMr)·sin(б_NMr-ц_Щ^a + /2) / sin(I_Mr)]. (38)

Индекс «а» означает, что дуга определена по данным наблюдений. После подстановки (34) в (32) угловое расстояние перигелия Марса от неподвижного земного экватора в работе [4] запишется в виде:

_рL = _L + D_0Mrd + _{0Mrd Mrd} + , (39)

где дуга D_0Mrd определяется из выражения (38) в начальную эпоху, т.е. при T_t = T_j = 0.

С помощью угла I_Mr = D_MrdK из работы [4], определим величину угла i_L = _MrdDK (угол между орбитой Марса и неподвижной плоскостью земного экватора А₀А?₀ (см. рис. 2б)). В треугольнике DK_Mrd известны два угла: угол I_Mr и из наблюдений угол K = -/2+д_NMr. Из наблюдений известны также две стороны: D_Mrd и DK, соответственно (см. соотношения (34) и (37)). Неизвестный угол i_L находится против неизвестной стороны K_Mrd. Применяя «теорему косинусов» [8] к величине угла i_L и стороны K_Mrd, после «исключения стороны» из двух равенств получаем формулу:

i_L = arcos{[sinI_Mr·cos(дN_Mr)·cos(DK)·cos(D_Mrd^a)+cosI_Mr·sin(дN_Mr)]/[1-sinI_Mr·cos(дN_Mr)·sin(DK)·sin(D_Mrd^a)]}. (40)

Входящие в выражение (40) приближенные представления для наблюдений (соотношения (13), (14), (24), (25), (35) и др.) справедливы на интервале времени порядка 1000 лет, поэтому формула (40) пригодна тоже на таком интервале.

Результаты, полученные нами и относящиеся к эволюции орбиты Марса на промежутке времени в 3 млн. лет, мы сравнили (см. рис. 7) с результатами, приведенными в работе [4]. Величины эксцентриситетов е совпадают с высокой точностью практически на всем этом интервале времени, поэтому можно, с большой вероятностью, допустить, что в вычислениях, приведенных в работе [22] и касающихся определения эксцентриситета, по-видимому, имеются погрешности.

Рис. 7. Сравнение численных решений (1) эволюции орбиты Марса на интервале времени (-3 млн. лет - 0) с расчетами (2) из работы [4] для эксцентриситета (e), угла наклона плоскости орбиты (i) к плоскости экватора (эпоха 1950 г.) и углового расстояния перигелия от неподвижного экватора (ц_р); I_Mr и _L - наклон орбиты и долгота перигелия в подвижных марсианских координатах ; I_Mr^a вычислялась по данным наблюдений с помощью (36); величина i_L - по данным Ж.Ляскара [4], вычисленным по формуле (40).

На рис. 7 мы сравнили значения угла наклона в двух интервалах времени: от 0 до - 100 тыс. лет и от 0 до -3 млн. лет. Из верхнего графика на рис. 7 видно, что на протяжении трех тысяч лет углы наклона орбиты Марса к подвижному экватору Марса совпадают с «наблюдаемыми» углами I_Mr^а. Это подтверждает справедливость выбора авторами статьи [4] начальных условий и корректность первых шагов интегрирования. Преобразованный угол i_L к плоскости неподвижного земного экватора также на протяжении 3-х тысяч лет совпадает с вычисленным нами углом i. Так как наши результаты совпадают с «наблюдениями» (см. рис. 3), это подтверждает правильность преобразования (40). Выражение (40) вне интервала аппроксимации не может быть использовано, поэтому на втором интервале времени мы сравнивали углы наклона i и I_Mr.

Вычисленные нами углы наклона плоскости орбиты i к неподвижной плоскости экватора качественно совпадают с углами, приведенными в статье [4] в «начальном» интервале времени (от 0 до 60 тыс. лет). Затем амплитуды углов I_Mr,, приведенные в [4], после 400 тыс. лет резко возрастают. Так как амплитуды колебаний наших углов наклона к плоскости орбиты i_e согласуются с амплитудами углов из работы К. Остервинтера и др. [22], (см. рис. 6), отсюда следует, что увеличение углов I_Mr, приведенное в работе [4], обусловлено изменением угла наклона плоскости экватора Марса. Мы отмечали, что вычисленный нами период колебаний Т_? = 1.15 млн. лет угла нутации ? плоскости орбиты Марса близок к периодам, приведенным в статье [4] и, можно сказать, что фактически совпадает (см. рис. 8) с периодами модуляции амплитуд углов I_Mr из работы [4]. Из этого можно заключить, что периоды изменения наклона орбиты Марса, приведенные в работе [4], вычислены корректно.

На рис. 7 сопоставлена динамика вычисленных нами угла перигелия ц_p и углов ц_pL из работы [4] по формуле (32). Так как у названных авторов, независимо от количества оборотов перигелия углы р_L изменяются от 0 до 2р, мы их привели к некоторому непрерывному углу. Видно, что эти углы практически совпадают на рассматриваемом интервале времени в 3 млн. лет. Совпадают также «возвратные» движения перигелия, которые имеют место при T = -1.35 млн. лет. На графиках также приведены значения долготы перигелия р_L относительно подвижного марсианского экватора.

Сравнение полученных нами результатов с результатами из работы [4], относящиеся к эволюции орбиты Марса за 21 млн. лет, приведено на рис. 8. Значения эксцентриситетов е качественно совпадают на всем этом интервале, т.е. совпадают периоды и амплитуды как короткопериодических, так и долгопериодических колебаний (см. рис. 7). Из этого вытекает, что теория, изложенная в статьях [1, 2, 4], позволяет с хорошей точностью вычислить эксцентриситет орбиты Марса. Амплитуда угла наклона I_Mr плоскости орбиты Марса к плоскости подвижного экватора Марса, начиная с -5 млн. лет, значительно возросла, т.е. произошло качественно новое изменение положения орбиты Марса относительно его подвижного экватора. Из графиков эволюции перигелиев Марса за период в 21 млн. лет (см. рис. 8, где представлены угловые положения перигелия ц_р и ц_pL) можно заключить, что наши численные исследования практически совпадают с результатами Ж.Ляскара и его соавторов [1,2,4 ].

Рис. 8. Сравнение численных решений (1) с вычислениями (2) из работы [4], показывающими эволюцию орбиты Марса за 21 млн. лет (-21 ч 0 млн. лет). Обозначения приведены на рис. 3.

6. Проблема устойчивости Солнечной Системы

На протяжении последних трех столетий различными аспектами проблемы устойчивости Солнечной Системы занимались выдающиеся ученые (И.Ньютон, Л.Эйлер, Ж.Лагранж, П.Лаплас, К.Гаусс, А.Пуанкаре, А.М.Ляпунов, О.Ю.Шмидт, А.Н.Колмогоров, Н.Д.Моисеев, Г.Н.Дубошин и другие), но, несмотря на их огромные достижения, эта фундаментальная проблема естествознания до конца до сих пор не решена.

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

Например, для интервала времени в 250 млн. лет Ж.Ляскар и его коллеги показали [4], что угол наклона орбиты к подвижной плоскости экватора I_Mr может изменяться по разному, в зависимости от принятой начальной скорости прецессии p_Mr и может достичь при этом больших величин. Значительное изменение угла наклона I_Mr плоскости орбиты Марса к плоскости подвижного экватора Марса (см. рис. 8) авторы интерпретируют как отсутствие устойчивости его орбиты и на этой основе делается вывод о неустойчивости Солнечной Системы и о хаотичности движений в ней.

Однако исследования для меньших периодов времени не подтверждают этот вывод. Например, эволюция наклона i_еL орбиты Марса к эклиптике 2000.0 г. (см. верхние четыре графика на рис. 9), в интервале времени в 40 миллионов лет, указывают на «устойчивость» его орбиты. На этих графиках результаты численного интегрирования (сплошная линия) были сравнены с их решением по вековым уравнениям (штриховая линия). Они отмечают, что заметное расхождение двух решений наблюдается после интервала времени в 32 млн. лет. Ниже, в таком же масштабе, выполнено сравнение угла наклона орбиты i_е0 к неподвижной эклиптике 1950г., вычисленного нами по формуле (30) по нашим решениям уравнений (3), с эволюцией угла i_еL, описанной в работе [4]. Их поведение весьма похоже. При наложении графиков i_е0 и i_еL, один на другой, они полностью совпадают на интервале времени в 4 миллиона лет, а существенные расхождения начинаются «после» -32 млн. лет. Это дает нам право сделать, по меньшей мере, два вывода:

Во-первых, совпадение углов i_е0 и i_еL подтверждает наши соображения в отношении решений К. Остервинтера с соавторами [22], представленных на рис. 6. Во-вторых, подтверждается наш вывод о том, что значительные изменения угла наклона I_Mr плоскости орбиты Марса к подвижной плоскости его экватора (см. работу [4]), вызваны движениями плоскости экватора.

Рис. 9. Сравнение эволюции угла наклона (i_eL) плоскости орбиты Марса к плоскости эклиптики 2000.0 г. по расчетам Ж. Ляскара. и др. [4], (четыре верхних графика), с эволюцией угла наклона (i_e0) плоскости орбиты Марса к плоскости эклиптики 1950.0 г. по нашим расчетам (четыре нижних графика) на интервале времени в 40 млн. лет.

Мы исследовали орбиты планет в неподвижных координатах. Из рис.8 видно, что эволюция вычисленного нами угла наклона i плоскости орбиты Марса неизменна в течение 21 млн. лет. Неизменна эволюция наклона i и других параметров орбиты на протяжении 50 млн. лет ( рис. 5 а), а также в течение просчитанного нами периода в 100 млн. лет ( рис. 5б). Подобная эволюция параметров орбиты Марса в неподвижных координатах описана и в работе К. Остервинтера и др. [22] (см. рис. 6), а также в работе [4] (см. рис. 9).

Из этого можно заключить, что большие изменения угла наклона I_Mr плоскости орбиты, приведенные в работе [4], вызваны, прежде всего, движением плоскости экватора.

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

Т.Куинн и его коллеги, в работе [3] справедливо отмечают, что численное интегрирование уравнений движения планет будет способствовать развитию новых аналитических теорий, специально разрабатываемых для исследования динамики Солнечной Системы на больших интервалах времени.

Выполненный анализ и сравнение с результатами других исследователей, дают основание считать, что, полученные нами решения, имеют достаточно высокую точность. Мы сравнили наши результаты, относящиеся к планете Земля, с результатами М. Миланковича [23], Ш.Г.Шараф и Н.А.Будниковой [24], А.Берже и М. Лутре [25]. Сравнение показывает, что более поздние по времени результаты лучше соогласуются с нашими результатами.

Считаем целесообразным также отметить и некоторые соображения, касающиеся результатов других исследователей. Из нашего анализа достоверно следует, что теория Д.Брауэера и А. Вуркома А. [21] достаточно хорошо описывает короткопериодические изменения эксцентриситета и долгопериодические колебания наклона плоскости орбиты, однако она менее точно описывает долгопериодические колебания эксцентриситета, а короткопериодические движения плоскости орбиты представлено такими колебаниями, которые, по-видимому, в действительности, отсутствуют. C другой стороны, теория Ж. Ляскара и его коллег [4] хорошо представляет эволюцию эксцентриситета, перигелия и угла наклона орбиты в неподвижном пространстве на интервале времени в несколько десятков миллионов лет.

Выводы

1. Мы вычислили вековые изменения параметров орбиты Марса на промежутке времени в несколько тысяч лет (-3.4 тыс. лет T 3.6) и сравнили их с аппроксимацией, построенной с помощью наблюдательных данных. Результаты совпадают с достаточно высокой точностью. Сравнение численного решения уравнений (3) на интервале времени в 50 лет (с 30.12.1949 г. по 30.12.1999 г.) с положениями планет, взятыми из Эфемерид DE406/LE406 [9], показывает, что в пределах точности, с которой определены массы планет, положения планет совпадают.

2. Были численно проинтегрированы уравнения движения девяти планет, Луны и Солнца на промежутке времени в 100 млн. лет «в прошлое» и исследована эволюция орбиты Марса. На основании этого, можно сделать следующие выводы:

2.1. Эксцентриситет его орбиты испытывает короткопериодические (с периодом Т_е1 = 95.2 тыс. лет) и долгопериодические (с периодом Т_е2 = 2.31 млн. лет) колебания.

2.2. Наклон плоскости орбиты Марса к неподвижной плоскости экватора и её восходящий узел совершают колебания с периодом в 73.1 тыс. лет.

2.3. Перигелий движется в плоскости орбиты в направлении орбитального движения, совершая один оборот в среднем за 78 тыс. лет.

2.4. На всём интервале времени в 100 млн. лет орбита Марса претерпевает малые изменения, т.е. обладает определенной устойчивостью.

3. Сравнение наших результатов с приближёнными теориями других авторов позволяет качественно оценить области применения этих теорий.

4. Движения в рассмотренной модели Солнечной Системы на интервале времени в 100 млн. лет имеют «стабильный» характер и не проявляется какая-либо тенденция к катастрофическим изменениям в геометрии и динамике Солнечной Системы.

Основные вычисления были выполнены на суперкомпьютерах МВС-1000 Сибирского Суперкомпьютерного Центра СО РАН (г. Новосибирск), Института Прикладной Математики им. М.В.Келдыша и Вычислительного Центра им. А.А.Дородницына РАН (г. Москва).

Авторы считают необходимым поблагодарить доктора наук М. Креславского за помощь в получении, прежде всего, материалов, изложенных в работах С. Когена, Е. Хьюббарда и С. Остервинтера, Ж. Ляскара с соавторами и других исследователей. Выполненное количественное сравнение наших расчетов с расчетами Ж. Ляскара с соавторами оказалось возможным благодаря тому, что свои результаты названные ученые разместили на сайте http://www.imcce.fr/Equipes/ASD/insola/mars/La2003-04/. Свободный доступ к результатам Лаборатории реактивного движения NACA на сайте http://ssd.jpl.nasa.gov/ (JPL Solar System Dynamics) позволил нам уточнить начальные условия и сравнить свои расчеты с современными эфемеридами.

В заключение отметим, что данная работа выполнена при поддержке грантов губернатора Тюменской области за 2003-й и 2004-й годы, интеграционной программы Президиума РАН № 13 за 2004г. и гранта РФФИ 04-01-00227.

Приложения

Приложение 1. Цитаты из работ Ж. Ляскара с соавторами

1. «Large-scale chaos in the solar system»

Abstract. Numerous integrations of the solar system have been conducted, with very close initial conditions, totaling an integration time exceeding 100 Gyr. The motion of the large planets is always very regular. The chaotic zone explored by Venus and the Earth is moderate in size. The chaotic zone accessible lo Mars is large and can lead to eccentricities greater than 0.2. The chaotic diffusion of Mercury is so large that its eccentricity can potentially reach values very close to 1, and ejection of this planet out of the solar system resulting from close encounter with Venus is possible in less than 3.5 Gyr.

2. “Long term evolution and chaotic diffusion of the insolation quantities of Mars.”

As the orbital motion is chaotic, even with a precise dynamical model, the computer roundoff numerical error alone will prevent obtaining a precise orbital solution for Mars over more than 60 Myr (see Fig.4.a). Moreover, the obliquity of Mars itself is chaotic, even more chaotic than its orbital motion (Laskar and Robutel 1993. Touma and Wisdom 1993). This will prevent even more drastically obtaining a precise solution for the obliquity over more than 10 to 20 Myr. with the present knowledge of the initial parameters (section 3.2.1).

Our goal in this section will thus be to obtain a solution for the insolation parameters of Mars as precise as possible over 10 to 20 Myr. for use in Mars paleoclimate studies. Then, with the same model, to explore the behavior of the solutions over 250 Myr and to derive a statistical vision of this chaotic system. Although no precise prediction is possible over this time interval, we will be able to derive a precise estimate of the density probability function for the evolution of the eccentricity of Mars and its obliq-uity. It is in fact paradoxical (see for example Lasota and Mackey 1994) that it is actually the chaotic behavior of the system that will allow us to make a precise prediction of the evolution of the density function of Mars' orbital and rotational parameters.

3. “A long-term numerical solution for the insolation quantities of the Earth…”

2. The orbital solution La90. The orbital solution La90 is obtained by the numerical integration of en extended averaged system, which represents the mean evolution of the orbits of the planets. All the 8 main planets of the solar system are taken into account, as well as the mean lunar and relativistic perturbations. The use of numerical integration for competing the solution of the secular system is one of the reasons for the good quality of this solution, which was checked by comparing with the available ephemeris over a short time scale (Laskar 1986, 1988). In (Laskar 1988), the solution La88 was represented in quasi-periodic form over 10 Myr, but these representations are slowly convergent, which prevents good accuracy of the solution.

Later on, the reason for this slow convergence was understood to be due to the presence of multiple resonances in the secular system of the inner solar system (Laskar 1990). Because of these resonances, the motion of the solar system is chaotic, and not quasi-periodic, as was first demonstrated by the computation of its Lуapunov exponents, which reaches l/5 Myr (Laskar 1989). This implies that it is not possible to give any precise solution for the motion of the Earth over more than about 100 Myr, and most probably, ephemerides can only be given with good precision for about 10 Myr to 20 Myr.

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

Табл. 1. Массы m_i и начальные условия (первый вариант): координаты x_i, y_i, z_i и скорости v_xi, v_yi, v_zi планет от Меркурия (Me) до Плутона (Pl) и Луны (Mo) на 30.0 дек. ET 1949 г., JD₀ = 2433280.5 в гелиоцентрической экваториальной системе координат на эпоху 1950.0, JDs = 2433282.4234 при G = 6.67259E-11 м³/(с²·кг) и массе Солнца: 1.98911765646E+30 кг. Положения и скорости - модифицированные данные DE19.

Тела	mi, кг	Координаты, м	Скорости, м/с
		xi	yi	zi	vxi	vyi	vzi
Me	3.32462E+23	51451371876.2197	6823568174.75167	-1634415569.98258	-14659.2853816947	44351.5132622098	25256.7096329313
Ve	4.869E+24	21388391208.9948	96791808863.7023	42259236150.9888	-34445.734417752	5390.50795859843	4604.86892148603
Ea	5.9742E+24	-20403349572.3061	133649779200.966	57962499498.9961	-29982.2390350433	-3894.92114820137	-1690.42257915082
Ma	6.4191E+23	-204926756004.512	126133319400.949	63405879798.2523	-12786.2796102665	-16409.8466816399	-7188.5117548751
Jp	1.8988E+27	501063754034.797	-519675392681.743	-235195315982.077	9671.45813229624	8592.04165989274	3449.57181755033
Sa	5.685E+26	-1342285170665.77	341044959330.524	199023100858.158	-3217.56963808403	-8629.44636855088	-3428.76469351675
Ur	8.6625E+25	-150035825137.752	2591594394295.33	1137683474058.44	-6848.45451198089	-650.865182644724	-188.460863014761
Ne	1.0278E+26	-4367462454310.03	-1154797346757.77	-363116492414.853	1421.18239257454	-4799.82638972212	-2001.79537726796
Pl	1.09775E+24	-3924270721583.65	3076004854807.31	2160776068098.85	-2278.14800627023	-4536.19858950397	-739.431338563464
Mo	7.35E+22	-20099352853.0859	133857559955.206	58067973377.9	-30562.0814206052	-3153.29849566654	-1270.85380799933

Табл. 2. Массы m_i и начальные условия (второй вариант): координаты x_i, y_i, z_i и скорости v_xi, v_yi, v_zi планет от Меркурия (Me) до Плутона (Pl) и Луны (Mo) на 30.0 дек. ET 1949 г., JD₀ = 2433280.5 в гелиоцентрической экваториальной системе координат на эпоху 2000.0, JD_S = 2451544 при G = 6.67259E-11 м³/(с²·кг) и массе Солнца: 1.98891948976803E+30 кг. Положения и скорости - данные DE406/LE406.

Тела	mi, кг	Координаты, м	Скорости, м/с
		xi	yi	zi	vxi	vyi	vzi
Me	3.30187842779737E+23	51378500712.7917	7398730258.73134	-1384413714.5863	-15276.8167285921	44183.2585096919	25183.7713637451
Ve	4.86855338156022E+24	20098437756.1693	97022456883.2422	42359548922.943	-34525.343689039	5004.73164455778	4437.16576036982
Ea	5.97369899544255E+24	-22177799063.1322	133409762224.645	57858151723.3985	-29927.5728090066	-4228.84376557956	-1835.38708532847
Ma	6.4185444055007E+23	-206626992480.617	123830934957.816	62404742310.1315	-12566.7488056756	-16551.3429850097	-7249.99987165166
Jp	1.89900429500553E+27	507971625555.693	-514027686416.492	-232739248300.242	9557.787142254	8699.41711715317	3496.24926688155
Sa	5.68604198798257E+26	-1346946949935.64	326006202166.041	192483689461.594	-3104.13240669983	-8664.66647884277	-3444.08286511961
Ur	8.68410787490547E+25	-184532065685.545	2589691250995.06	1136852552974.41	-6839.65008976224	-727.382755054996	-221.718557104065
Ne	1.02456980223201E+26	-4352452287291.9	-1203523662980.49	-384291106982.115	1484.35872636776	-4783.5269899047	-1994.70645836478
Pl	1.65085753263927E+22	-3969057713524.92	3031889154729.01	2141676186885.16	-2223.92294605229	-4561.78400082869	-750.399172063762
Mo	7.34767263035645E+22	-21860138201.1275	133632383977.255	57970881141.0125	-30536.9357540009	-3571.10974187999	-1466.38653570499

Табл. 3. Сопоставление результатов интегрирования уравнений (3) за 50 лет, с 30.0 дек. 1949 г. JD0 = 2433280.5 до 30.0 дек. 1999 г. JDf = 2451542.5 (верхние числа) с данными DE406/LE406 на ту же конечную дату (нижние числа). Координаты x_i, y_i, z_i и скорости v_xi, v_yi, v_zi планет от Меркурия (Me) до Плутона (Pl) и Луны (Mo) в гелиоцентрической экваториальной системе координат на эпоху 2000.0, JD = 2451544; Дц_i - проекция на плоскость экватора разности углового положения тел на орбите.

Тела	Координаты, м	Скорости, м/с	Дцi, радианы
	xi	yi	zi	vxi	vyi	vzi
Me	-27244931161.36578 -27250751819.76618	-57547047079.05123 -57544485667.38114	-27913749550.38322 -27911777753.74239	35013.83456495187 35011.77049279869	-13510.56095560906 -13515.3275344842	-10848.24752917202 -10850.57935386931	9.851E-5
Ve	-107488775402.6213 -107488764214.6746	4508043.433882821 4745110.048553557	6804610208.31456 6804716122.693854	-1083.587133251849 -1083.670285579284	-32097.14673367943 -32097.14769414636	-14371.49148962517 -14371.48662231141	2.171E-6
Ea	-20039777471.25795 -20040383501.21171	133711555559.6049 133711493201.8992	57970416797.15794 57970388317.83772	-29998.02022870762 -29998.03736247704	-3814.898490683077 -3815.022803925509	-1654.426543729306 -1654.476364250177	-4.501E-6
Ma	207725381980.5454 207725522951.1265	-4958199493.666488 -4956233639.810995	-7890749392.650456 -7889850986.827501	1824.915025362562 1824.655481281149	23910.81526617862 23910.82245919431	10917.78330229408 10917.793483078	-9.518E-6
Jp	600272372864.9594 600271538043.9522	407182549280.8615 407183656391.9888	159908186034.7253 159908680579.1977	-7867.863761424506 -7867.884321515599	10212.08818437833 10212.07391886902	4568.95256774887 4568.946941038761	-1.910E-6
Sa	959990440191.2358 959989658746.129	922394198483.991 922394901617.7257	339686323126.5457 339686647213.0083	-7421.340204237804 -7421.345521813606	6119.467955244539 6119.462360989111	2846.724756339436 2846.722675206618	-7.870E-7
Ur	2157973036652.258 2157973612539.625	-1871830749562.907 -1871830390715.251	-850356722393.0418 -850356575438.5552	4639.399347233582 4639.397946737375	4260.85932689923 4260.860147703253	1800.516264157553 1800.516647445246	-2.269E-7
Ne	2514082274539.627 2514083052562.011	-3438399709289.256 -3438399345461.22	-1469944545800.924 -1469944416163.849	4466.730071503929 4466.72940525403	2887.755850345153 2887.756455291418	1071.003275988534 1071.003539817953	-1.979E-7
Pl	-1478463862173.61 -1478463046087.412	-4185154424340.236 -4185154855001.015	-860265513672.4745 -860265895224.6543	5243.729378621804 5243.729531171102	-1953.631698552901 -1953.630796975506	-2190.768497175936 -2190.768268488547	-2.060E-7
Mo	-20421412252.06083 -20421819304.72665	133620881993.3932 133619777434.6888	57967399445.62585 57967006904.80833	-29839.60237015856 -29836.74314646615	-4732.942291554387 -4732.402371602139	-2010.441015676699 -2010.541969780811	-4.276E-6

Приложение 3. Исследование достоверности результатов интегрирования уравнений (3).

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

1. Контроль сохранения общего момента количества движения M Солнечной Системы. В силу существования десяти первых интегралов системы дифференциальных уравнений (3) [ 8 ], например, проекция момента М_z на ось z имеет вид

, (41)

где n - количество тел( в нашей модели n= 11).

Если наблюдается изменение М, это свидетельствует о погрешности вычислений. При шаге интегрирования t=1·10^-4 года, относительное изменение проекции момента на ось z после 10 тыс. лет составляет дM_z = 1.4·10^-14, а после 100.08 млн. лет - дM_z = 7.9·10^-11. Мы установили, какими могут быть погрешности в движении тел при вышеупомянутом изменении момента. Например, при шаге счета t = 1·10^-4 года и дM = -4.5·10^-14 смещение тела по орбите за 10 тыс. лет изменяется в радианах (см. табл. 3) от значения -4.94·10^-10 для Плутона до значения 3.44·10^-4 для Луны. Если изменение параметров орбиты за этот период значительно превышают эти погрешности, то такие результаты нельзя считать достоверными.

2. Контроль количества движения Солнечной Системы Р. Аналогично его проекция , например, на ось х имеет вид:

. (42)

Отличие величины Р от нуля также свидетельствует о накоплении погрешностей вычислений. В начальный момент, из-за округления при переходе к двоичным кодам, величина Р отличается от нуля и равна 5·10^-18. После реализации счета на 10 тыс. лет Р= 2.4·10^-15, а после 100.08 млн. лет - Р= 3.6·10^-13. Здесь и в дальнейшем приведены безразмерные величины.

3. Учет влияния разных точностей счета. Если при повышении точности счета результаты повторяются, то их можно считать достоверными.

4. Вычисление параметров движения тел в « прошлое» и в « будущее». Сравнивая два решения (в «прошлое» и в «будущее») при «безошибочных» вычислениях не должны появиться какие-либо скачки в изменениях параметров в начальный момент времени.

5. Расчеты в отдаленную по времени эпоху и возвращение в исходную. Это позволяет выявить все погрешности вычислений, их величину и структуру. Было установлено, что основная погрешность заключается в смещении тела на орбите (например, от 1.5° для Меркурия до 4.5°·10^-9 для Плутона при счете на 10 тыс. лет с шагом t=2·10^-4 года) без существенного изменения орбиты. Наибольшее смещение в этом случае наблюдается у Луны, и оно достигает 180° относительно Земли. При шаге счета t = 2·10^-5 года смещение Луны по ее орбите равняется 2'. Следует отметить, что этот тест дает завышенные погрешности из-за разрыва начальных производных при обратном счете.

6. Использование тестовых задач, имеющих точное аналитическое решение. Такой задачей является, например, осесимметричное взаимодействие n-тел, находящихся в одной плоскости [15]. Этот метод позволяет выявить вычислительные погрешности, прежде всего, в самой тестовой задаче и было установлено, что погрешность увеличивается с увеличением эксцентриситета орбиты. При эксцентриситете, равном эксцентриситету Меркурия, в задаче 12-ти тел смещение их по орбите за 10 тыс. лет равняется 0.16" при шаге t=2·10^-4 года.

7. Сравнение с данными наблюдений.

7.1. Сравнение координат и скоростей тел на интервале в несколько десятков лет. Например, мы вычислили положение тел на дату 30.12.1999,0 г по начальным данным, известным на 30.12.1949,0 г. и получили, что средняя «по всем телам» относительная погрешность в положениях небесных тел за 50 лет составила ?r_?? = 1.14·10^-5.

Такое сопоставление позволяет выявить погрешность в начальных условиях, в исходных данных (например, в массах тел), а при высокой точности счета и погрешность в наблюдениях. 7.2. Сравнение расчетной траектории Солнца с траекторией, определенной по эфемеридам планет за 60 лет. Мы получили практически совпадение расчетов [6] с данными эфемерид.

Б. Смещение по орбите. Существует достаточно богатое разнообразие способов проверки результатов интегрирования уравнений. Наши исследования и исследования других авторов, показали, что наибольшая погрешность в решениях проявляется в смещении тела по его орбите. В табл. 3 представлены результаты наших исследований о зависимости смещений тел по орбите Дц_i (точнее, проекция смещения на плоскость экватора) от относительной погрешности в моменте количества движения Солнечной Системы дM. Величины ?_i определяются по формуле

Дц_i = arctg(y_i/x_i) - arctg(y_0i/x_0i), (43)

где y_0i, x_0i, y_i, x_i - координаты тел в барицентрической экваториальной системе координат;

y_0i, x_0i - координаты тел в решении с наименьшим значением дM = 1.46?10^-14 при Дt = 1?10^-4 года и с удвоенной точностью;

Дц_i - проекция смещения i-того тела по орбите на плоскость экватора.

В последней строке таблицы 3 дано смещение Луны относительно Земли (MoE), которое определяется по формуле:

Дц_MoE = arctg[(y₁₀- y₃)/(x₁₀- x₃)] - arctg[(y₀₁₀- y₀₃)/(x₀₁₀- x₀₃)], (44)

где индекс 10 относится к Луне, а индекс 3 - к Земле.

В расчетах, представленных в столбцах 1 - 5, последовательно увеличивался шаг dt, вследствие чего погрешность момента ?M увеличивалась. При увеличении шага на порядок погрешность момента возросла на 5 порядков. Погрешность положения ??_i, как видно из таблицы 4, для разных тел изменяется по-разному. Для внешних планет она практически не изменяется, а для внутренних планет и, особенно для Луны, погрешность положения растет с увеличением относительной погрешности момента дM. Например, для Меркурия ??_i возрастает на 8 порядков, поэтому величину дM можно считать интегральной характеристикой точности вычислений. На этом основании мы составили таблицу смещений по орбитам в зависимости от величины изменения момента для одиннадцати тел плюс смещение Луны относительно Земли.

Табл. 4. Смещение тел по их орбитам Dj_i за -10 тыс. лет в решениях с разным изменением момента ?M по отношению к решению с расширенной длиной числа (34 десятичных знака) при dt =1?10^-4 года и ?M= 1.46?10^-14: расчет в прошлое от эпохи 1949, дек. 30.0 ET, JD₀ = 2433280.5;--Dj_I - в радианах; планеты от Меркурия (Me) до Плутона (Pl), Mo - Луна, Su - Солнце; MoE - смещение Луны относительно Земли.

Тело	i	1	2	3	4	5
		?M =-4.5?10-14 dt = -1?10-4	?M = 2.3?10-13 dt = -4/3?10-4	?M = 3.4?10-13 dt = = -2?10-4	?M = 4.7?10-11 dt = -5?10-3	?M = 1.5?10-9 dt = -1?10-3
Me	1	-4.78076E-8	0.000813003	0.00784215	0.880401	6.234269
Ve	2	-4.41872E-8	2.13205E-7	2.93696E-6	0.000296793	0.00947411
Ea	3	-1.21898E-7	-1.00427E-5	-6.11468E-5	-3.43795E-5	-0.000247816
Ma	4	-3.21748E-8	8.11215E-9	5.47077E-10	1.2199E-7	4.19873E-6
Jp	5	-2.85178E-9	-8.8565E-10	-1.57537E-9	-1.95009E-10	3.13616E-12
Sa	6	5.65911E-10	1.46593E-9	3.15613E-9	-1.40978E-9	1.44393E-10
Ur	7	2.02969E-11	-5.09915E-10	7.50903E-10	-6.51694E-11	-7.28375E-11
Ne	8	-7.06222E-10	3.01113E-10	-3.27041E-10	-6.95799E-12	-2.10711E-10
Pl	9	-4.93649E-10	8.25153E-11	-8.58744E-11	-2.79003E-11	-1.09287E-10
Mo	10	-2.52436E-7	0.000822199	0.00496378	0.00238666	0.00119156
Su	11	-1.79509E-9	3.66315E-9	4.20909E-08	7.5401E-6	2.03371E-6
MoE	12	0.000344309	0.941297	3.07668	7.84055	11.41+2?k??

В. Сравнение результатов интегрирования уравнений (3) с эфемеридами. Уравнения (3) были проинтегрированы с шагом с dt = 1·10^-4 года со вторым вариантом начальных условий по DE406/LE406 (см. табл. 2) в эпоху JD₀ = 2433280.5 за 50 лет до 00 h ET 30.12.1999 г. с JD_f = 2451542.5 и сопоставлены с данными DE406/LE406 на ту же дату. В табл. 4 приведены результаты сравнения декартовых координат и скоростей тел, а также угловых смещений ц_i, рассчитанных нами положений тел относительно положений тел, данных в «эфемериде» DE406/LE406. Величина ц_i была вычислена по формуле (43) при y_{0i =} y_ji и x_0i = x_ji - координаты тел- по данным из « DE406/LE406» в барицентрической экваториальной системе координат. Для Земли это смещение ц₃ = -0.93ґґ. Среднее относительное отклонение положений тел ?r_?? = 1.14·10^-5 (по 11 небесным телам), где

dr_??=

- барицентрический радиус-вектор рассчитанного положения i-того тела;

- такой же вектор i-того тела, полученный из эфемерид DE406/LE406.

Время счета в юлианских столетиях определялось по формуле

T =(JD_f - JD₀)/(k_T ·36525),

где k_T = 1.000017537231278 - коэффициент приведения расчетного времени к эфемеридному времени. Отличие k_T от 1 объясняется различиями принятых при интегрировании исходных данных от действительных.

Аналогично мы проинтегрировали уравнения на промежутке времени в 50 лет со вторым вариантом начальных условий, с шагом интегрирования dt = 1·10^-5 года и с расширенной длиной числа. Погрешность момента в этом случае получилась равной ?M_z = 7.39·10^-22 и по сравнению с погрешностью в предыдущем решении, ?M_z = 1.23·10^-14. Таким образом, точность улучшилась на восемь порядков, и при этом результаты отклонений не изменились.

Отметим также, что сравнение результатов интегрирования уравнений (3), выполненное с первым вариантом данных по каталогу DE19, с данными на дату 8.01.1994 г. по Астрономическому ежегоднику, дало среднее отклонение положения ?r_? = 1.62·10^-3 и k_T = 1.0001337. Таким образом, второй вариант исходных данных и начальных условий привел к уменьшению отклонения положения на два порядка и уменьшению отличия коэффициента k_T от 1 на порядок. Отсюда можно заключить, что второй вариант начальных условий и данные, взятые из эфемерид DE406/LE406, являются более точными.

Как было отмечено ранее, отличие между массами планет значительно больше достигнутой величины среднего значения (по всем телам) относительного отклонения положений тел ?r_?? = 1.14·10^-5. Например, неопределенность гравитационной постоянной G согласно Конвенции IERS2003 (см. сайт IERS Conventions Center. IERS Conventions (2003) - http://tai.bipm.org/iers/conv2003/conv2003.html) составляет ?G = 1.5·10^-4. Поэтому можно считать, что в пределах точности задания масс результаты расчетов совпадают с эфемеридами DE406/LE406.

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

Литература

1. Laskar J. Large-scale chaos in the solar system // Astron. Astrophys. - 1994. - Vol. 287, L9 - L12.

2. Laskar J. Marginal stability and chaos in the Solar system / Ferraz Mello S. et al. (eds.), Dynamics, ephemerides and astrometry of the Solar System. - IAU: Netherlands. - 1996. - P. 75 - 88.

3. Quinn T. R., Tremaine S., Duncan M. A three million year integration of the earth's orbit // Astronomical Journal. - 1991. - V. 101. - P. 2287-2305.

4. Laskar J., Correia A. C. M., Gastineau M., Joutel F., Levrard B. and Robutel P. Long term evolution and chaotic diffusion of the insolation quantities of Mars // Icarus.- 2004.- Vol. 170, Iss. 2.- P. 343-364.

5. Смульский И.И. Теория взаимодействия. - Новосибирск: Из-во Новосиб. ун-та, НИЦ ОИГГМ СО РАН, 1999. - 294 с.

6. Мельников В.П., Смульский И.И., Кротов О.И., Смульский Л.И. Орбиты Земли и Солнца и возможные воздействия на криосферу Земли (постановка проблемы и первые результаты// Криосфера Земли. - 2000 - Т. IV, №3, с. 3-13.

7. Мельников В. П., Смульский И.И. Астрономические факторы воздействия на криосферу Земли и проблемы их исследования// Криосфера Земли. - 2004. - Т. VIII, № 1, с. 3-14.

8. Справочное руководство по небесной механике и астродинамике / Под ред. Г. Н. Дубошина. Изд. 2-е, доп. и перераб. М., Наука, 1976, 862 с.

9. Стэндиш Е.М. (Standish E.M.) JPL Planetary and Lunar Ephemerides, DE405/LE405.// Interoffice memorandum: JPL IOM 312. F - 98-048. August 26. 1998. (ftp://ssd.jpl.nasa.gov/pub/eph/export/DE405/).

10. Труды ИПА РАН. Вып. 10. Эфемеридная астрономия. - Санкт-Петербург: ИПА РАН, 2004 - 488 с.

11. Бордовицына Т.В. Современные численные методы в задачах небесной механики. М.: Наука, 1984, 136 с.

12. Стеффенсон К.Ф. (Steffensen K.F.) On the problem of Three Bodies in the Plane. - Kong. Danske Videnskab. Seskab., Mat. Fys. Medd. 1957, v. 13, N 3, 18 p.

13. Яров-Яровой М.С. О применении уточненных методов численного интегрирования в небесной механике. Труды ГАИШ. - М., 1974, т. 45, с. 179-200.

14. Гребеников Е.А., Козак-Сковородкина Д. и Якубяк М. Методы компьютерной алгебры в проблеме многих тел. - Москва: Из-во Российского университета дружбы народов. - 2001 г. - 213 с.

15. Смульский И.И. Осесимметричная задача гравитационного взаимодействия N-тел// Математическое моделирование. - 2003, а, т. 15, № 5, с. 27-36.

16. Гребеников Е.А. Существование точных симметричных решений в плоской ньютоновской проблеме многих тел//Математическое моделирование.-1998.-Т.10, N8. - С.74-80.

17. Смульский И.И. Новая геометрия эволюции орбит // Proceeding of Joint International Scientific Conference “New Geometry of Nature, August 25 - September 5, 2003, б, Kazan State University, 2003. - с. 192-195.

18. Смульский И.И. Расчет взаимодействий в Солнечной системе за 50 млн. лет для изучения эволюции климата // «Большая Медведица». Журнал проблем защиты Земли. - Межрегиональный общественный фонд им. Ломоносова. Новосибирск, Россия.- 2005 г. No. 1, с. 44-56.

19. Newcomb S. The elements of the fourt inner planrts a nd the fundamental constants of astronomy. Washington: Government printing office. 1895. -202 p.

20. Симон Дж.Л. и др. (Simon J.L., Bretagnon P., Chapront J. et. al.) Numerical Expression for Precession Formulae and Mean Elements for the Moon and the Planets // Astron. Astrophys. - 1994, vol. 282, p. 663-683.

21. Браувер Д. и Вурком А. (Brouwer D., Van Woerkom A. J. J.) The secular variation of the orbital elements of the principal planets// Astr. Pap. - 1950. - 13, 2.

22. Коген, Хьюббард и Остервинтер (Cohen C.J., Hubbard,E. C.; Oesterwinter C.) Planetary Elements for 10,000,000 Years//Celestial Mechanics. - 1973. - No. 3. -Pp. 438-448.

23. Миланкович М. Математическая климатология и астрономическая теория колебаний климата. - М.-Л. -ГОНТИ. - 1939. -207 с.

24. Шараф Ш. Г. и Будникова Н. А. О вековых изменениях элементов орбиты Земли, влияющих на климаты геологического прошлого. // Бюл. ИТА АН СССР. 1967, вып. 11, № 4 - С. 231 - 261.

25. Берже А. и Лоутре М.Ф. (Berger A. and Loutre M. F.) Insolation values for the climate of the last 10 million years// Quaternary Science Reviews. - 1991. - № 10, Р. 297 - 317.