Метод симплектического интегрирования уравнений движения малых тел Солнечной системы Текст научной статьи по специальности «Математика»

УДК 521.52,17

МЕТОД СИМПЛЕКТИЧЕСКОГО ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ ДВИЖЕНИЯ МАЛЫХ ТЕЛ СОЛНЕЧНОЙ СИСТЕМЫ

ЕЛ. Чегодаева

В статье предложен метод симплектичекого интегрирования. На основе его создан симплектические интеграторы разных порядков и проведено их тестирование.

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

Рассмотрим основы теории симплектических интеграторов для канонических уравнений движения с гамильтонианом Н для системы N тел [1, 12].

с1х} _ дН йрх _ дН ^

Л др1 ' Л дх1 '

где / = х ~(хих27х3) - обобщенная координата, р~(рир2,р3) - обобщенный импульс.

Используя (1), скорость изменения любой динамической величины д{х,р,/) может быть записана в виде

с1д ^ ( дд йхг дд с1рг

N

дд дН дд 8Н

& ТГДйх, Л др1 (31 где {,} - скобки Пуассона, ^ - дифференциальный оператор вида

N

* ( д дН д дНЛ

^ чйх, дрг дрг дхг Общее решение уравнения (2) имеет вид

(3)

т = * (4)

где д^-т) значение д в предыдущий момент времени.

Предположим, что Я = Нл+еНВ9 где НА и Нв являются интегрируемыми гамильтонианами. Симплектический интегратор второго порядка может быть записан в виде

д(()^етететлет£в12д^-т), (5)

где г - шаг интегрирования, а А и В ~ дифференциальные операторы вида ^ для гамильтонианов Нл и Нв соответственно.

Каждый шаг интегрирования для схемы (2) состоит из трех подшагов.

1) полшага решение канонических уравнений движения с гамильтонианом Нв.

2) шаг решение канонических уравнений движения с гамильтонианом НА .

3) полшага решение канонических уравнений движения с гамильтонианом Нв. Применение симплектических интеграторов равносильно точному решению уравнений движения (1) с гамильтонианом Нтщт [2]:

Яю1еёГ=Я + Яе1Т(г'), (6)

где Яегг - член, возникающий из-за некоммутативности операторов А и В, р ~ порядок интегратора. Интеграторы порядка выше второго намного более громоздкие, а так же имеют неудобства в виде отрицательных подшагов, поэтому чаще всего используются интеграторы второго порядка как достаточно эффективные и быстрые. Однако есть другие способы повысить точность симплектического интегратора - интеграторы псевдовысокого порядка [2, 5] и корректоры [7, 9, 11]. Интеграторы псевдовысокого порядка получают из 32(т), разделяя экспоненты специальным образом [5]. Например,

(7)

2_ 1' 75' * 2

s6(T) = V2rVirV2rV2rV>rS,

л/5 J 12 2 12

cl -~>С2

Интегрируемый гамильтониан для этой схемы Я,теёг = Я + е2т2 (132^) {{НА, Нв}, Нв} + 0(ете).

Ласкар и Робутель [5] представляют общий вид коэффициентов через полиномы Бернулли, позволяя найти схему интегратора для любого псевдовысокого порядка. Эти интеграторы тратят больше компьютерного времени, но имеют более высокую точность, поэтому можно увеличить шаг интегрирования. Чемберс и Морисон [2] считают, что интегратор псевдошестого порядка наиболее эффективен.

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

Л-1 N тпг

7 (8)

1=0 ZЩг

(=0 у=/+1 I Яг ~Я} I

где дг - обобщенные координаты, рг - обобщенные импульсы, тг - масса г-го тела, г = 0 соответствует Солнцу. Висдом и Холман [8] в своем интеграторе использовали координаты Якоби. Можно также использовать гелиоцентрическую систему координат, барицентрическую [4], смешанную [3]. Левисон, Дункан, Ли [3] предложили использовать смешанную систему координат. Такая система координат удобна для исследования объектов имеющих тесные сближения. Координаты <2{ - гелиоцентрические, а импульсы Р1 - барицентрические:

Q,

Яг ~#о>еСЛИ 1 *

1И Р

— V /я, <7,, если / = 0, 1

mi - л

М 7=0

7-0

(9)

где Л^ХНо^' Тогда (17) имеет вид #02,/)) - ЯКер1 + Я3ип н-Я1тег, где

Я

п (\т> (К у

Kepi

i=l

п-1

m

о J

О?п1т0

1ЙТ

я,

Бип

I I Щ.

т0 /=1 J=\,J*I

/7-1

Я

GmlmJ

Inter X X I п - П 1=1 у=/+1 I И И j

(10)

(11)

(12)

Для такого гамильтониана, разделенного на три части, схема симплектического интегрирования второго порядка имеет вид [3]

= (13)

где А,В,С - операторы соответственно для ЯКер|,Я1тег,Я5ип . При симплектическом интегрировании шаг интегрирования по независимой переменной является постоянным. Так в симплекти-ческом интеграторе Висдома и Холмана [8] временной шаг интегрирования является постоян-

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

Интегратор в смешанных координатах для кометных орбит, разработанный в нашей работе, основан на методе Левисона, Дункана и Ли [3], упомянутом выше, и методе переменного временного шага Емельяненко [4].

В программе производится симплектическое интегрирования отдельно для планет и отдельно для частиц, что экономит компьютерное время. Система смешанных координат, гамильтониан и схема интегрирования для планет описываются формулами (9)-(13). Симплектическое интегрирование планет идет с постоянным временным шагом.

Найдем гамильтониан для частицы в смешанных координатах. Для этого перепишем (8) в виде

* п 2 п2 N-1 N т т N

7=0 ¿щ 1=0 I чу I 1=0 I Ят Чх I

где т - индекс частицы. Переходим к смешанным координатам, используя (9) для планет, а для частицы

<2т = Чт-Яо> (15)

?т = ттГт = ттУт~ — ^Р1- (16)

Центр масс (/ = 0) движется как свободная частица, в гамильтониане его не учитываем. Выпишем из гамильтониана системы части, относящиеся к частице:

(17)

2 2т0 i=0|ßm-ß

т,

о /=/

На тт можно сократить, после чего считаем тш = 0 . В результате чего получаем гамильтониан

Н -НКер1 -#Sun -tfInter

(18)

Kepi 2 ,

1 " #Sun =--VmEßiit),

(19)

т0 i=1

" lr2m

Нш = У _КЩ-. (20)

tt\Qm-QM

Для данного гамильтониана применим преобразования Емельяненко [4]. Расширяем фазовое пространство, введем канонические переменные q4~t и р4 =-Н . Затем, следуя Микколе [6], определим новую независимую переменную г соотношением, dt = dt/r . Тогда уравнения движения имеют вид:

^ = = (21) dt dpj' dr dqi

где / = 1,2,3,4, а К = r(#Kepl + р4)~ г(#Inter - Я5ш). Добавим выражение гВ0+Вх> 0 к обеим частям гамильтониана К (BQ9Bi - некие малые константы). Тогда К = К0 -Кх, где

= '-(якер, + Р4 + В0 + Уг),К] = r(tfInter + tfSun + В0+ Y)' (22)

и - 0 на траектории движения. Далее применяем преобразование, предложенное Микколой и Таникавой [6],

ds = Kxdr = K0dt =—dt9 (23)

г

после которого уравнения движения можно записать в следующем виде:

dQm ЗГ df:

дГ

dQ

т

дТ

с новым гамильтонианом

ds dfm ds

dq4 _ дТ dp4 ds ф4' ds dq4

г = г0+г„

(24)

(25)

(26)

где Г0 = 1п^-)Г1=-1п^ r r

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

Ех{к5/2)ЕМЕх(к5/2\ (27)

где Е0 и Ех - операторы эволюции Г0 и Г{ соответственно, к3 - длина шага для новой переменной 5 . Во время симплектического интегрирования Г0 « Г\ и, соответственно, К0 « Кх. Тогда уравнение (25) имеет вид

ds

dt

К, К

kepi

Р4+В0 +

В\ '

Для возмущенной части решается система дифференциальных уравнений

dQm 5Г,

ds дут Ст0 j

dym_ агг

ds

dt ds

¿i \йт-Ш

дР4

В

Q,

1 • ' «3

ej ,

0 =>t = const.

dp4 дГ}

ds

где С = -у- константа

dt

mn dt ^ \q-qM dt

i=i

ы

(28)

(29)

(30)

(31)

(32)

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

1) решаем уравнения для кеплеровой части планет с шагом по времени 81~ТраП ~Тр1, где

Трап ~ вРемя для которого должны быть найдены новые координаты планет, Тр1 - время в которое координаты уже найдены,

2) зная шаг интегрирования планет &, находим АТ~ТрггХ -Гр1 + Д/72 и решаем уравнения

движения для гамильтонианов (12)и(11)с шагом по времени А Т, находя таким образом нужные нам координаты и импульсы планет:

dQ(t) 1

dt

m

Ё Pi+Pt

0 j=\,j*i

Kmi

1 1

— + —

Щ

dptf) m*mAQt-Q? ¿Чад,

dt

(33)

(34)

QrQ

Q,

Для возмущенной части (29) - (32) выражения получаются громоздкими, на которые тратится довольно много компьютерного времени - это один из минусов смешанного интегратора. Второй минус интегратора в смешанных координатах - это довольно большое изменение части гамильтониана Я5ип. В некоторых случаях выражение под знаком логарифма может стать даже отрицательным. Чтобы этого избежать, надо определить интервал изменения Я8ип и подобрать коэффициенты В0,В1. Принцип выбора параметров В0,В1 следующий: если частица далеко от планет, то

\И\«В0+^-,\П\>В0+^ (35)

г г

при тесных сближениях с планетами. Значит, если взять большое значение В^+Вх/г. то это скомпенсирует изменение ТУ5ип, шаг интегрирования по времени при тесных сближениях будет меняться незначительно. Это означает, что для полученного интегратора нужно очень тщательным образом контролировать параметры интегратора для каждого класса орбит. В этом состоит главный недостаток интегратора в смешанных координатах.

Для того чтобы разработать интегратор с большей областью применения будем использовать барицентрическую систему координат для частиц. Для интегрирования уравнений движения планет будем использовать смешанную систему координат (9)-(11). В данном случае эта система координат наиболее удобна. Для частиц опять используется метод симплектического интегрирования с переменным временным шагом Емельяненко [3]. Гамильтониан имеет следующий вид Я = Я0-Я1, где

#0 = --(36)

2 г

1т1

( 1 1Л

,=0 \\4т~Чх^)\ Г)

(37)

где г =| дт |.

Для решения кеплеровой части переходим от сЬ к Л и решаем задачу двух тел с новым гамильтонианом

Г0=^--^ + В0 + р4, (38)

2 г

где С0 константа К0.

Численные эксперименты

Для тестирования интеграторов были взяты две орбиты с большими полуосями ах - 60,91 а.е. и а2 =97,12 а.е., перигилийными расстояниями дх =30,01 а.е. и д2 =30,2 а.е., эксцентриситетами ех - 0,50, е2- 0,68, наклонами орбит \ = 16°, ¿2 = 13°. Частицы в начальный момент времени находятся на расстоянии 91 и 96 а.е. от Солнца соответственно. В интегрировании участвуют только 4 внешние планеты: Юпитер, Сатурн, Уран и Нептун. Шаг интегрирования для планет равен 365 дней.

Тестированию были подвергнуты симплектический интегратор второго порядка, комбинирующий смешанные и барицентрические координаты, описанный выше и интегратор псевдошес-того порядка. Для построения последнего использовался интегратор, комбинирующий смешанные и барицентрические координаты. Схема интегратора псевдошестого порядка задана уравнением (7). Интегратор основан на уравнениях (36)-(38).Для обеих частиц В0,В1 подобраны на основе соотношений (35): В0 ~0,ВХ =0,21-10"6 . Дяя интегратора второго порядка начальные шаги по времени для частиц составляли 298 и 310 дней соответственно, для интегратора псевдошестого порядка начальные шаги интегрирования равны 910 дней и 953 дня соответственно. Тестирование интеграторов было проведено на интервале 100 млн лет. Время, затраченное процессором на интегратор второго порядка примерно на 5% меньше времени затраченного на интегратор псевдошестого порядка.

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

дн

4 0*10 6 0*10

t лет

3 0x10 1 Ох 10

80

'>10 3x10 4x10 ЭХ10

дН

4 0x10' 6 0x10

1 лет

00 2 0* 10 4 0* 0 6 0*10 а 0x10 1 0x10®

I лет

<и

П!

та 70

2 0x10 4 0x10 5 0x10 3 0x10 1 0x10

t лет

Зависимости ошибок АН и большой полуоси а от времени для первой и второй частицы симплектического интегратора второго порядка и симплектического интегратора псевдошестого порядка

таким образом, чтобы перигелийное расстояние частиц удовлетворяло условию q> 26 а.е. Для орбит с q <26 а.е. должны быть другие коэффициенты В0уВ}. Как видно на рисунке, точность интегратора второго порядка сохраняется на протяжении всего времени интегрирования частиц. Для псевдошестого порядка видно, особенно для второй частицы, накопление ошибок в гамильтониане. Очевидно, это связано с очень большим шагом интегрирования. Известно, что шаг интегрирования не может превышать 400 дней, ограничение, связанное с периодом обращения Юпитера. Кроме того, интегратор псевдошестого порядка требует больше компьютерного времени, чем интегратор второго порядка.

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

Данная работа была поддержана грантами РФФИ-Урал (04-02-96042) и ИНТАС (00-240%

Литература

1. Chambers J JE. A hybrid symplectic integrator that permits close encounters between massivebod-ies// Mon. Not. Roy. Astron. Soc. - 1999. - V. 304. - P. 793-799.

2. Chambers J.E., Murison M.A. Pseudo-high-order symplectic integrators// The Astronomical Journal. - 2000. - V. 119. - P. 425-433.

3. Duncan M., Levison H. and Lee M.A multiple time step symplectic algorithm for integrating close encounters// The Astronomical Journal. - 1998. - P. 2067-2077.

4. Emel'yanenko V. An explicit symplectic integrator for cometary orbits// Celestial Mechanics and Dynamical Astronomy. - 2001. - V. 74. - P. 287-295.

5. Laskar J., Robutel P. High order symplectic integrators for perturbed hamiltonian systems// Celestial Mechanics and Dynamical Astronomy. - 2001. - V. 80. - P. 39-62.

6. Mikkola S., Tanikawa K. Explicit symplectic algorithms for time-transformed Hamiltonians // Celestial Mechanics and Dynamical Astronomy. - 1999. - P. 287-295.

7. Mikkola S., Palmer H. Simple derivation of symplectic integrators with first order correctors// Celestial Mechanics and Dynamical Astronomy. - 2000. - V. 77. - P. 305-307.

8. Wisdom J., Holman M. Symplectic maps for the N-body problem// The Astronomical Journal. -1991.-P. 1528-1538.

9. Wisdom J., Holman M., Touma J. Symplectic correctors. Integration algorithms and classical mechanics// Fields Instituts Community. - 1996. - P. 217-244.

10. Wu X., Huang T.Y., Wan X.S. On correctors of symplectic integrators// Chinese Astronomy and Astrophisics. - 2003. - V. 27. - P. 114-125.

11. Wu X., Huang T.Y., Zhang H., Wan X.S. A note on the algorithm of symplectic integrators// Astrophisics and Space Science. -2003. -V. 283. - P. 53-65.

12.Yohida H. Construction of higher order symplectic integrators// Physics letters A. - 1990. -V. 150.-P. 262-268.

Поступила в редакцию 10 сентября 2004 г.