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

[sharp]

проблема только в численном методе

Какой посоветуете?

[xd]

Предлагаю начать с прочтения вот этого: http://www.sai.msu.ru/neb/pcm/pcm04_11.pdf и http://crecs.ru/ru/numlabs/help/Integr/Integr.html#_Toc190585558
После этого можно будет по существу поговорить.

[Ajaja]

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

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

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

[Tarakan]

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

Можно, но если эти данные использовать, то получается у меня масса Плутона и Харона не достаточно точно известны, а это уже другой вопрос)

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

Я рассматриваю их скорости относительно Барицентра. Я уже взял относительные скорости, так как если-бы я взял абсолютные, то их нормализованные версии ни как не были-бы равны по углу и различны по направлению. А относительно своего Барицентра, как мне кажется, они не должны убегать.
Но я могу ошибаться в выводах и поэтому промоделировал процесс подлинее, на промежутке в год это число не изменилось, а уж за 1/248 орбита не хило так изменяется

Там были неплохие упрощенные модели расчета эфемерид, которые не "уплывали" на промежутках в несколько десятков лет

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

Моделируем движение, и система начинает "уползать".

Я не против уползания, но одно дело - когда эффект проявляется спустя месяцы, годы, а другое секунды и дни

Предлагаю начать с прочтения вот этого

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

Но все в итоге упертся в погрешность исходных данных.

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

[Ajaja]

Можно, но если эти данные использовать, то получается у меня масса Плутона и Харона не достаточно точно известны, а это уже другой вопрос)

https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/plu055.cmt
GM - гравитационные параметры тел

[Tarakan]

https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/plu055.cmt
GM - гравитационные параметры тел

Спасибо, теперь у меня с достаточно большой точностью сошлось! Скорость убегания центра масс реального от расчётного миллиметры в месяц!

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

[Ajaja]

https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/plu055.cmt
GM - гравитационные параметры тел

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

Ну так сравните с тем же Horizons, он умеет и элементы орбиты генерировать.

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

Ephemeris Type [change] : ELEMENTS
Target Body [change] : Charon (PI) [901]
Center [change] : Pluto (body center) [500@999]
Time Span [change] : Start=2015-01-01, Stop=2015-01-10, Step=1 d
Table Settings [change] : output units=KM-S
...
Keplerian GM    : 9.7559000499039507E+02 km^3/s^2
...
2457023.500000000 = A.D. 2015-Jan-01 00:00:00.0000 TDB
 EC= 1.477380420046948E-04 QR= 1.959396192601740E+04 IN= 1.128963726104548E+02
 OM= 2.274017192362429E+02 W = 1.898761621437299E+02 Tp=  2457022.905704734381
 N = 6.523439875612244E-04 MA= 3.349597911873988E+01 TA= 3.350532360219091E+01
 A = 1.959685712731884E+04 AD= 1.959975232862027E+04 PR= 5.518560864580865E+05

   EC     Eccentricity, e                                                   
      QR     Periapsis distance, q (km)                                       
      IN     Inclination w.r.t XY-plane, i (degrees)                           
      OM     Longitude of Ascending Node, OMEGA, (degrees)                     
      W      Argument of Perifocus, w (degrees)                               
      Tp     Time of periapsis (Julian Day Number)                             
      N      Mean motion, n (degrees/sec)                                     
      MA     Mean anomaly, M (degrees)                                         
      TA     True anomaly, nu (degrees)                                       
      A      Semi-major axis, a (km)                                           
      AD     Apoapsis distance (km)                                           
      PR     Sidereal orbit period (sec)         


[Tarakan]

Ну так сравните с тем же Horizons, он умеет и элементы орбиты генерировать.

Я знаю какие числа должны получится, но у меня даже самое элементарное - большая полуось - получилась в два раза меньше. Единственное что совпало - наклонение орбиты. А так как все элементы рассчитываются только из двух констант, то логично предположить что я их неправильно посчитал. (У меня наибольшие подозрения на счёт постоянной энергии, но я не знаю чему она должны быть равна)

[Ajaja]

Ну так сравните с тем же Horizons, он умеет и элементы орбиты генерировать.

Я знаю какие числа должны получится, но у меня даже самое элементарное - большая полуось - получилась в два раза меньше. Единственное что совпало - наклонение орбиты. А так как все элементы рассчитываются только из двух констант, то логично предположить что я их неправильно посчитал. (У меня наибольшие подозрения на счёт постоянной энергии, но я не знаю чему она должны быть равна)

Вы же считаете кеплеровские элементы от тела, а не от барицентра? С константой энергии там все просто, как и с векторной константой площади , но нужно считать все в координатах одного тела относительно другого. Т.е. берете положение Харона относительно Плутона, берете его скорость относительно Плутона, умножаете эти два вектора векторным произведением, получаете векторную константу. Берете квадрат скорости Харона относительно Плутона, вычитаете из него гравтационный  параметр системы  ( Keplerian GM    : 9.7559000499039507E+02 km^3/s^2)  умноженный на два и деленный на расстояние между Плутоном и Хароном - получаете константу энергии.

[Tarakan]

берете его скорость относительно Плутона

А разве скорость относительно Плутона - это не скорость относительно Барицентра? Ведь для скорости не имеет значения где начало координат.

Для Плутона

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

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

Откуда
\[r_{P} = (-1.469439551574074E+03; -7.555931024424377E+02; 1.350215119844029E+03)\]
\[r_{H} = (1.202424140083403E+04; 6.183415842751257E+03; -1.104788638477244E+04)\]
\[V_{P} = (6.238419402226713E-03; 1.693834777870799E-02;  1.627254565976924E-02)\]
\[V_{H} = (-5.104149256575340E-02; -1.385895457549670E-01; -1.331441087647647E-01)\]
\[r_{P-H} = r_P - r_H = (X = -13493.680952408104;Y = -6939.0089451936947;Z = 12398.101504616469)\]
Постоянная энергия h:
\[h = V_P^2 - 2*GM/r_{P-H} = 0.0243027^2 - 2*9.7559000499039507E+02/19594.4 = -0.09898783026\]
А большая полуось соответственно
\[a = - GM / h = 9.7559000499039507E+02 / 0.09898783026 = 9855.65601779\]
А в реальность она только 2133.84.
Если-бы всё было так просто, я-бы и сам справился, но у меня слишком сильно не сходятся результаты.
Правда я не догадался раньше брать точку отсчёта относительно другого тела, а брал относительно Баррицентра и по поводу гравитационного параметра у меня вопрос есть - он-же высчитывается как G*M, У G порядок е-11, а масса Плутона е+22, это без учёта других масс, значит GM должна получится в районе порядка е+11, где-то так.

[Ajaja]

Vp_h =  (0.572799119679801156e-1, 0.155527893533674966, 0.149416654424533946)
h = Vp_h^2 - 2*GM/Rp_h = -0.04978298298
a = 9.7559000499039507e02 / 0.04978298298 = 19596.85713

Сравниваем:

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

Ephemeris / WWW_USER Sat Mar 11 14:36:10 2017 Pasadena, USA      / Horizons   
*******************************************************************************
Target body name: Charon (901)                    {source: plu055l_merged}
Center body name: Pluto (999)                     {source: plu055l_merged}
Center-site name: BODY CENTER
*******************************************************************************
Start time      : A.D. 2015-Jan-01 00:00:00.0000 TDB
Stop  time      : A.D. 2015-Jan-02 00:00:00.0000 TDB
Step-size       : 1440 minutes
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii    : 1195.0 x 1195.0 x 1195.0 km     {Equator, meridian, pole}   
Keplerian GM    : 9.7559000499039507E+02 km^3/s^2
Output units    : KM-S, deg, Julian Day Number (Tp)                           
Output type     : GEOMETRIC osculating elements
Output format   : 10
Reference frame : ICRF/J2000.0                                                 
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch                 
*******************************************************************************
JDTDB
   EC    QR   IN
   OM    W    Tp
   N     MA   TA
   A     AD   PR
*******************************************************************************
$$SOE
2457023.500000000 = A.D. 2015-Jan-01 00:00:00.0000 TDB
 EC= 1.477380420046948E-04 QR= 1.959396192601740E+04 IN= 1.128963726104548E+02
 OM= 2.274017192362429E+02 W = 1.898761621437299E+02 Tp=  2457022.905704734381
 N = 6.523439875612244E-04 MA= 3.349597911873988E+01 TA= 3.350532360219091E+01
 A = 1.959685712731884E+04 AD= 1.959975232862027E+04 PR= 5.518560864580865E+05
A = 1.959685712731884E+04

[Ajaja]

У G порядок е-11

В HORIZONS все в километрах, а не в метрах, а  в GM  м^3  и км^3 , соответственно разница в числах на 9 порядков.

[Tarakan]

Vp_h =  (0.572799119679801156e-1, 0.155527893533674966, 0.149416654424533946)

Что-то я не пойму как вы это получили. Взять хотя-бы часть Z - у Харона и Плутона там е-1 и е-2, а у вас получается е+0, даже если складывать не получается.

[Ajaja]

Vp_h =  (0.572799119679801156e-1, 0.155527893533674966, 0.149416654424533946)

Что-то я не пойму как вы это получили. Взять хотя-бы часть Z - у Харона и Плутона там е-1 и е-2, а у вас получается е+0, даже если складывать не получается.

1.627254565976924E-02  - (-1.331441087647647E-01)  = 0.14941665442453394

[Tarakan]

Сравниваем:

А теперь для Плутона:
Vh-p = (X = -0.057279911967980113;Y = -0.15552789353367499;Z = -0.14941665442453394;l = 0.223148)
h =Vh_p^2 - 2*GM/Rh_p = -0.04978342159, а соответственно
a = 9.7559000499039507e02 / 0.04978342159= 19596.6844751.
И это логично, ведь скорость Плутона относительно Харона и Харона относительно Плутона - одно и то-же, а следовательно и остальное совпадет. С Хароном у меня всё сходится, а Плутон теперь не попадает(

[Ajaja]

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

[Tarakan]

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

Можете пояснить, а то я не понимаю откуда это число надо брать.

[Ajaja]

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

Можете пояснить, а то я не понимаю откуда это число надо брать.

Если не ошибаюсь, находим на большой полуоси барицентр, расстояние до него - соответственно большие полуоси орбит вокруг барицентра, из них и периода орбиты (вычисленного раньше) получаем эти GM (по формуле, связывaющей период, полуось и GM).
Но я особого смысла в таких вычислениях не вижу. Достаточно одной нормальной кеплеровской орбиты, например Харона вокруг Плутона,  а барицентр и положения тел относительно него можно уже потом вычислять.

[Tarakan]

Но я особого смысла в таких вычислениях не вижу.

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

(по формуле, связывaющей период, полуось и GM)

Я, к сожалению, ещё не сталкивался с такой формулой. Расскажите пожалуйста поподробнее про этот момент и про расчёт в таком случае.

[Дрюша]

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

Ну, Рунге-Кутт четвёртого порядка, я полагаю, тут его хватит.

Но одно непонятно. Вот, у Плутона с Хароном - какие-то проблемы. А у Земли с Луной - что, никаких? По идее, Земля с Луной должны быть круче.