Задача двух тел. Плутон и Харон

[Tarakan]

Здравствуйте.
У меня проблема с моделированием процессов вращения небесных тел.
Я создаю программу планетарий, в которой движение планет определяется итерациями. На каждой итерации я высчитываю по закону ньютона все силы взаимодействия и затем рассчитываю новое положение тел.
И для задачи двух тел и несоизмеримыми массами всё работает на отлично - Меркурий, Венера, Земля... летают как надо (хоть Меркурий и без процессии). Но вот движение Плутона и Харона у меня не получается вычислить.

Моя программа за основу берёт Кеплеровские элементы орбиты и превращает их в положение тела. Затем считают то-же самое через секунду. Из первой итерации получается радиус-вектор положения из второй итерации - скорость (V=S/t).
После чего начинается цикл просчёта.
За каждый цикл (5 секунд для примера) я считаю расстояние между телами, потом рассчитываю силу взаимодействия по закону Ньютона. После чего накладываю силу на вектор соединяющий тела и передаю обоим телам с разными знаками.
Сами тела силу делят на свою массу и получают ускорение. Потом вычисляется скорость = старая скорость плюс ускорение умноженное на время. А потом и новое положение - старое положение плюс скорость умноженная на время.

Такими циклами программа считает все положения.
Для Меркурия, Венеры и т.д. эти результаты сходятся с реальностью, а вот для Барицентра не получается.
Что я смог придумать что-бы найти где ошибка:
Первое - длины орбит соотносятся так-же как соотносятся и скорости. То есть а_п/а_х = V_п / V_х Это соотношение сохраняется.
Второе - закон сохранения импульса - m_п*V_п + m_х*V_х = const это тоже выполняется, но как-то странно - на сколько я понимаю импульс Плутона и Харона должен быть одинаковым, а у меня эта сумма хоть и постоянна, но далеко не ноль (причём в зависимости от того, на какой момент рассчитывать начальное положение и скорость и сумма получается разной, хоть и постоянной для каждого раза).
Третье - вектора скорости направлены друг на встречу другу и отличаются длиной.
Собственно это всё, что я смог придумать что-бы понять где ошибка, но так и не понял.

Данные:

Плутон Масса Средняя аномалия Эпоха Орбитальный период, дни Эксцентриситет Большая полуось, км Долгота восходящего угла Аргумент перицентра Наклонение 1.303e+22 93.4395 2000 6.3872304 0.0002 2269.8886541 273.5582 133.7606 0.0048 Харон Масса Средняя аномалия Эпоха Орбитальный период, дни Эксцентриситет Большая полуось, км Долгота восходящего угла Аргумент перицентра Наклонение 1.586e+21 93.4395 2000 6.3872304 0.0002 19635.715 273.5582 313.7606 0.0048 Расчётные данные на 2015 год 1 января полночь Плутон Средняя аномалия Эксцентрическая аномалия Истинная аномалия Вектор положения, км Вектор скорости, м/с 93.4395 93.4509383771 93.4623766854 Х=-1758.59, У=1435.23, Z=-0.139581, Dist=2269.92 X=-16.3444,Y=-20.0188, Z=-0.00147071,Leight=25.8436 Харон Средняя аномалия Эксцентрическая аномалия Истинная аномалия Вектор положения, км Вектор скорости, м/с 93.4395 93.4509383771 93.4623766854 Х=15212.7, У=-12415.5, Z=1.20745, Dist=19636.0 X=141.387,Y=-173.173, Z=-173.173,Leight=223.560 Для подтверждения: ап/ах = Vп / Vх = 2269.92/19636.0 = 25.8436/223.560 = 0.1155999 mп*Vп + mх*Vх = X=1.12731e+22, Y=1.38074e+22, Z=1.01438e+18, Leight=1.78249e+22, я так понял, что это значит что Барицентр движется равномерно и прямолинейно, но разве Барицентр не должен быть неподвижен? И третье - нормализованные вектора скорости: Плутон - X=-0.632434, Y=-0.774614, Z=-5.69080e-05 Харон -   X= 0.632434, Y= 0.774614, Z= 5.69080e-05 Я не могу понять где ошибка. С известными данными это сходится - я имею в виду модули скоростей. Плутон и Харон, надо отметить, всё равно совершают круги друг во круг друга: Если рассматривать Барицентр неподвижным, то они просто улетают в бесконечность, но если Барицентр - это центр масс и его положение каждый раз пересчитывать, то они вращаются вокруг него с периодом около 16 дней, но Плутон имеет перицентр - 2269.92км, а апоцентр около 15000км. Подозреваю, что если увеличить период обращения, то скорости уменьшаться, апоцентр уменьшится и в какой-то момент Плутон и Харон будут летать с эксцентриситетом около 0, но это-же не относится к действительности. И ещё мне не понятно почему импульс Плутона и Харона не равны, но это тоже вроде решиться увеличением периода обращения...

[xd]

То есть Вы интегрируете положение методом правых прямоугольников.
Каждый метод интегрирования имеет погрешности. Тут Вы их наблюдаете во всей красе.

[Крупин]

Не стоит ли отправить в "Астрономию и компьютеры"?

[Tarakan]

То есть Вы интегрируете положение методом правых прямоугольников.
Каждый метод интегрирования имеет погрешности. Тут Вы их наблюдаете во всей красе.

Но чем больше разбиение, тем больше точность и меньше погрешность, так?
Я разбиваю по 5 секунд, а это 0.00906032%, это достаточно маленькая величина, а погрешность уж слишком большая. К тому-же интегрирование по 50 секунд приводит к абсолютно такому-же результату (с большой точностью)
Начальные условия вычисляются с точностью до 17 знаков мантиссы, а сами расчёты до 190 знаков мантиссы, так что в точности расчётов на не большой срок я уверен.

Не стоит ли отправить в "Астрономию и компьютеры"?

Я-же не рассказываю как я считал, хоть на бумаге, я просто говорю что при таких начальных условиях получается странная вещь)

[xd]

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

[Tarakan]

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

Но разве ошибка при 5с и 50с не должна различаться?
К тому-же - разве импульсы Плутона и Харона не должны быть равны?

[xd]

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

[Ajaja]

Можно попробовать взять исходные данные с
http://ssd.jpl.nasa.gov/horizons.cgi

Код: [Выделить]

Target body name: Pluto (999)                     {source: plu055l_merged}
Center body name: Pluto Barycenter (9)            {source: plu055l_merged}
Center-site name: BODY CENTER
*******************************************************************************
Output units    : KM-S                                                         
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0                                                 
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch                 
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
$$SOE
2457023.500000000 = A.D. 2015-Jan-01 00:00:00.0000 TDB
  -1.469439551574074E+03 -7.555931024424377E+02  1.350215119844029E+03
   6.238419402226713E-03  1.693834777870799E-02  1.627254565976924E-02
   7.117709438900199E-03  2.133835608017691E+03  2.792263306690860E-06

Код: [Выделить]

Target body name: Charon (901)                    {source: plu055l_merged}
Center body name: Pluto Barycenter (9)            {source: plu055l_merged}
Center-site name: BODY CENTER
*******************************************************************************
Output units    : KM-S                                                         
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0                                                 
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch                 
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
$$SOE
2457023.500000000 = A.D. 2015-Jan-01 00:00:00.0000 TDB
   1.202424140083403E+04  6.183415842751257E+03 -1.104788638477244E+04
  -5.104149256575340E-02 -1.385895457549670E-01 -1.331441087647647E-01
   5.824231666426191E-02  1.746060727239344E+04  1.540413229915741E-05



т.к. скорей-всего проблема в

Моя программа за основу берёт Кеплеровские элементы орбиты и превращает их в положение тела.

[Astro_Coder]

я так понял, что это значит что Барицентр движется равномерно и прямолинейно, но разве Барицентр не должен быть неподвижен?

Барицентр системы Плутон-Харон движется вокруг Солнца. Из-за того что размер орбиты большой, вот и получилось на небольшом отрезке времени почти равномерное и прямолинейное движение. По идее, нужно брать систему 3-х тел Солнце-Плутон-Харон и ее интегрировать в координатах общего барицентра.

И, как выше сказали, исходные данные брать с horizons'а.

[Mikola]

Для численного расчета положений небесных тел законов Ньютона и Кеплера недостаточно, они, увы, "грубоваты" и слишком быстро накапливают ошибки, что приводит к "полному улету" траекторий от реальных за достаточно короткое время. Для приближенных к реальности расчетов спросите ОКгугл по теме "расчет эфемерид". Еще более глубоко можно запросить модели DE, VSOP, EPM. Не удивляйтесь, что в этих моделях не увидите ничего похожего на законы Ньютона и Кеплера (все таки небесная механика несколько сложнее полета ядра в поле тяготения Земли).

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

[VLG]

«Астрономические формулы для калькуляторов» Меёс Ж. 1988г.
http://эквм.рф/biblioteka/0020/0020.htm
у  меня была такая книга живьем, продал тут. Но для эфемерид там приводились разложения в ряд по времени первой-второй степени. .

[Крупин]

Для численного расчета положений небесных тел законов Ньютона и Кеплера недостаточно, они, увы, "грубоваты" и слишком быстро накапливают ошибки,

Законы Кеплера, конечно, неточны, но Ньютон - дело другое. Релятивистские поправки почти несущественны, тем более для дальних объектов - Плутона и Харона. накопление ошибок никак с Ньютоном не связано, а определяется вычислительным методом и длиной шага интегрирования.

[Mikola]

Для численного расчета положений небесных тел законов Ньютона и Кеплера недостаточно, они, увы, "грубоваты" и слишком быстро накапливают ошибки,

Законы Кеплера, конечно, неточны, но Ньютон - дело другое. Релятивистские поправки почти несущественны, тем более для дальних объектов - Плутона и Харона. накопление ошибок никак с Ньютоном не связано, а определяется вычислительным методом и длиной шага интегрирования.

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

[Ajaja]

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

[xd]

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

В контексте исходной задачи либо усложняете, либо говорите ерунду.

Основной источник погрешностей, которые надо устранить - именно погрешности численного метода.

[Ajaja]

В контексте исходной задачи либо усложняете, либо говорите ерунду.

Основной источник погрешностей, которые надо устранить - именно погрешности численного метода.

В контексте исходной задачи (система Плутона-Харона) основная проблема не в численных методах (там точность можно настраивать), а в недостатке исходных данных. Хотя Новые Горизонты ситуацию немного и исправил.

[xd]

В контексте задачи двух тел (система Плутон-Харон взята лишь в качестве примера, можно спокойно заменить на Тело А и Тело Б) проблема только в численном методе.

[Ajaja]

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


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

[xd]

Предполагаем начальные условия известными с достаточной точностью. Пускай будет идеализированная система, во вселенной, где нет иных тел.  Моделируем движение, и система начинает "уползать". Это указывает на косяки метода, независимо от данных. Конечно, это не отменяет неточность исходных данных и всякие физические эффекты второго порядка малости, но в данном случае всё указывает именно на очень грубый метод.

[Mikola]

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

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