Simulating the Outer Solar System

Data

The chosen units are masses relative to the sun, meaning the sun has mass $1$. We have taken $m_0 = 1.00000597682$ to take account of the inner planets. Distances are in astronomical units, times in earth days, and the gravitational constant is thus $G = 2.95912208286 \cdot 10^{-4}$.

planet	mass	initial position	initial velocity
Jupiter	$m_1 = 0.000954786104043$	<ul><li>-3.5023653</li><li>-3.8169847</li><li>-1.5507963</li></ul>	<ul><li>0.00565429</li><li>-0.00412490</li><li>-0.00190589</li></ul>
Saturn	$m_2 = 0.000285583733151$	<ul><li>9.0755314</li><li>-3.0458353</li><li>-1.6483708</li></ul>	<ul><li>0.00168318</li><li>0.00483525</li><li>0.00192462</li></ul>
Uranus	$m_3 = 0.0000437273164546$	<ul><li>8.3101420</li><li>-16.2901086</li><li>-7.2521278</li></ul>	<ul><li>0.00354178</li><li>0.00137102</li><li>0.00055029</li></ul>
Neptune	$m_4 = 0.0000517759138449$	<ul><li>11.4707666</li><li>-25.7294829</li><li>-10.8169456</li></ul>	<ul><li>0.00288930</li><li>0.00114527</li><li>0.00039677</li></ul>
Pluto	$ m_5 = 1/(1.3 \cdot 10^8 )$	<ul><li>-15.5387357</li><li>-25.2225594</li><li>-3.1902382</li></ul>	<ul><li>0.00276725</li><li>-0.00170702</li><li>-0.00136504</li></ul>

The data is taken from the book “Geometric Numerical Integration” by E. Hairer, C. Lubich and G. Wanner.

import Plots, OrdinaryDiffEq as ODE
import ModelingToolkit as MTK
using ModelingToolkit: t_nounits as t, D_nounits as D, @mtkbuild, @variables
Plots.gr()

G = 2.95912208286e-4
M = [
    1.00000597682,
    0.000954786104043,
    0.000285583733151,
    0.0000437273164546,
    0.0000517759138449,
    1 / 1.3e8
]
planets = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"]

pos = [0.0 -3.5023653 9.0755314 8.310142 11.4707666 -15.5387357
       0.0 -3.8169847 -3.0458353 -16.2901086 -25.7294829 -25.2225594
       0.0 -1.5507963 -1.6483708 -7.2521278 -10.8169456 -3.1902382]
vel = [0.0 0.00565429 0.00168318 0.00354178 0.0028893 0.00276725
       0.0 -0.0041249 0.00483525 0.00137102 0.00114527 -0.00170702
       0.0 -0.00190589 0.00192462 0.00055029 0.00039677 -0.00136504]
tspan = (0.0, 200_000.0)

(0.0, 200000.0)

The N-body problem's Hamiltonian is

\[H(p,q) = \frac{1}{2}\sum_{i=0}^{N}\frac{p_{i}^{T}p_{i}}{m_{i}} - G\sum_{i=1}^{N}\sum_{j=0}^{i-1}\frac{m_{i}m_{j}}{\left\lVert q_{i}-q_{j} \right\rVert}\]

Here, we want to solve for the motion of the five outer planets relative to the sun, namely, Jupiter, Saturn, Uranus, Neptune, and Pluto.

const ∑ = sum
const N = 6
@variables u(t)[1:3, 1:N]
u = collect(u)
potential = -G *
            ∑(
    i -> ∑(j -> (M[i] * M[j]) / √(∑(k -> (u[k, i] - u[k, j])^2, 1:3)), 1:(i - 1)),
    2:N)

\[ \begin{equation} - 0.00029591 \left( \frac{3.3636 \cdot 10^{-13}}{\sqrt{\left( - u\_{1,4}\left( t \right) + u\_{1,6}\left( t \right) \right)^{2} + \left( - u\_{2,4}\left( t \right) + u\_{2,6}\left( t \right) \right)^{2} + \left( - u\_{3,4}\left( t \right) + u\_{3,6}\left( t \right) \right)^{2}}} + \frac{1.2488 \cdot 10^{-8}}{\sqrt{\left( - u\_{1,3}\left( t \right) + u\_{1,4}\left( t \right) \right)^{2} + \left( - u\_{2,3}\left( t \right) + u\_{2,4}\left( t \right) \right)^{2} + \left( - u\_{3,3}\left( t \right) + u\_{3,4}\left( t \right) \right)^{2}}} + \frac{2.1968 \cdot 10^{-12}}{\sqrt{\left( - u\_{1,3}\left( t \right) + u\_{1,6}\left( t \right) \right)^{2} + \left( - u\_{2,3}\left( t \right) + u\_{2,6}\left( t \right) \right)^{2} + \left( - u\_{3,3}\left( t \right) + u\_{3,6}\left( t \right) \right)^{2}}} + \frac{3.9828 \cdot 10^{-13}}{\sqrt{\left( - u\_{1,5}\left( t \right) + u\_{1,6}\left( t \right) \right)^{2} + \left( - u\_{2,5}\left( t \right) + u\_{2,6}\left( t \right) \right)^{2} + \left( - u\_{3,5}\left( t \right) + u\_{3,6}\left( t \right) \right)^{2}}} + \frac{4.3728 \cdot 10^{-5}}{\sqrt{\left( - u\_{1,1}\left( t \right) + u\_{1,4}\left( t \right) \right)^{2} + \left( - u\_{2,1}\left( t \right) + u\_{2,4}\left( t \right) \right)^{2} + \left( - u\_{3,1}\left( t \right) + u\_{3,4}\left( t \right) \right)^{2}}} + \frac{4.175 \cdot 10^{-8}}{\sqrt{\left( - u\_{1,2}\left( t \right) + u\_{1,4}\left( t \right) \right)^{2} + \left( - u\_{2,2}\left( t \right) + u\_{2,4}\left( t \right) \right)^{2} + \left( - u\_{3,2}\left( t \right) + u\_{3,4}\left( t \right) \right)^{2}}} + \frac{5.1776 \cdot 10^{-5}}{\sqrt{\left( - u\_{1,1}\left( t \right) + u\_{1,5}\left( t \right) \right)^{2} + \left( - u\_{2,1}\left( t \right) + u\_{2,5}\left( t \right) \right)^{2} + \left( - u\_{3,1}\left( t \right) + u\_{3,5}\left( t \right) \right)^{2}}} + \frac{2.7267 \cdot 10^{-7}}{\sqrt{\left( - u\_{1,2}\left( t \right) + u\_{1,3}\left( t \right) \right)^{2} + \left( - u\_{2,2}\left( t \right) + u\_{2,3}\left( t \right) \right)^{2} + \left( - u\_{3,2}\left( t \right) + u\_{3,3}\left( t \right) \right)^{2}}} + \frac{4.9435 \cdot 10^{-8}}{\sqrt{\left( - u\_{1,2}\left( t \right) + u\_{1,5}\left( t \right) \right)^{2} + \left( - u\_{2,2}\left( t \right) + u\_{2,5}\left( t \right) \right)^{2} + \left( - u\_{3,2}\left( t \right) + u\_{3,5}\left( t \right) \right)^{2}}} + \frac{0.00095479}{\sqrt{\left( - u\_{1,1}\left( t \right) + u\_{1,2}\left( t \right) \right)^{2} + \left( - u\_{2,1}\left( t \right) + u\_{2,2}\left( t \right) \right)^{2} + \left( - u\_{3,1}\left( t \right) + u\_{3,2}\left( t \right) \right)^{2}}} + \frac{7.6924 \cdot 10^{-9}}{\sqrt{\left( - u\_{1,1}\left( t \right) + u\_{1,6}\left( t \right) \right)^{2} + \left( - u\_{2,1}\left( t \right) + u\_{2,6}\left( t \right) \right)^{2} + \left( - u\_{3,1}\left( t \right) + u\_{3,6}\left( t \right) \right)^{2}}} + \frac{0.00028559}{\sqrt{\left( - u\_{1,1}\left( t \right) + u\_{1,3}\left( t \right) \right)^{2} + \left( - u\_{2,1}\left( t \right) + u\_{2,3}\left( t \right) \right)^{2} + \left( - u\_{3,1}\left( t \right) + u\_{3,3}\left( t \right) \right)^{2}}} + \frac{1.4786 \cdot 10^{-8}}{\sqrt{\left( - u\_{1,3}\left( t \right) + u\_{1,5}\left( t \right) \right)^{2} + \left( - u\_{2,3}\left( t \right) + u\_{2,5}\left( t \right) \right)^{2} + \left( - u\_{3,3}\left( t \right) + u\_{3,5}\left( t \right) \right)^{2}}} + \frac{2.264 \cdot 10^{-9}}{\sqrt{\left( - u\_{1,4}\left( t \right) + u\_{1,5}\left( t \right) \right)^{2} + \left( - u\_{2,4}\left( t \right) + u\_{2,5}\left( t \right) \right)^{2} + \left( - u\_{3,4}\left( t \right) + u\_{3,5}\left( t \right) \right)^{2}}} + \frac{7.3445 \cdot 10^{-12}}{\sqrt{\left( - u\_{1,2}\left( t \right) + u\_{1,6}\left( t \right) \right)^{2} + \left( - u\_{2,2}\left( t \right) + u\_{2,6}\left( t \right) \right)^{2} + \left( - u\_{3,2}\left( t \right) + u\_{3,6}\left( t \right) \right)^{2}}} \right) \end{equation} \]

Hamiltonian System

NBodyProblem constructs a second order ODE problem under the hood. We know that a Hamiltonian system has the form of

\[\dot{p} = -H_{q}(p,q)\quad \dot{q}=H_{p}(p,q)\]

For an N-body system, we can symplify this as:

\[\dot{p} = -\nabla{V}(q)\quad \dot{q}=M^{-1}p.\]

Thus, $\dot{q}$ is defined by the masses. We only need to define $\dot{p}$, and this is done internally by taking the gradient of $V$. Therefore, we only need to pass the potential function and the rest is taken care of.

eqs = vec(@. D(D(u))) .~ .-MTK.gradient(potential, vec(u)) ./
                         repeat(M, inner = 3)
@mtkbuild sys = MTK.System(eqs, t)
prob = ODE.ODEProblem(sys, [vec(u .=> pos); vec(D.(u) .=> vel)], tspan)
sol = ODE.solve(prob, ODE.Tsit5());

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 862-element Vector{Float64}:
      0.0
      0.43456238545920783
      4.780186240051286
     48.23642478597206
    195.27337773985937
    404.37880043365203
    641.2461951827266
    917.4024243418721
   1259.375074434556
   1620.723519074113
      ⋮
 199469.8984301734
 199550.292229572
 199620.51170967382
 199697.87840262795
 199765.84076245074
 199832.36143102907
 199887.04316168063
 199959.3068873699
 200000.0
u: 862-element Vector{Vector{Float64}}:
 [-3.1902382, -0.00136504, -25.2225594, -0.00170702, -15.5387357, 0.00276725, -10.8169456, 0.00039677, -25.7294829, 0.00114527  …  -3.8169847, -0.0041249, -3.5023653, 0.00565429, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 [-3.190831391666314, -0.0013650244789412218, -25.223301179998256, -0.001706897186454997, -15.537533140796315, 0.002767325672882192, -10.816773167663484, 0.0003968207122091521, -25.728985182525822, 0.0011453906312636451  …  -3.818776551780467, -0.004121794983806699, -3.4999075394383845, 0.005657137663985553, -2.77776298862327e-10, -1.2785303085738865e-9, -6.745773496345314e-10, -3.104871144001842e-9, -5.098303779711312e-10, -2.346105131364384e-9]
 [-3.196762937128868, -0.0013648691152799539, -25.230716044387762, -0.0017056688965533717, -15.52550574063629, 0.002768082053772159, -10.815047632292558, 0.0003973277927111847, -25.724005124769256, 0.0011465968224111933  …  -3.8366207222927584, -0.004090656421676488, -3.4752620277231188, 0.005685511868486027, -3.36696487911887e-8, -1.4100666088910384e-8, -8.175135524038143e-8, -3.4233554324164556e-8, -6.152988189104976e-8, -2.5706941518695228e-8]
 [-3.2560410257920296, -0.0013633001902609797, -25.304570881838067, -0.0016933706555559782, -15.405051486626984, 0.0027756110286434707, -10.797671145111675, 0.0004023944185903396, -25.673916441625973, 0.0011586465368538337  …  -4.007487516495138, -0.003770599600091882, -3.2221835308829867, 0.0059587932915218115, -3.4874659022100094e-6, -1.4594000158084996e-7, -8.45230700457103e-6, -3.53353735173139e-7, -6.100402511077371e-6, -2.4909025937027336e-7]
 [-3.4560950438194866, -0.001357786536892723, -25.55048864881513, -0.001651557531681458, -14.995086973641925, 0.002800613918439817, -10.737247074610087, 0.00041948022029475646, -25.500564258550423, 0.0011992504322559396  …  -4.4763282740831025, -0.0025806560822637585, -2.2862427084241728, 0.006730897201673277, -6.024100672070627e-5, -6.369559649765036e-7, -0.0001451482502017574, -1.529138371536945e-6, -9.029656259744727e-5, -8.562666292142268e-7]
 [-3.7391503253370364, -0.0013494060257915718, -25.88957839011003, -0.001591589478503688, -14.405852122492028, 0.002834905191574055, -10.647004315296812, 0.00044361966174052704, -25.24379557801365, 0.001256535830385216  …  -4.819472301720937, -0.0006664593990847885, -0.8022848285526655, 0.007365006227690945, -0.0002747692410211364, -1.430096198203768e-6, -0.0006569678433674034, -3.396085915620486e-6, -0.0003224635321908885, -1.2700600202222442e-6]
 [-4.057582856734454, -0.0013391609150743518, -26.25846463383041, -0.001523008173343417, -13.729932034469996, 0.002871931193437516, -10.53870984206983, 0.0004707259464700574, -24.938543951977763, 0.0013207457633296043  …  -4.703474908001999, 0.0016548599609000731, 0.9542708808238635, 0.007324487906032102, -0.0007276022576520447, -2.398512185960335e-6, -0.0017256604570546608, -5.631384329495253e-6, -0.0006083763727699927, -1.0075133892324712e-6]
 [-4.425639265202051, -0.0013262311425467392, -26.667923865275117, -0.0014422802284886843, -12.931149248426918, 0.0029126298323277314, -10.404390682372131, 0.0005019874376483292, -24.563583836576818, 0.001394641388147916  …  -3.882175023379029, 0.004236424139516364, 2.845652151911805, 0.006174492073878071, -0.0015432638606036754, -3.487397690399148e-6, -0.00362761106169569, -8.082240979433938e-6, -0.0007246641716727011, 3.554671854327595e-7]
 [-4.876236753538832, -0.001308794260336332, -27.143912620248173, -0.001341340187806172, -11.927010904521099, 0.002959304010544335, -10.226181232787367, 0.0005401516040701686, -24.071220358929434, 0.0014846157954746797  …  -1.9952056202290902, 0.006575790764247332, 4.495930608252941, 0.003230725705352898, -0.002923904354474405, -4.497046006563047e-6, -0.006799227714296791, -1.0242905047747271e-5, -0.00010701718486927201, 3.4901821399911367e-6]
 [-5.345590814017202, -0.0012887220283255018, -27.60919076445552, -0.0012337563665169732, -10.849435759690564, 0.0030041096081688755, -10.023818065895448, 0.0005797663178001924, -23.517872437082897, 0.0015777192883279383  …  0.5606535208229736, 0.007201496082832887, 4.92017444073576, -0.0009678270857685428, -0.004632370364695403, -4.805134987913165e-6, -0.010648108542765499, -1.069645302652259e-5, 0.0019228578780514536, 7.821335996287854e-6]
 ⋮
 [-15.123524109751513, 0.00013284020319727106, -14.882998244513848, 0.0020747246888518423, 35.67769176029667, 0.0017136284822634343, 7.436620121868759, 0.0008806864102202099, 19.539224015936842, 0.002013531478450076  …  -1.1643976782207888, 0.01134288637029172, 2.3730907200049622, 0.008108454173305517, -0.2473389548602936, -5.093603697278523e-6, -0.49304821487901596, -1.185699607645447e-5, 1.2328033605875286, -1.3756693948065007e-6]
 [-15.112625423639356, 0.00013828798207319754, -14.715992016328638, 0.0020799658018381775, 35.81494863277498, 0.0017009850441485415, 7.507142974535007, 0.0008737381154511724, 19.700372699422818, 0.001995427209127771  …  -0.12493043015207193, 0.012846264984274715, 2.541228724228247, -0.004187633084377313, -0.24781180439440806, -5.99276530986655e-6, -0.49412351319702025, -1.329959030311957e-5, 1.2331516669364948, 1.0290996628378757e-5]
 [-15.102748478644132, 0.00014302510354387141, -14.569779165569413, 0.0020844764593477808, 35.934003211244324, 0.0016899424648723028, 7.568281794283585, 0.0008676173706800154, 19.83993185715745, 0.0019794958633450083  …  0.5915803013161267, 0.006546782898839691, 1.9160297314474934, -0.012770349595343568, -0.24816439512662802, -3.6152094928204766e-6, -0.49488054215873767, -7.293774303388143e-6, 1.2341882364715295, 1.8422204096992478e-5]
 [-15.091481930070932, 0.0001482215247894014, -14.408320171689876, 0.0020893741126009864, 36.06427786307255, 0.0016777770342673783, 7.635143840047986, 0.0008608182155292477, 19.99239586237846, 0.0019618165802084606  …  0.675032566286632, -0.004589223538620178, 0.813224447350581, -0.014001165263383377, -0.2482737623072602, 9.096653957515305e-7, -0.4950413772692069, 3.3263048234364494e-6, 1.2357204412008833, 1.9528391601416645e-5]
 [-15.081253895317323, 0.0001527665612371363, -14.266176958634356, 0.002093614408274302, 36.177940438353644, 0.0016670915913214366, 7.693442752448569, 0.0008547979350096965, 20.125194697411878, 0.001946177878148452  …  0.07533781104715866, -0.01224547721394641, 0.06405171436535234, -0.006942679994640394, -0.24808961687670428, 4.202101508766178e-6, -0.49454096133041286, 1.06232374929685e-5, 1.2368523742826898, 1.2729356638888073e-5]
 [-15.070944294976355, 0.00015719732180050376, -14.126771839281062, 0.0020977088789344085, 36.28848866089309, 0.0016566341977526057, 7.750107303014459, 0.0008488625269245839, 20.254143926393027, 0.0019307734270636034  …  -0.817011865694014, -0.013316857464432426, -0.04371633207359502, 0.0038863394000273107, -0.24777020061184934, 4.887272127162588e-6, -0.4937605230399095, 1.1631352784062142e-5, 1.2373591654262124, 2.332407338051759e-6]
 [-15.062249199909228, 0.00016082628675612338, -14.011974422384405, 0.002101033421118864, 36.37884128435904, 0.0016480391398844246, 7.796390444864761, 0.0008439519562307502, 20.359374079521185, 0.0019180389122714194  …  -1.4404523397585942, -0.00875982851352661, 0.39138425411842737, 0.011515694611077484, -0.24754075210021056, 3.1960236113747064e-6, -0.493224878656107, 7.266774386690968e-6, 1.237272902913443, -4.998629839283249e-6]
 [-15.050454548041067, 0.00016560378353686677, -13.859988838214676, 0.002105370115605544, 36.49752437667235, 0.0016366823258962035, 7.857141812876675, 0.0008374192140909384, 20.497367771613536, 0.0019011115661871172  …  -1.7118857496101654, 0.0016317426524223222, 1.3986922296917232, 0.014801736432386369, -0.2474538537389069, -9.844336242577215e-7, -0.49304569527143566, -2.674826944228626e-6, 1.236742262800202, -8.196748588569976e-6]
 [-15.043661040110429, 0.00016828493415562527, -13.77426559455976, 0.00210778385883905, 36.563995971518366, 0.0016302880834588448, 7.891143772357861, 0.0008337189626961726, 20.574535097530337, 0.0018915305767110075  …  -1.5256130867342341, 0.007352366100350984, 1.9645916784360877, 0.012543989239067051, -0.24754386361021266, -3.380358119435055e-6, -0.4932692391451687, -8.148941151253767e-6, 1.2364428080041467, -6.0747229013532665e-6]

plt = Plots.plot()
for i in 1:N
    Plots.plot!(plt, sol, idxs = (u[:, i]...,), lab = planets[i])
end
Plots.plot!(plt; xlab = "x", ylab = "y", zlab = "z", title = "Outer solar system")