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

Рассматривается астрономическая теория ледниковых периодов на базе решения неупрощенных дифференциальных уравнений движения современными численными методами. Проанализированы результаты астрономических теорий ледниковых периодов и обоснована целесообразность решения этой проблемы более точными методами. Изменение инсоляции Земли зависит от эволюции ее двух движений: орбитального и вращательного. Разработан численный метод повышенной точности для решения задачи о первом движении и приведены результаты эволюции орбиты Земли в течение несколько миллионов лет. Получены периоды изменения параметров орбиты Земли и независимым методом подтверждены результаты предшествующих исследователей.

Ледниковые периоды, Солнечная система, уравнения движения, точность решения, периоды колебаний орбиты.

The astronomical theory of glacial ages is considered on the basis of the decision of the not simplified differential equations of motion by modern numerical methods. The results of astronomical theories of glacial ages are analyzed and the expediency of the investigation of this problem by more exact methods is proved. The change of insolation of the Earth depends on evolution of its two movements: orbital and rotary. The numerical method of the higher accuracy for the solution of the problem on the first movement is developed and results of evolution of the Earth's orbit for some millions years are given. The periods of parameters change of the Earth's orbit are received and the independent method confirms results of previous researchers.

Для объяснения колебаний климата Земли в геологические эпохи привлекались разнообразные астрономические факторы, в том числе прямое падение внешних тел на Землю или Солнце. Однако в настоящее время существуют только два фактора, которые поддаются строгому научному анализу. Это изменения собственного вращательного движения Земли и движения ее по орбите под воздействием тел Солнечной системы. Такая астрономическая теория ледниковых периодов получила достаточно широкое обоснование в работах М. Миланковича в 1920-1930 гг. [Миланкович, 1939]. При расчете эволюции параметров орбиты Земли М. Миланкович и последующие исследователи: Брауэр и Вурком [Brouwer and Woerkom, 1950], Шараф и Будникова [Шараф и Будникова, 1967], Берже и Лоутре [Berger and Loutre, 1991] и др. основывались на теории вековых возмущений У. Леверье (1843г.). Это метод приближенного аналитического решения уравнений взаимодействия тел Солнечной системы, при котором среднее движение планеты аппроксимируется движением по эллиптической орбите, а для разности с действительной орбитой получают упрощенные дифференциальные уравнения. В настоящее время решение уравнений в виде эволюции эксцентриситета орбиты, наклона ее плоскости и смещения перигелия получены для периода до 30 миллионов лет [Шараф и Будникова, 1969].

Этими исследованиями установлен характер воздействия параметров орбиты Земли на ее климат. При эксцентриситете орбиты, отличном от нуля, Земля в перигелии наиболее близко расположена к Солнцу, поэтому она будет получать тепла больше, чем в афелии. Например, при современном эксцентриситете e=0,0167 инсоляция Земли в перигелии на 7% больше, чем в афелии. Кроме того, за счет меньшей скорости в афелии период года, приходящийся на более холодный афелий, удлиняется. В результате с увеличением эксцентриситета уменьшается общее количество тепла, получаемое всей Землей.

При увеличении угла наклона орбиты Земли к плоскости ее экватора количество тепла, приходящее на полярные широты, увеличивается. Например, увеличение в современную эпоху угла наклона на 1^о приводит к увеличению инсоляции на широте 65^о северного полушария на 1,4%.

Положение перигелия относительно восходящего узла орбиты Земли оказывает влияние на среднюю температуру сезона, приходящегося на перигелий. В современную эпоху Земля находится в перигелии во время зимы в северном полушарии. Поэтому на широте 65^о количество тепла в северном полушарии поступает на 2,2% больше, чем в южном на той же широте. В связи с этим зима северного полушария теплее зимы южного.

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

На рис. 1 показано изменение инсоляции на широте 65^о северного полушария за миллион лет по расчетам разных авторов. За этот период количество тепла изменялось от величин, получаемых Землей на широте 75^о, до количества тепла, идентичных широте 60^о. Изменение инсоляции происходит с нерегулярными периодами и амплитудой, при которых возможны выбросы до очень сильных похолоданий и потеплений.

Из графиков видно, что максимум последнего потепления в северном полушарии произошел 10 тыс. лет назад. И с тех пор идет непрерывное похолодание, которое, как следует из результатов исследований [Шараф и Будникова, 1969] будет продолжаться еще 10 тыс. лет. В современную эпоху Земля находится посредине между крайними изменениями инсоляции.

Эти результаты согласуются с известными данными по истории Земли. Потепление от современного уровня (см. рис. 1), начавшееся 20 тыс. лет назад, привело к таянию ледника в Европе, и к моменту наступления максимума инсоляции ледник полностью исчез. В северной части Азии в период максимума инсоляции наблюдались благоприятные условия для существования крупных млекопитающих: мамонтов, носорогов, бизонов, оленей, лошадей и др. [Сулержицкий и Романенко, 1997]. Находки останков наиболее поздних останков животных датируются 3-4 тысячами лет, т.е. еще до этого времени сохранялись необходимые для их существования жизненные условия. Возможно, что дальнейшее похолодание и вызванные им природные изменения привели к полному исчезновению этих животных.

Рис. 1. Изменение инсоляции летом для широты 65^о северного полушария по данным разных исследователей: а - М. Миланковича (1939 г.); б - А. Вуркома (1950 г.); в - Ш. Г. Шараф и Н. А Будниковой (1969 г.); г - А. Берже и М. Ф. Лоутре (1991 г.). По оси абсцисс отложено время в тысячелетиях от 1950 г.; по оси ординат: а, б, в - инсоляция в эквивалентных широтах в течение летнего полугодия [Шараф и Будникова, 1969], г - среднемесячная инсоляция в июле в Вт/м² [Berger and Loutre, 1991].

Из рис. 1 видно, что расчеты инсоляции разных авторов в пределах первой сотни тысяч лет согласуются, затем между ними начинают нарастать вначале количественные изменения, а затем и качественные. В связи с этим можно сделать следующие выводы. Во-первых, рассмотренные астрономические факторы могут приводить к изменению климата. Во-вторых, так как результаты разных авторов отличаются, то необходимо установить действительные климатические изменения, обусловленные этими факторами. Тогда их можно будет сопоставлять с историей Земли, интерпретировать ее и датировать с высокой точностью.

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

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

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

При аналитическом приближенном решении этих двух задач были введены упрощения в физическом и математическом планах. Поэтому авторы проведенных исследований остерегаются продолжать расчеты инсоляции Земли на большие периоды времени. Так М. Миланкович [Миланкович, 1939] полагает, что надежные результаты могут быть получены на период 600 тыс. лет. Шараф и Будникова [Шараф и Будникова, 1969], уточнив исходные данные и прецессионные решения, рассчитали инсоляцию на 30 млн. лет, Берже и Лоутре [Berger and Loutre, 1991] - на 5 млн. лет. Ж. Ляскар и др. (Laskar at al., 1993) усовершенствовали теорию вековых возмущений, получили решения до 200 млн. лет, но за их хаотического поведения ограничились результатами в пределах -20 млн. лет до +10 млн. лет.

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

1) какие изменения орбит потеряны в результате упрощения уравнений?

2) не обусловлены ли полученные изменения орбит упрощениями уравнений?

При численном решении задачи можно избежать многих упрощений физического плана и решать задачу практически без математических упрощений, т.е. получить все параметры движения и для всех тел. Так что первая проблема здесь устраняется. Вторая проблема, проблема достоверности результатов, может быть решена повышением точности счета. В связи с развитием вычислительной техники в последние десятилетия рядом авторов предприняты численные решения задачи орбитального движения. В наиболее полном виде этим способом получили решения Т. Кинн и др. [T.Quinn at all, 1991] на период времени -3.056 млн. лет. Приведем дополнительные два аргумента в целесообразности выполнения численных решений, которые приводят эти авторы, заключаются в необходимости: 1) проверки результатов, полученных в рамках теории вековых возмущений выводов о хаотичности движения внешних планет или всей Солнечной системы; 2) дальнейшего развития аналитических теорий по динамике Солнечной системы за большие периоды времени.

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

Взаимодействие тел Солнечной определяется законом тяготения Ньютона, который на протяжении трех столетий проходит непрерывную проверку. В 20^ом веке начали применять дополнительные слабые воздействия на тела, среди которых следует упомянуть следующие: воздействие астероидов, комет и спутников планет; релятивистскую добавку тяготения; приливное трение; давление света; сопротивление межпланетной среды; изменение масс тел; действие инерциальных сил, обусловленных орбитальным движением Солнечной системы и др. Однако, величины этих сил не определены с достаточной точностью, а обоснованность многих из этих воздействий далека от завершения. Поэтому целесообразно получить достоверную эволюцию Солнечной системы, которая обусловлена законом тяготения Ньютона. В последующем, по мере уточнения дополнительных воздействий, их можно будет накладывать на основное воздействие. Таким образом можно достоверно определить значимость каждого из этих дополнительных слабых воздействий.

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

Рассматривается [Смульский, 1999; Мельников, Смульский и др., 2000] гравитационное воздействие 9 планет, Солнца и Луны. Согласно закону всемирного тяготения тело с номером k притягивает тело с номером i силой

Если тел n, то на тело i-e остальные будут воздействовать силой

Под воздействием этой силы, в соответствии со вторым законом механики , i-ое тело будет двигаться относительно инерциальной системы отсчета с ускорением

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

Соотношение (3) представляет систему 3n нелинейных дифференциальных уравнений, где n= 11. Для ее решения мы задаем 3n значений координат и 3n значений компонент скорости на определенную дату, например, 30 декабря 1949 г. [Справочное руководство…, 1976]. Задача решается в барицентрической экваториальной системе координат.

Существует большое разнообразие численных методов. В результате численных экспериментов и их анализа мы пришли к выводу, что конечно-разностные методы интегрирования не обеспечивают необходимую точность. Для решения задачи мы разработали алгоритм и программу Galactica, которая постоянно совершенствуется. Значение функции в следующий момент времени t=t₀ + Dt определяется с помощью ряда Тейлора, который, например, для координаты x имеет вид:

Значение скорости x' определяется по аналогичной формуле, а ускорение x₀'' - по формуле (3). Более высокие производные x₀^(k) определены аналитически в результате дифференцирования уравнений (3). Сейчас используется расчетная схема шестого порядка точности, т.е. при K=6. Используются числа с двойной точностью, при которой они выражается 17 десятичными знаками. Решение за большие сроки выполняется на суперкомпьютерах RM-600 и МВС-1000 в Новосибирском ВЦ СО РАН. Через каждые 10 тыс. лет данные о положении тел записываются в файл. Для анализа параметров орбиты по этим данным просчитывается один оборот планеты вокруг Солнца [Мельников, Смульский и др., 2000], и по результатам расчетов определяются два параметра плоскости орбиты: угол наклона i и положение восходящего узла j_W и восемь параметров орбиты в том числе: эксцентриситет e и угол положения перигелия j_р.

В процессе исследований результаты решений сопоставлялись с данными наблюдений разного уровня и обнаруженные несоответствия устранялись повышением точности счета. Сопоставления проводились для Земли, Солнца и других планет. Все известные особенности движения тел подтверждаются результатами расчетов. Различные проверки показали, что методика численного решения, алгоритм и реализованная программа обеспечивают необходимую точность счета с большим запасом. Основные результаты получены при шаге счета Dt=2*10^-4 года и двойной длиной числа. Проверочные и уточняющие расчеты выполнялись с меньшими шагами, а также с расширенной длиной числа, при которой оно выражается 34-я десятичными знаками.

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

3.1. Контроль момента количества движения системы M, проекция которого, например, на ось z имеет вид:

Так как величина момента для изолированной системы не должна изменяться, то изменение М свидетельствует о погрешности расчетов. Сейчас относительное изменение момента за 10 тыс. лет составляет 0.6*10^-12. В результате расчета при разных точностях счета установлено, какие могут быть погрешности в движении тел при вышеупомянутом изменении момента. Например, при шаге счета Dt=2*10^-4 года относительная погрешность положения изменяется от 0.7*10^-10 для Плутона до 0.8*10^-3 для Луны. Если изменение параметров орбиты за этот период значительно превышают погрешности, то такие результаты могут являться достоверными.

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

Величина Р должна быть равна нулю. Поэтому отличие Р от нуля свидетельствует о накоплении погрешностей округления. В начальный момент из-за округления при переходе к двоичным кодам величина Р отличается от нуля и равна 5·10^-18 . После 10 тыс. лет счета Р= 5·10^-15, после 2.56 млн. лет — Р= 2·10^-14 . Здесь и в дальнейшем приведены безразмерные величины.

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

3.4. Расчет движения тел в прошлое и будущее. Наличие излома в изменении параметра в начальный момент времени свидетельствует о погрешности результатов.

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

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

3.7. Просчет тестовых задач, имеющих точное аналитическое решение. Такой задачей является осесимметричное взаимодействие n-тел, находящихся в одной плоскости [Смульский, 1999]. Этот метод позволяет также выявить все вычислительные погрешности, но только в отношении к тестовой задаче. Было установлено, что погрешность увеличивается с увеличением эксцентриситета орбит. При эксцентриситете, равном эксцентриситету Меркурия, в задаче 12 тел смещение их по орбите за 10 тыс. лет равняется 0.16".

3.8.1.Сопоставление координат и скоростей тел на интервале несколько десятков лет. Например, рассчитывалось положение на 7.01.94 г по начальным данным на 30.12.49г. Средняя по всем телам погрешность положения за 44 года составила 0.16%. В дальнейшем точность расчетов была увеличена на 5 порядков, но погрешность положения не изменилась.

Такое сопоставление позволяет выявить погрешность начальных данных, исходных данных (например, масс тел), а при высокой точности счета и погрешность данных наблюдений.

3.8.2.Сопоставления изменения параметров орбит планет с вековыми возмущениями. В настоящее время выполнено для орбит большинства планет и получено совпадение.

3.8.3.Сопоставление расчетной траектории Солнца с траекторией, определенной по эфемеридам планет за 60 лет. Получено [Мельников, Смульский и др., 2000] совпадение расчетных результатов с данными наблюдений.

3.9. Сопоставление с вычислительными результатами других авторов. Выполнено сопоставление с результатами приближенных аналитических решений и получено согласие с ними.

Уравнения (3) многократно интегрировались за разные периоды времени, с разными алгоритмами счета и разными шагами. Ниже представлены результаты расчетов по последней версии программы Galactica с шагом Dt=2·10^-4 года. Расчеты были выполнены на 1.2 млн. лет в будущее 2.56 млн. лет в прошлое.

На рис. 2 представлено изменение эксцентриситета е, долготы восходящего узла j_W, угла наклона плоскости орбиты і, и положения перигелия j_р за 3.7 млн. лет. Время Т представлено в юлианских столетиях. Эксцентриситет испытывает короткопериодические изменения со средним периодом Те = 94 тыс. лет вокруг его среднего значения е = 0.028. Кроме того, наблюдаются более длительные колебания с периодом 4Те, которые приводят к достижению крайних значений эксцентриситета е = 0.0022 и е = 0.062. Долгота восходящего узла j_W отсчитывается в плоскости экватора 1950 г. от точки весеннего равноденствия на ту же эпоху. Как видим из графика, она меняется со средним периодом Т_W = 68.4 тыс. лет вокруг среднего значения j_Wm = 0.068 радиан. На основной период накладываются колебания с большей длительностью.

Угол наклона плоскости орбиты і к плоскости экватора 1950 г. испытывает колебания с таким же периодом Т_і = 68.4 тыс. лет вокруг среднего значения і_m = 0.402 радиан. Колебания угла i происходит в пределах от 0.36 до 0.45 радиан. Диапазон колебаний составляет 5°.

Рис. 2. Эволюция параметров орбиты Земли за 3.76 млн. лет по результатам численного решения уравнений (3): е - эксцентриситет; j_W - угловое расстояние восходящего узла орбиты в плоскости экватора 1950 г. от точки равноденствия; i - угол наклона плоскости орбиты к плоскости экватора 1950 г.; j_р - угловое расстояние по орбите от восходящего узла до перигелия; Т - время в тысячелетиях от 1950 г.; Т_е, Т_W ,Т_i - наименьшие периоды в тыс. лет колебаний эксцентриситета, восходящего узла и наклона орбиты; Т_р - длительность в тыс. лет одного оборота перигелия; углы даны в радианах.

Угол положения перигелия j_р отсчитывает в плоскости орбиты от ее восходящего узла. Как видно из графика, угол монотонно увеличивается со временем. Перигелий со средней за 3.76 млн. лет угловой скоростью w_р = 970" в столетие перемещается в направлении обращения Земли вокруг Солнца, совершая один оборот за 133.7 тыс. лет.

После работ М. Миланковича последующие исследователи проводили расчеты, уточняя приближенный аналитический метод и исходные данные. Из сопоставления эволюции орбиты Земли, рассчитанной разными авторами видно, что вначале результаты их совпадают, а затем с увеличением времени в несколько сот тысяч лет нарастают рассогласования и при достижении миллиона лет изменения становятся существенными. В качестве иллюстрации этому рассмотрим сопоставления наших данных с данными Ш.Н. Шараф и Н.А. Будниковой [Шараф и Будникова, 1969] (см. рис.3а, б, в), также А. Берже и М.Ф. Лоутре [Berger and Loutre, 1991] (см. рис.3г, д, е) для эксцентриситета е. Последние авторы выполняли свои расчеты спустя два десятилетия. Как видно из графиков, их результаты лучше согласуются с результатами наших численных решений уравнений движения (3).

С дальнейшим увеличением времени сверх 2 млн. лет в прошлое наши результаты начинают отличаться от аналитических результатов А. Берже и М.Ф. Лоутре. Это отличие, естественно, приведет и к отличиям в рассчитанных климатических изменениях. Поэтому, если ставить вопрос о сопоставлении рассчитанных климатических изменений с палеоданными давностью больше 2-х миллионов лет, необходимо быть уверенным, что результаты расчетов не содержат ошибочных изменений. Достоверность численных решений может быть установлена и проверена рассмотренными методами. Что же касается проверки достоверности аналитических решений, то она всегда будет оставаться проблематичной.

Угловые параметры орбиты j_W, j_р и і мы определяем относительно неподвижной плоскости экватора на эпоху 1950 г. В аналитических методах положение перигелия и наклон плоскости орбиты рассчитываются относительно подвижной плоскости экватора, которая в свою очередь, как известно из наблюдений, прецессирует со скоростью около 50" в год в направлении обратном движению Земли. Пересчитаем эволюцию перигелия Земли, представленную на рис. 2, в подвижную систему координат.

Восходящий узел g перемещается против движения Земли в плоскости неподвижной эклиптики с угловой скоростью, отнесенной к одному году [Справочное руководство…, 1976]

где T₀ - время в столетиях, отсчитываемое от фундаментальной эпохи 1900 г

Рис. 3. Сопоставление эволюции эксцентриситета (точки) за 2.5 млн. лет по результатам численных решений уравнений (3) (Т - время в тысячелетиях) с данными: а, б, в - Ш. Г. Шараф . и Н. А Будниковой (1969 г.); г, д, е - А. Берже и М. Ф. Лоутре (1991 г.).

Рис. 4. Сопоставление результатов численного решения (точки) с аналитическими результатами А. Берже и М.Ф. Лоутре (j_рB): а, б - положение перигелия относительно восходящего узла в подвижной плоскости экватора j_рn; в, г - угла наклона орбиты i к неподвижной плоскости экватора с углом наклона e к подвижной плоскости экватора ) (Т - время в тысячелетиях).

Расчет движения Солнечной системы мы проводим в неподвижной экваториальной системе координат, ось x которой направлена на точку весеннего равноденствия g₀ на эпоху 1950 г. Так как этой эпохе соответствует величина времени Т₀ = 0.5, то, согласно (7), угловая скорость будет 5037",3305 угловых секунд в столетие. Тогда в любой другой момент времени Т смещение точки весеннего равноденствия будет определяться дугой gg₀ :

где T - время в юлианских столетиях от 30.12.49 г. Тогда положение перигелия от подвижного восходящего узла будет

Следует отметить, что при выводе (8) мы пренебрегли зависимостью p₁ от времени, так как она теряет смысл при распространении ее на большие сроки.

На рис.4а, б сопоставлены наши изменения перигелия j_pn относительно подвижной плоскости экватора с изменениями j _pB, рассчитанными А. Берже и М.Ф. Лоутре. Результаты расчетов в пределах первого миллиона лет совпадают с высокой точностью, затем начинают расходиться. Для сравнения, на рис. 4а приведены изменения перигелия j_p относительно неподвижной плоскости экватора. Как видим, основная часть изменения перигелия обусловлена движением восходящего узла (gg₀), т. е. прецессией плоскости экватора. Поэтому имеющиеся расхождения j_pn и j _pB обусловлены точностью описания движения плоскости экватора.

На рис.4в, г сопоставлен угол наклона плоскости орбиты і с результатами тех же авторов, но для угла наклона плоскости орбиты e к подвижной плоскости экватора. Величина угла e включает изменения, обусловленные прецессией плоскости экватора. Из графика видно, что в начальный момент (T = 0) изменения углов і и e совпадают. За счет более короткого периода прецессии земной оси, равного 25.7 тыс. лет, изменение угла e не следует за изменением і в течение всего периода Ті = 68.4 тыс. лет, а изменяется на обратное. Как видим, на всем интервале Т градиенты изменения і и e совпадают, т. е. они обусловлены движением плоскости орбиты относительно неподвижного пространства. А на периоды колебаний e и на их амплитуду накладывает влияние движение плоскости экватора.

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

В работе Ж. Ляскара и др. (Laskar at al., 1993) решения получены в результате численного интегрирования упрощенных уравнений движения в рамках теории вековых возмущений. В этой работе также рассматривается численное решение уравнений прецессионного движения. Авторы получили решение за интервал времени -200 млн. лет. Однако решения имеют небольшое схождение и после 100 млн. лет носят хаотический характер. Поэтому Ж. Ляскар и др. ограничились периодом от -20 млн. лет до +10 млн. лет.

В работе Кинна и др. [T.Quinn at all, 1991] численно интегрируются неупрощенные уравнения движения орбитального движения. Авторы применяет конечно-разностный метод Штёрмера 12 и 13 порядков. В обеих этих работах накладываются дополнительные слабые воздействия: релятивистская добавка к закону тяготения Ньютона и приливное трение. Ж. Ляскар и др. сопоставили свои результаты по наклону орбиты e c результатами Кинна и др. Последние имеют небольшое отставание по времени, которое, как доказали Ляскар и др., обусловлено влиянием приливного трения.

Кинн и др. сопоставили в пределах -1 млн. лет свои решения по эксцентриситету с результатами Берже и Лоутре [Berger and Loutre, 1991]. Характер отличия такой же как на рис. 3а. Отсюда следует, что наши результаты по эксцентриситету в пределах -1млн. лет лучше согласуются с численными решениями Кинна и др. чем с аналитическими Берже и Лоутре. Этот вывод качественен. Получение количественных сопоставлений нецелесообразно, так как эти авторы ввели дополнительные воздействия. В связи с этим сопоставим постановки задачи в нашей работе и в статье Кинна и др.

Эти авторы решают уравнения относительного движения, т.е. движения тел по отношению к Солнцу. Мы решаем в барицентрической системе координат и дополнительно получаем также решения для движения Солнца. Во-первых, движение Солнца представляет самостоятельный интерес для короткопериодических изменений климата [Мельников, Смульский и др., 2000]. Во-вторых, его движение позволяет понять ряд особенностей движения планет. В-третьих, в барицентрической системе координат остаются неизменные интегралы движения. Поэтому погрешность координат и скоростей тел легко определяется и в процессе решения, и при многочисленных переработках алгоритма и программы решения уравнений.

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

Рассмотрим теперь погрешность в работе Кинна и др. С этой целью в табл. 1 моменты количества движения Солнечной системы исходной эпохи (30.12.1949 г., юлианский день которой JD₀ = 2433280.5) сопоставлены с моментами в конце счета до времени 3.056 млн. лет назад (JD_k = -11133787075.5). Как видно из таблицы, относительное изменение модуля момента составило величину dM = 1.55·10^-7.

В табл. 2 сопоставлены моменты начальной и конечной эпох, согласно расчетам, представленным в нашей статье. Как видим, относительное изменение момента составляет величину dM = 1,3·10^-10, которая на 3 порядка меньше соответствующей величины Кинна и др [2].

Эту задачу мы решили также с меньшим шагом Dt=10^-4 года. В этом случае погрешность момента dM = 5.5·10^-12, т.е. на 5 порядков меньше чем в работе Кинна и др. Кроме того, мы также выполнили расчет с расширенной длиной числа (34 десятичных знака) и Dt = 10^-5 года на период от -2.5 тыс. лет до 2.5 тыс. лет. Погрешность момента растет линейно со скоростью dM/dt = 1.48·10^-23 1/год. Это даст погрешность момента на эпоху -3,056 млн. лет dM= 4.4·10^-17, что на 10 порядков меньше, чем в работе Кинна и др.

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

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

2. В задаче интегрируются уравнения движения Луны. Это позволяет, во-первых, более точно определить изменения орбиты Земли. Во-вторых, полученное движение Луны позволит более точно рассчитать прецессию оси Земли.

3. Задача решается более совершенным методом, который позволяет интегрировать уравнение движения Луны. В статье представлены первые результаты за -2.56 млн. лет - +1.2 млн. лет. Метод позволяет достичь периодов сотен млн. лет, которые предыдущими исследователями, заявлены как предел стабильности Солнечной системы

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

1. Разработана методика численного интегрирования уравнений, которая позволяет с приемлемой точностью рассчитывать эволюцию Солнечной системы за сотни миллионов лет.

2. Рассчитана эволюция орбиты Земли за несколько миллионов лет:

2.1. Эксцентриситет орбиты изменяется в пределах е = 0.0022 - 0.062. Наименьший период колебаний равен 94 тыс.лет. А перигелий за время 133.7 тыс. лет совершает один оборот в направлении движения Земли.

2.2. Плоскость орбиты относительно фиксированной плоскости экватора колеблется с наименьшим периодом 68.4 тыс. лет и размахом колебаний 5°.

3. Независимым методом подтверждены результаты предшествующих исследователей об эволюции орбиты Земли и о возможном ее влиянии на криосферу Земли.

Авторы считают необходимым выразить свою признательность заведующему отделом небесной механики ГАИШ МГУ, доктору физ.-мат. наук Н.В. Емельянову за детальный анализ рукописи, полезные замечания и ценную информацию; а также благодарят руководство и сотрудников Сибирского суперкомпьютерного центра СО РАН за безотказную помощь и консультации по работе на суперкомпьютерах.

1. Мельников В.П., Смульский И.И., Кротов О.И., Смульский Л.И. Орбиты Земли и Солнца и возможные воздействия на криосферу Земли (постановка проблемы и первые результаты)// Криосфера Земли. - 2000, т. IV, №3, с. 3-13.
    2. Миланкович М. Математическая климатология и астрономическая теория колебаний климата. - М.-Л. -ГОНТИ. - 1939. -207 с.
    3. Смульский И.И. Теория взаимодействия. - Новосибирск: Из-во Новосиб. ун-та, НИЦ ОИГГМ СО РАН, 1999 г. - 294 с.
    4. Справочное руководство по небесной механике и астродинамике / Под ред. Г. Н. Дубошина. Изд. 2-е, доп. и перераб. М., Наука, 1976, 862 с.
    6. Сулержицкий Л.Д., Романенко Ф.А. Возраст и расселение "мамонтовой" фауны азиатского Заполярья (по радиоуглеродным данным)// Криосфера Земли. - 1997, т. I, №4, с. 12-19.
    5. Шараф Ш. Г. и Будникова Н. А. О вековых изменениях элементов орбиты Земли, влияющих на климаты геологического прошлого. // Бюл. ИТА АН СССР. 1967, вып. 11, № 4 - С. 231 - 261.
    6. Шараф Ш. Г. и Будникова Н. А. Вековые изменения элементов орбиты Земли и астрономическая теория колебаний климата // Тр. Инст. теоретич. астрономии. - Вып. XIV. - Л.: Наука. - 1969 г. - с. 48 - 109.
    7. Berger A. and Loutre M. F. Insolation values for the climate of the last 10 million years// Quaternary Science Reviews. - 1991. - № 10, Р. 297 - 317.
    8. Brouwer D., Van Woerkom A. J. J. The secular variation of the orbital elements of the principal planets // Astr. Pap. - 1950. - 13, 2.
    9. Laskar J., Joutel F., Boudin F. (1993) Orbital, precessional, and insolation quantities for the Earth from -20 Myr to +10 Myr. Astronomy and Astrophysics. V. 270. N. 1-2. P. 522-533.
    10. Quinn T. R., Tremaine S., Duncan M. (1991) A three million year integration of the earth's orbit. Astronomical Journal. V. 101. P. 2287-2305.