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$	[-3.5023653, -3.8169847, -1.5507963]	[0.00565429, -0.00412490, -0.00190589]
Saturn	$m_2 = 0.000285583733151$	[9.0755314, -3.0458353, -1.6483708]	[0.00168318, 0.00483525, 0.00192462]
Uranus	$m_3 = 0.0000437273164546$	[8.3101420, -16.2901086, -7.2521278]	[0.00354178, 0.00137102, 0.00055029]
Neptune	$m_4 = 0.0000517759138449$	[11.4707666, -25.7294829, -10.8169456]	[0.00288930, 0.00114527, 0.00039677]
Pluto	$m_5 = 1/(1.3 \cdot 10^8 )$	[-15.5387357, -25.2225594, -3.1902382]	[0.00276725, -0.00170702, -0.00136504]

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}\]

where each $p_i$ and $q_i$ is a 3-dimensional vector describing the planet's position and momentum, respectively.

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 simplify 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")