Математическая модель солнечной системы

И.И. Смульский

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

Рис. 1. Параметры орбиты планеты в неподвижной экваториальной гелиоцентрической системе координат Sxyz (а) и в плоскости орбиты Sx_oy_o (б):

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

Отметим, что вышеупомянутая модель Солнечной системы (назовем ее моделью Кеплера), в отличие от предшествующих моделей, представляет движение планет вокруг Солнца, в то время, как в более ранних моделях движение небесных светил было отнесено к небосводу и представляло движение светил вокруг Земли. Одна из самых знаменитых, модель Птолемея, служила человечеству более 14 веков. Эта модель описана в книгах Клавдия Птолемея [1], однако известно [2], что она основана на более ранних работах Александрийской и Греческой астрономических школ (к которой, например, относится и великий древнегреческий астроном Гиппарх), на трудах вавилонских астрономов и, по-видимому, ученых древнего государства Шумер. Кроме птолемеевской системы, древним народам были известны и другие модели Солнечной системы [2], например, одна из них - это Дагестанская модель, в которой движение Солнца и Луны, а возможно и других светил, моделировалось перемещением метки по лункам, выдолбленным на камне. М.И. Исрапилов в своей книге [3] приводит много вариантов таких моделей, возраст которых исчисляется десятками тысяч лет.

Будем рассматривать движение планеты с массой m при ньютоновском воздействии на нее со стороны Солнца, масса которого равна М. Дифференциальные уравнения ее движения относительно Солнца в безразмерных величинах имеют вид:

где – относительный радиус положения планеты относительно Солнца;

_G(M+m) – параметр взаимодействия;

Из уравнения (1) видно, что взаимодействие двух тел, а также их движение определяется только параметром α₁, который мы назовем параметром траектории [4]. Его использование вместо эксцентриситета e оправдано тем, что уравнения взаимодействия и их решение приобретают общность как в случае сил, зависящих только от расстояния, так и в случае сил, которые зависят также и от скорости [5]. В таком виде оказалось удобным рассмотрение задачи многих тел, «осесимметрично» расположенных на плоскости [6].

Уравнение траектории тела в безразмерых переменных в полярной системе координат (рис. 1б), с началом в Солнце, имеет вид [4]:

, (2)

где угол φ отсчитывается от отрезка SB = .

Уравнение (2) при ₁ = -1 представляет окружность, при -1 < ₁ < -0.5 – эллипс, при ₁=-0.5 – параболу, при -0,5<₁<0 – гиперболу, а при ₁=0 – прямую.

Легко показать, что в безразмерных величинах радиус афелия R_a, скорость в афелии v_a, время движения t_a от перигелия до афелия и период обращения T_tr определяются следующими равенствами:

Кроме того, запишем выражения для радиальной и трансверсальной скоростей

а также для эксцентриситета е орбиты и ее большой полуоси a,

В результате интегрирования уравнений (1) мы получили [4], что величина промежутка «безразмерного времени» движения планеты по эллиптической орбите от точки перигелия до любого положения на орбите с радиусом может быть вычислена по формуле:

Если известны в момент времени T координаты (x,y,z) и скорости (v_x,v_y,v_z) планеты в неподвижной гелиоцентрической экваториальной системе координат Sxyz ( рис.1 а), можно определить ее положение на орбите в полярной системе координат r, ( рис. 1б). С этой целью сначала найдем модули радиуса и скорости планеты

Угол β между радиусом планеты и вектором ее скорости (рис. 1б) определяется из выражения

Тогда компоненты скорости в плоскости орбиты равны

Значение угла β можно также выразить с помощью выражений (4)

Используя формулы (8) и (10), можно найти величину , и, далее, радиус (расстояние) перигелия, который будет равен

Верхний индекс “d” указывает на то, что величина R_p вычислена с использованием декартовых координат. Если известен параметр орбиты ₁, тогда формула (11) дает возможность определить радиус перигелия R_p^d.

Вместе с тем, очевидно, что реальная орбита планеты отличается от эллиптической. Декартовые координаты планеты x, y, z и компоненты ее скорости v_x, v_y, v_z могут принадлежать точке орбиты, которая «несколько» удалена от усредняющего ее эллипса, и тогда, определенный по формуле (11) радиус (расстояние) перигелия может соответствовать орбите, отличающейся от этого эллипса. В результате численных экспериментов было установлено, что лучше использовать значение радиуса перигелия R_p, который определен для эллипса, «усредняющего» данную орбиту, но при его отсутствии можно воспользоваться и выражением (11).

Если известно значение R_p, тогда находим относительный радиус и из уравнения траектории (2) определяем полярный угол  планеты по формуле

Скорость в перигелии по декартовым координатам с учетом (4) определяется с помощью трансверсальной скорости (9)

а по известному периоду обращения T_tr ,более точно , с учетом соотношения (3), она запишется в виде

Далее, по формуле (6) можно определить время движения от перигелия до точки орбиты с радиусом :

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

Таким образом, используя известные декартовы компоненты радиуса и скорости , по формулам (12) и (15) можем определить угол положения и время движения t_p от точки перигелия.

3. Определение положения планеты на орбите в новый момент времени.

3.1. Алгоритм . Прогнозирование величины радиуса.

Найдем теперь положение планеты в новый момент времени, T₁=T+h. Так как выражение (6) для момента времени – трансцендентно и не выражается в явной компактной форме через радиус , будем решать эту задачу методом последовательных приближений, для определения которого потребуется радиальные и угловые проекции ее ускорения. После дифференцирования по времени величины , согласно (4), получаем величину радиального ускорения

а после дифференцирования угловой скорости по времени получим угловое ускорение:

Итерационный процесс целесообразно реализовать в двух вариантах: по радиусу r и по полярному углу φ. В первом варианте прогнозируемое значение радиуса в новый момент времени T₁ определим с учетом второй производной

а прогнозируемое значение угла φ - с учетом только первой производной:

где h_i – значение приращения времени в i-той итерации, которое в первой итерации i= 1 принимается равным h₁=h.

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

. (20)

Далее можно определить отклонение между этой величиной и заданным приращением времени h

и ее относительное значение по очевидной формуле

Пусть точность вычислений характеризуется величиной EPS. Если , то вычисляем новое приращение времени

и повторяем процесс итераций, начиная с формулы (18). Итерационный процесс прекращается при достижении . Аналогично вычисляем величину φ₁. Если φ₁>2π, то сначала вычисляется количество прохождений планеты через перигелий по формуле K_tr'=K_tr+1. Затем угол φ₁ вычисляется более точно по формуле (12) с учетом «окончательного» значения r₁. Таким образом, при завершении итерационного цикла мы определяем для нового момента времени T₁ полярные координаты φ₁ и r₁ планеты в плоскости ее орбиты с заранее заданной точностью EPS.

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

Согласно (4), в момент T в точке орбиты с радиусом , мы вычисляем трансверсальную скорость

, (24)

а также квадрат и первую степень радиальной скорости

Если , то принимается . Когда планета находится в нижней полуплоскости (π<φ<2π) (рис. 1б), ее радиальная скорость – отрицательна, т.е. . Затем по формуле (17) определяется угловое ускорение ε.

В новый момент времени T₁ мы прогнозируем значение угла φ₁ с учетом второй производной:

. (26)

Если в выражение для времени (6) подставим значение величины из уравнения (2) , тогда получим выражение для времени движения планеты от перигелия до точки с полярным углом по формуле

и далее, по этой же формуле находим время движения планеты от перигелия до положения с новым углом . Затем, используя формулу, аналогичную выражению (20) , находим разность в значениях времени

Далее, необходимо повторить алгоритм, описанный соотношениями (22), (23) , включая сюда и расчет нового приращения времени

до тех пор, пока относительное значение невязки не станет меньше EPS.

Таким образом, в этом случае итерационный цикл начинается с формулы (24) и заканчивается формулой (29). При этом в момент времени T₁, с точностью порядка EPS, мы определяем полярный угол планеты и ее время движения от момента прохождения через перигелий . Полярный радиус r₁ находится из уравнения траектории (2) в зависимости от значения , а по формулам (4) находим радиальную скорость , а также ее трансверсальную составляющую и полную скорость в зависимости от величины . При количество прохождений K_tr через перигелий увеличивается на единицу, а угол пересчитывается, согласно формуле .

4. Дополнительные условия.

Из-за того, что время прохождения планетой нижней и верхней частей эллиптической орбиты (рис. 1б) определяется одними и теми же зависимостями (6) и (27), а угол является циклической переменной, в «апсидных» точках с радиусами R_p и R_a возникают особенности, поэтому необходимо использовать некоторые дополнительные условия. Рассмотрим основные из них.

Из рис. 1б видно, что при нахождении планеты в нижней части эллипса (т. P₂) угол β>90º, а радиальная скорость v_r отрицательна (v_r<0 ). В этих точках ,по формулам (6) и (27) определяется время движения планеты на участке BP₂, а не на участке BP₁P₂. Угол для участка BP₂ также вычисляется по формуле (12), поэтому после определения угла и с учетом условия cosβ<0, найденное значение необходимо заменить значением . Далее пересчитывается время движения до перигелия, используя формулу (15). Иными словами, имеем, что, если cosβ<0, тогда t_p'=2t_a-t_p. В дальнейшем вычислительном процессе новые значения величины осуществляются с учетом условий и . Радиальная скорость v_r также меняется свой знак из-за этого условия. Подобные вычисления проводятся также и при использовании второго алгоритма.

Нестандартная ситуация возникает в нижней полуплоскости, когда шаг по времени h превосходит промежуток времени движения до перигелия , т.е. при условии, что . В этом случае положение планеты в момент T «приводится» к точке перигелия r=R_p, и тогда все величины принимают свойственные перигелию значения: ; , . Приращение времени задается в виде , а абсолютная и относительная невязка вычисляются по формулам

Аналогичная ситуация возникает и при прохождении противоположной точки орбиты, афелия.

При реализации второго алгоритма ситуация в перигелии зависит от неравенства . В этом случае, после вычисления разности времени по формуле (28), она пересчитывается согласно соотношению .

Вблизи точек перигелия и афелия полярный радиус r₁ изменяется слабо, поэтому приближенное его значение, вычисленное по формуле (18), может оказаться меньшим R_p или большим R_a. В этом случае величина r₁ пересчитывается по формулам:

при r₁<R_p , r₁'=0.5(r+R_p), при r₁>R_a , r₁'=0.5(r+R_a). (32)

Для второго алгоритма такая ситуация не возникает, т.к. значение радиуса r₁ вычисляется из уравнения траектории (6).

Зная полярные координаты планеты r,φ в плоскости ее орбиты, необходимо определить ее декартовы координаты в неподвижной системе координат Sxyz (рис. 1а). Для этого спроектируем координаты планеты в плоскости ее орбиты (см. рис 1б) на взаимно перпендикулярные прямые:

Отрезок SL_o (см. рис. 1а) спроектируем на ось z и на проходящую в плоскости Sxy прямую SL

Итак, в плоскости Sхy координаты планеты в системе координат с осями SK и SL, выражаются формулами

С учетом (33) и (35) полярные координаты планеты r,φ на подвижной плоскости орбиты преобразуются в ее «неподвижные» декартовые координаты по формулам:

Реальные параметры орбит планет подвержены, в частности, вековым изменениям, т.е. их параметры , и i(t), e(t), R_p(t) и T_tr(t) являются функциями времени. Поэтому, если на каждом шаге вычислений полярные координаты планеты преобразовывать в декартовые координаты, тогда последние будут представлять движение планеты практически с той же точностью, что и вековые изменения параметров орбиты.

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

6. Начальные данные.

Под начальными данными здесь мы понимаем всю совокупность параметров, которые необходимы для реализации описанных алгоритмов. Для задачи двух тел (2) – (6) начальными данными являются параметр орбиты α₁ или эксцентриситет e, радиус перигелия R_p, период обращения T_tr и момент T₀ прохождения через перигелий или другую точку орбиты, которая характеризуется полярным углом . Эти параметры могут быть определены по начальным данным планеты различными способами. В результате проведения вычислений и сопоставлений мы остановились на следующих начальных данных планеты. В момент T₀ в гелиоцентрической экваториальной системе координат задаются компоненты координат планеты x, y, z и скорости v_x, v_y, v_z. Кроме того, задается 6 параметров орбиты i(t), , , (или ), R_p(t) и T_tr(t). Эти параметры мы получили в результате численного интегрирования уравнений движения Солнечной системы за разные отрезки времени, в том числе за 100 млн. лет. На интервале в нескольких тысячелетий их можно также получить из вековых возмущений С. Ньюкомба [7], которые приведены также на стр. 487 [8], или из средних элементов орбит в работе Дж.Л. Симона с соавторами [9].

Могут быть упрощения при задании параметров орбиты. Так как полуось орбиты a и период обращения T_tr испытывают небольшие колебания, то ими в этом варианте можно пренебречь и считать их постоянными, а параметры i, , , e – линейными функциями времени. Радиус перигелия R_p выражается через полуось a по формуле . С такими значениями параметров можно вычислить положение планет в пределах нескольких столетий.

Если пренебречь зависимостями от времени всех параметров i, , , e, R_p и T_tr, то тогда можно рассчитать движения планет в интервалах времени в несколько лет.

Рассматриваемая модель была разработана для исследования эволюции вращательного движения Земли и была реализована в компьютерной «среде» MathCad, а также на языке программирования Фортран. Она является составной частью общей программы интегрирования уравнений вращательного движения планеты [11].

В вышеприведенном алгоритме все упомянутые координаты и параметры для каждой планеты свои, т.е. должны индексироваться по планетам. Только время T и приращение времени h являются общими. Начальные положения девяти планет Солнечной системы, от Меркурия до Плутона, задавались на момент 1949, декабрь 30.0 ET в гелиоцентрической экваториальной системе координат эпохи 1950.0.

Мы исследовали, в частности, реальную эффективность итерационного процесса. Орбиты больших планет имеют небольшие эксцентриситеты, поэтому для них алгоритм №2 оказался более эффективным, то есть вычислительный процесс лучше сходится нежели в алгоритме №1. В табл. 1 представлены рассчитанные по формуле (22) относительные погрешности времени δh_i после первых трех итераций для всех планет. Из таблицы видно, что после каждой итерации относительная погрешность во времени уменьшается не менее чем на три порядка.

Таблица 1. Относительные погрешности времени δh_i после первых трех итераций для планет от Меркурия (1) до Плутона (9) при счете по алгоритму 2 математической модели Солнечной системы.

В алгоритме №1 хорошо «прогнозируется» значение полярного радиуса, поэтому он будет иметь меньшие погрешности для сильно вытянутых орбит, например, для кометных орбит.

Для проверки точности расчета по данной модели движения планет были проинтегрированы дифференциальные уравнения движения Солнца, планет и Луны [10] и получены их законы движения, т.е. зависимости от времени их координат: x(t), y(t), z(t) за одно обращение планеты вокруг Солнца. Такие же «законы движения» этих планет были получены для математической модели Солнечной системы. Расхождение между этими «законами» очень мало.

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

где x_E, y_E и z_E выражают «закон движения» Земли. На рис. 2 представлены проекции орбит планет относительно Земли на экваториальную плоскость, т.е. y_r = f(x_r). Орбита Земли дана относительно Солнца. Другими словами, на рис. 2 приведено также сравнение проекциий орбит планет и Солнца при их движении относительно Земли. Результаты численного интегрирования изображены жирной пунктирной линией, а рассчитанные по нашей модели орбиты – тонкой сплошной. В вычислениях использовался нами шаг интегрирования h = 1·10^-4 года. Интервал времени изменялся от 0.24 года для Меркурия до 250 лет для Плутона. Как видно из графиков на рис. 2, орбиты планет и Солнца относительно Земли, рассчитанные двумя методами, фактически не различаются.

Рис. 2. Сравнение проекций орбит планет на плоскость экватора x_ry_r, рассчитанных двумя способами: 1 – по результатам численного интегрирования дифференциальных уравнений движения тел Солнечной системы; : 2 – по математической модели Солнечной системы. Орбита Земли приведена относительно Солнца, а планет – относительно Земли; время счета равно периоду обращения каждой планеты; координаты x и y выражены в астрономических единицах; для Плутона графики приведены в двух видах: а – в обычном; б – в увеличенном.

Как уже упоминалось ранее, по этой модели можно рассчитывать движение малых планет, астероидов, комет и спутников планет. Однако, для спутников, в отличие от других тел, все орбитальные расчеты необходимо выполнять в планетоцентрической системе. Еще раз отметим, что для расчета движения всех этих тел необходимо иметь начальные их координаты и вековые изменения параметров орбит: _р(t), _(t), i(t), e(t), R_p(t), T_tr(t).

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

2. Вычисление движения планет Солнечной системы и сравнение его с движением, рассчитанным другим методом, говорит в пользу достоверности модели, полученной более эффективными и экономичными методами.

В заключение автор выражает признательность профессору Е.А. Гребеникову за ряд существенных корректировок текста статьи, формулировку цели статьи и ее выводов.

Литература

2. Лаплас П.С. Изложение системы мира. – Л.: Наука, 1982. - 376с.

3. Исрапилов М.И. Наскальные рисунки Дагестана и колебания полюсов и наклона оси Земли в голоцене / Махачкала: Изд-во «Юпитер», 2003г. - 432с.

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

5. Смульский И.И. Траектории при взаимодействии двух тел, зависящем от относительного расстояния и скорости//Математическое моделирование. - 1995. - Т.7. - N7. - С.117-126.

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

7. Newcomb S. The elements of the four inner planets and the fundamental constants of astronomy. Washington: Government printing office. 1895. –202 p.

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

9. 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.

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