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 × 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 × 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 OrdinaryDiffEq as ODE
using ModelingToolkit: System, t_nounits as t, D_nounits as D, @mtkbuild, @variables
using Symbolics: gradient
using Plots: plot, plot!
G = 2.95912208286e-4
M = [
1.00000597682,
9.54786104043e-4,
2.85583733151e-4,
4.37273164546e-5,
5.17759138449e-5,
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} ∑_{i=0}^N \frac{p_i^T p_i}{m_i} - G ∑_{i=1}^N ∑_{j=0}^{i-1} \frac{m_i m_j}{\left\| q_i - q_j \right\|}\]
where each $q_i$ and $p_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{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{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{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}}} + \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.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{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{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{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{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{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.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{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{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{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}}} \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} = -\frac{∂H}{∂q}, \quad \dot{q} = \frac{∂H}{∂p}\]
For an N-body system, we can simplify this as:
\[\dot{p} = -∇ 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))) .~ .-gradient(potential, vec(u)) ./ repeat(M, inner = 3)
@mtkbuild sys = 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.4345623854592078
37.90065823152096
178.68563700235498
385.53178881782515
618.3405413763061
889.0455659758934
1226.3237947453283
1586.7960088544032
1961.2423186888175
⋮
199470.04374603156
199541.0448501004
199616.5773636085
199686.17196112013
199761.18431651435
199827.5708158215
199900.45121765035
199969.80962167826
200000.0
u: 862-element Vector{Vector{Float64}}:
[-3.1902382, -10.8169456, -7.2521278, -1.6483708, -1.5507963, 0.0, -25.2225594, -25.7294829, -16.2901086, -3.0458353 … 0.00137102, 0.00483525, -0.0041249, 0.0, 0.00276725, 0.0028893, 0.00354178, 0.00168318, 0.00565429, 0.0]
[-3.190831391666314, -10.816773167663484, -7.251888638015288, -1.6475343823037847, -1.5516242540172815, -2.777762988623252e-10, -25.223301179998256, -25.728985182525822, -16.289512746418094, -3.0437339894169093 … 0.0013712954942851254, 0.004835677090231295, -0.004121794983806699, -3.1048711440018414e-9, 0.002767325672882192, 0.0028892462117755886, 0.0035416394156575887, 0.0016819059109041002, 0.005657137663985553, -2.3461051313643855e-9]
[-3.241948330485552, -10.801823976210718, -7.231068987293383, -1.575049498744649, -1.6209113188793949, -2.1444623584402983e-6, -25.287053467383586, -25.685877142405612, -16.2376913235701, -2.8618835103113174 … 0.0013949986875277222, 0.004871442971061505, -0.003848129945124167, -2.7618358170406887e-7, 0.00277382605956878, 0.002884586502450322, 0.0035294264124931, 0.0015715064759500212, 0.00589554669961706, -1.9767103667688806e-7]
[-3.4335671375376173, -10.744189356298108, -7.149320167071274, -1.2965336982997142, -1.8414928954402379, -5.0162228641070864e-5, -25.52305378686038, -25.520419235571733, -16.035071941551116, -2.1675813315616264 … 0.0014831919402959173, 0.004986773626543438, -0.0027224352250950416, -1.3897820189116103e-6, 0.002797829859691834, 0.0028666833025901904, 0.0034819138583236254, 0.0011474857815938649, 0.00665635178601216, -7.999757819769807e-7]
[-3.7137106858976683, -10.65534479409368, -7.019274196072119, -0.8725710519906974, -2.036809451249237, -0.00024852295196801786, -25.85953052636742, -25.26742908845482, -15.715092100302337, -1.1232573047741323 … 0.0016101753659002887, 0.00509923134392997, -0.0008467895445745209, -3.2212170367279017e-6, 0.002831875867247711, 0.0028392573253226325, 0.003407547772275999, 0.0005021095088762221, 0.007332344335878154, -1.2564208091642303e-6]
[-4.026896745518497, -10.549462243214293, -6.859062677097583, -0.38084579491990617, -2.0504190378786276, -0.0006737431640637975, -26.22350284768089, -24.96872577084524, -15.323961230110575, 0.07045730112338254 … 0.0017492293699314775, 0.005140141607894802, 0.001429549535311931, -5.415656009329909e-6, 0.0028684356953043036, 0.0028068053263882895, 0.003317543566139734, -0.0002469855375663619, 0.007366042731032005, -1.0690517355334715e-6]
[-4.388011985023564, -10.418580215386676, -6.654860202702344, 0.19904224096715833, -1.7790928353632138, -0.0014458506901734443, -26.626907285009917, -24.603024630917634, -14.8291509828353, 1.4550859101647695 … 0.0019054864985540553, 0.0050677075335554935, 0.003989956549850301, -7.85018532822309e-6, 0.002908574032740688, 0.0027669859919848836, 0.0032048119437642145, -0.0011316130682388927, 0.006347103551167493, 1.6343592556388924e-7]
[-4.832950427206354, -10.243973454828646, -6.374480510182532, 0.9124536344093833, -1.054081113773402, -0.0027764381243728445, -27.099417689173645, -24.12014642586469, -14.154823892445176, 3.123448249887229 … 0.002091503073580616, 0.004790495824728955, 0.0064130548505908665, -1.0097282662914537e-5, 0.002954973888973382, 0.0027142818106974755, 0.0030527928355762075, -0.002221124201936867, 0.0035754220503216583, 3.1303036697196277e-6]
[-5.301834487820669, -10.043425570222864, -6.044470356813267, 1.6324696388286903, 0.014754094034316748, -0.004469124405445839, -27.567160549523994, -23.571253807794363, -13.366738532818701, 4.7624850496591105 … 0.0022790096060926733, 0.004264275065985961, 0.00723851734414768, -1.0748157585763421e-5, 0.003000100652640787, 0.0026542211850299263, 0.002876953390018217, -0.0033188564871101797, -0.0005605108121298306, 7.403502171637199e-6]
[-5.780992436524783, -9.820142737820397, -5.670384330893398, 2.295996439740908, 1.1131630735447262, -0.00619414408504289, -28.011954137560274, -22.965970847707634, -12.478927717766307, 6.219458924468294 … 0.002460705588213133, 0.0034793543138389787, 0.005757881090494939, -9.123150370842618e-6, 0.0030420717113689094, 0.0025878234238336426, 0.0026806441180786695, -0.004321542545469004, -0.0047199182666720595, 1.1673193478401132e-5]
⋮
[-15.123504880335279, 7.436746717146681, 7.0761907396791255, -0.3237599663918287, -0.24463120803472516, -0.24764237032684164, -14.882701102443601, 19.539513342467565, 16.072033745949224, -1.6713112197225537 … -0.001160251667326252, -0.005128909514173758, 0.013406024965720098, -1.3826895772078644e-5, 0.0017136059028290632, -0.0022727749008594126, -0.003806445734746988, 0.00035827273609225163, -0.0008070565008280845, 7.136130605900503e-6]
[-15.113901419297289, 7.499058102575787, 7.043085974149579, -0.4753174980723805, 0.12498254477801703, -0.2480407769352984, -14.7352278407547, 19.68190731669791, 15.987794510743505, -2.034359619653406 … -0.0012126206126505483, -0.0050962296801012425, 0.00903042241181972, -9.655371622955402e-6, 0.0017024396342444237, -0.0022891500359236067, -0.003790064818694251, 0.000585847255365571, -0.010752958467877403, 1.6567423276266954e-5]
[-15.103310765058385, 7.564866246309667, 7.006065780345468, -0.6361170109102849, 0.2848818208407543, -0.24824187589935637, -14.577984093448466, 19.832138852478565, 15.89410396364736, -2.41768811362579 … -0.0012681468530740126, -0.005052225319344712, -0.0012080846199018069, 1.1090198599622733e-7, 0.0016905614536774445, -0.0023064279701084, -0.003771662221979792, 0.0008251948984017012, -0.014709895113032258, 2.0277165409791478e-5]
[-15.093212592225893, 7.6250593243969975, 6.970315552665813, -0.7836254589316607, 0.12588778459639344, -0.24813478550939833, -14.432779380897545, 19.969411068950667, 15.804072814445728, -2.7676429555305884 … -0.0013191256148484029, -0.005003424009126401, -0.010367586697101, 8.845330136683076e-6, 0.0016796179901478485, -0.0023222166539209953, -0.0037538170305137708, 0.001042869233914279, -0.009867077739329011, 1.559121167009565e-5]
[-15.081964639971586, 7.689460127250504, 6.930026130572176, -0.9416427549587274, -0.27744312009764, -0.24779806463362763, -14.275929480639178, 20.116126762395417, 15.703067756145312, -3.1407194351246654 … -0.0013738610165446427, -0.004942165596253046, -0.013701465533688091, 1.2014244382035964e-5, 0.001667823942347297, -0.0023390928066408002, -0.003733628716205003, 0.0012740708038924335, 0.0017079480543546791, 4.473568017181597e-6]
[-15.071696732554562, 7.74603836925127, 6.892855000957991, -1.0804215728216584, -0.6263998108018394, -0.24750791391709984, -14.136824853947955, 20.24488853101855, 15.61025945652639, -3.4668080997230453 … -0.001422106671047666, -0.00488065320257526, -0.008934450850791283, 7.448130599511361e-6, 0.0016573875617627715, -0.002353905038257173, -0.003714936701361599, 0.0014754344430469447, 0.01135485828477465, -4.7946689575205085e-6]
[-15.060087018082095, 7.807696777437353, 6.850415937231508, -1.2313799093871334, -0.7733415861109256, -0.24741515562963104, -13.983802629936768, 20.385067145201305, 15.504692395290489, -3.8198195638342836 … -0.0014748476318909832, -0.004805441643873802, 0.0015099929442110973, -2.542313377690709e-6, 0.0016459321000008114, -0.002370031924568646, -0.0036935253632895794, 0.001692653966850431, 0.014830260403125675, -8.17504995263642e-6]
[-15.048711768820109, 7.8659306857835585, 6.808447606661263, -1.373462790665074, -0.6070450696230812, -0.2476195329514708, -13.83787783371418, 20.51731859222903, 15.400665353234457, -4.150423249259988 … -0.0015248090874105844, -0.004726590470355707, 0.010466680863388722, -1.1113475172623207e-5, 0.0016350322014704528, -0.0023852479540122937, -0.0036722851602993974, 0.0018953701844962368, 0.009529652186270844, -3.172166109151845e-6]
[-15.043661195419459, 7.891142465336795, 6.7896998135214375, -1.4347720735878677, -0.46136831765110276, -0.24777859657528684, -13.774270020226949, 20.574532028027424, 15.354303509533002, -4.292573293212216 … -0.0015464825683999225, -0.004690103854536925, 0.012686428866911669, -1.3241951067406622e-5, 0.0016302883546701367, -0.002391830913303292, -0.0036627771872230768, 0.0019823252000823006, 0.005302535202939572, 8.38894652418241e-7]plt = plot(xlab = "x", ylab = "y", zlab = "z", title = "Outer solar system")
for i in 1:N
plot!(plt, sol, idxs = (u[:, i]...,), lab = planets[i])
end
plt