GitHub

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}$.

planetmassinitial positioninitial 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.

using Plots, OrdinaryDiffEq, ModelingToolkit
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 t u(t)[1:3, 1:N]
u = collect(u)
D = Differential(t)
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{7.6924 \cdot 10^{-9}}{\sqrt{\left( - u(t)_{1}ˏ_1 + u(t)_{1}ˏ_6 \right)^{2} + \left( - u(t)_{2}ˏ_1 + u(t)_{2}ˏ_6 \right)^{2} + \left( - u(t)_{3}ˏ_1 + u(t)_{3}ˏ_6 \right)^{2}}} + \frac{2.1968 \cdot 10^{-12}}{\sqrt{\left( - u(t)_{1}ˏ_3 + u(t)_{1}ˏ_6 \right)^{2} + \left( - u(t)_{2}ˏ_3 + u(t)_{2}ˏ_6 \right)^{2} + \left( - u(t)_{3}ˏ_3 + u(t)_{3}ˏ_6 \right)^{2}}} + \frac{0.00028559}{\sqrt{\left( - u(t)_{1}ˏ_1 + u(t)_{1}ˏ_3 \right)^{2} + \left( - u(t)_{2}ˏ_1 + u(t)_{2}ˏ_3 \right)^{2} + \left( - u(t)_{3}ˏ_1 + u(t)_{3}ˏ_3 \right)^{2}}} + \frac{4.3728 \cdot 10^{-5}}{\sqrt{\left( - u(t)_{1}ˏ_1 + u(t)_{1}ˏ_4 \right)^{2} + \left( - u(t)_{2}ˏ_1 + u(t)_{2}ˏ_4 \right)^{2} + \left( - u(t)_{3}ˏ_1 + u(t)_{3}ˏ_4 \right)^{2}}} + \frac{0.00095479}{\sqrt{\left( - u(t)_{1}ˏ_1 + u(t)_{1}ˏ_2 \right)^{2} + \left( - u(t)_{2}ˏ_1 + u(t)_{2}ˏ_2 \right)^{2} + \left( - u(t)_{3}ˏ_1 + u(t)_{3}ˏ_2 \right)^{2}}} + \frac{1.4786 \cdot 10^{-8}}{\sqrt{\left( - u(t)_{1}ˏ_3 + u(t)_{1}ˏ_5 \right)^{2} + \left( - u(t)_{2}ˏ_3 + u(t)_{2}ˏ_5 \right)^{2} + \left( - u(t)_{3}ˏ_3 + u(t)_{3}ˏ_5 \right)^{2}}} + \frac{7.3445 \cdot 10^{-12}}{\sqrt{\left( - u(t)_{1}ˏ_2 + u(t)_{1}ˏ_6 \right)^{2} + \left( - u(t)_{2}ˏ_2 + u(t)_{2}ˏ_6 \right)^{2} + \left( - u(t)_{3}ˏ_2 + u(t)_{3}ˏ_6 \right)^{2}}} + \frac{1.2488 \cdot 10^{-8}}{\sqrt{\left( - u(t)_{1}ˏ_3 + u(t)_{1}ˏ_4 \right)^{2} + \left( - u(t)_{2}ˏ_3 + u(t)_{2}ˏ_4 \right)^{2} + \left( - u(t)_{3}ˏ_3 + u(t)_{3}ˏ_4 \right)^{2}}} + \frac{2.7267 \cdot 10^{-7}}{\sqrt{\left( - u(t)_{1}ˏ_2 + u(t)_{1}ˏ_3 \right)^{2} + \left( - u(t)_{2}ˏ_2 + u(t)_{2}ˏ_3 \right)^{2} + \left( - u(t)_{3}ˏ_2 + u(t)_{3}ˏ_3 \right)^{2}}} + \frac{5.1776 \cdot 10^{-5}}{\sqrt{\left( - u(t)_{1}ˏ_1 + u(t)_{1}ˏ_5 \right)^{2} + \left( - u(t)_{2}ˏ_1 + u(t)_{2}ˏ_5 \right)^{2} + \left( - u(t)_{3}ˏ_1 + u(t)_{3}ˏ_5 \right)^{2}}} + \frac{2.264 \cdot 10^{-9}}{\sqrt{\left( - u(t)_{1}ˏ_4 + u(t)_{1}ˏ_5 \right)^{2} + \left( - u(t)_{2}ˏ_4 + u(t)_{2}ˏ_5 \right)^{2} + \left( - u(t)_{3}ˏ_4 + u(t)_{3}ˏ_5 \right)^{2}}} + \frac{4.175 \cdot 10^{-8}}{\sqrt{\left( - u(t)_{1}ˏ_2 + u(t)_{1}ˏ_4 \right)^{2} + \left( - u(t)_{2}ˏ_2 + u(t)_{2}ˏ_4 \right)^{2} + \left( - u(t)_{3}ˏ_2 + u(t)_{3}ˏ_4 \right)^{2}}} + \frac{3.3636 \cdot 10^{-13}}{\sqrt{\left( - u(t)_{1}ˏ_4 + u(t)_{1}ˏ_6 \right)^{2} + \left( - u(t)_{2}ˏ_4 + u(t)_{2}ˏ_6 \right)^{2} + \left( - u(t)_{3}ˏ_4 + u(t)_{3}ˏ_6 \right)^{2}}} + \frac{4.9435 \cdot 10^{-8}}{\sqrt{\left( - u(t)_{1}ˏ_2 + u(t)_{1}ˏ_5 \right)^{2} + \left( - u(t)_{2}ˏ_2 + u(t)_{2}ˏ_5 \right)^{2} + \left( - u(t)_{3}ˏ_2 + u(t)_{3}ˏ_5 \right)^{2}}} + \frac{3.9828 \cdot 10^{-13}}{\sqrt{\left( - u(t)_{1}ˏ_5 + u(t)_{1}ˏ_6 \right)^{2} + \left( - u(t)_{2}ˏ_5 + u(t)_{2}ˏ_6 \right)^{2} + \left( - u(t)_{3}ˏ_5 + u(t)_{3}ˏ_6 \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))) .~ .-ModelingToolkit.gradient(potential, vec(u)) ./
                         repeat(M, inner = 3)
@named sys = ODESystem(eqs, t)
ss = structural_simplify(sys)
prob = ODEProblem(ss, [vec(u .=> pos); vec(D.(u) .=> vel)], tspan)
sol = solve(prob, Tsit5());
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 862-element Vector{Float64}:
      0.0
      0.4345623854592078
      4.780186240051286
     48.23642478597206
    195.2733726487084
    404.37878073611245
    641.246171048518
    917.4023638439503
   1259.3749858717613
   1620.723435061644
      ⋮
 199469.90760955325
 199550.3013901666
 199620.52093184175
 199697.8876414425
 199765.85003240773
 199832.37074896548
 199887.0525301258
 199959.31627976737
 200000.0
u: 862-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.00565429, -0.0041249, -0.00190589, 0.00168318, 0.00483525, 0.00192462, 0.00354178  …  -1.6483708, 8.310142, -16.2901086, -7.2521278, 11.4707666, -25.7294829, -10.8169456, -15.5387357, -25.2225594, -3.1902382]
 [-2.3461051313643855e-9, -3.1048711440018414e-9, -1.2785303085738863e-9, 0.005657137663985553, -0.004121794983806699, -0.0019046284357956216, 0.0016819059109041002, 0.004835677090231295, 0.0019248511933338966, 0.0035416394156575887  …  -1.6475343823037847, 8.31168109382013, -16.289512746418094, -7.251888638015288, 11.472022169413352, -25.728985182525822, -10.816773167663484, -15.537533140796315, -25.223301179998256, -3.190831391666314]
 [-2.570694151869523e-8, -3.423355432416456e-8, -1.4100666088910388e-8, 0.005685511868486027, -0.004090656421676488, -0.0018919724055948748, 0.0016691568418040644, 0.004839932607591813, 0.0019271571241916761, 0.003540232218112257  …  -1.6391646887336664, 8.327068669865051, -16.283547627110032, -7.249494087175589, 11.48457657746139, -25.724005124769256, -10.815047632292558, -15.52550574063629, -25.230716044387762, -3.196762937128868]
 [-2.4909025937027336e-7, -3.533537351731393e-7, -1.4594000158084998e-7, 0.0059587932915218115, -0.003770599600091882, -0.0017614403625624616, 0.0015408613020808024, 0.004880940100833716, 0.0019496118325323336, 0.003526025124018189  …  -1.554925794138472, 8.480606040821046, -16.22323923167503, -7.225255950862118, 11.609991524392697, -25.673916441625973, -10.797671145111675, -15.405051486626984, -25.304570881838067, -3.2560410257920296]
 [-8.562666124366701e-7, -1.5291383285282015e-6, -6.369559468830282e-7, 0.006730897179308313, -0.002580656126047471, -0.0012701900733253118, 0.0010966413711339354, 0.004998334414871359, 0.0020172033027033515, 0.003476148756678132  …  -1.263130363916593, 8.995429151468823, -16.01038370336294, -7.139325016951497, 12.032571237495207, -25.50056426465599, -10.737247076745723, -14.995086987900272, -25.5504886404068, -3.4560950369067904]
 [-1.270060008529579e-6, -3.3960857323731713e-6, -1.4300961196408601e-6, 0.007365006196212584, -0.000666459588156778, -0.0004651305514942468, 0.0004422141061446354, 0.005105982369438332, 0.002089810001970621, 0.003400506948947264  …  -0.8332345946735483, 9.714497324515916, -15.684637553718064, -7.006844903981199, 12.62867344820616, -25.243795602764315, -10.647004324035025, -14.405852178332685, -25.889578358759636, -3.7391502987570577]
 [-1.007513458232923e-6, -5.631384102662558e-6, -2.3985120866280533e-6, 0.007324487954155879, 0.0016548597238463258, 0.0005308765916475159, -0.0003215351780469653, 0.005139085078150229, 0.002136329561824399, 0.003308336255405168  …  -0.33194282247033946, 10.509184073050166, -15.283739886453963, -6.842521993262041, 13.296698993611704, -24.938543983852917, -10.538709853430426, -13.729932103781781, -26.25846459707381, -4.057582824414865]
 [3.554667634218679e-7, -8.082240490859979e-6, -3.48739746950257e-6, 0.006174492455065497, 0.004236423620030281, 0.00166546646021687, -0.0012242830831770985, 0.005052449172558458, 0.002139373108865423, 0.0031925159343643857  …  0.25974229996959247, 11.407012060467682, -14.774890228539132, -6.6323806318048275, 14.065323791019793, -24.56358392094972, -10.404390712741327, -12.931149424634967, -26.667923778020157, -4.425639184967824]
 [3.4901811625643595e-6, -1.0242904678740642e-5, -4.4970458226438195e-6, 0.003230726642949325, 0.006575790349577193, 0.002739938161380959, -0.002325513928515102, 0.00475206744994242, 0.002062677619382529, 0.0030372322298507144  …  0.980818022484472, 12.472577563460694, -14.085404632753704, -6.345502797745267, 15.000995668497582, -24.071220490411154, -10.2261812806247, -11.927011166605329, -27.143912501455333, -4.876236637628356]
 [7.821334963992643e-6, -1.0696453179221607e-5, -4.805135026272916e-6, -0.0009678260792008035, 0.007201496199593491, 0.003110415487860888, -0.0034164816398609067, 0.004202822632183789, 0.0018827832856910369, 0.002859724549035003  …  1.696709357432359, 13.538407580658335, -13.289128533479984, -6.0118582932866556, 15.969038388790834, -23.51787256963099, -10.023818114603046, -10.849436012073225, -27.609190660804597, -5.345590705748484]
 ⋮
 [-1.3756899865923538e-6, -1.1856988849768463e-5, -5.093599926230342e-6, 0.008108466947515833, 0.011342878131347564, 0.004663313543039903, 0.0003563582555440711, -0.005129121965515959, -0.002136109117456437, -0.0038064753898782414  …  -0.322525095220091, -3.8008812735055177, 16.07219105571034, 7.076252247480188, 21.9120036360489, 19.53924249954513, 7.436628206296885, 35.677707490565055, -14.882979199390995, -15.123522890358654]
 [1.0290969703860403e-5, -1.3299595876632719e-5, -5.992766866093057e-6, -0.00418761356476289, 0.012846269782764255, 0.005605335920948736, 0.0006138768679876507, -0.005091557164576367, -0.0021317369675808973, -0.0037878631617113342  …  -0.49411027036413474, -4.106157704287349, 15.97653753859869, 7.038648508503974, 21.72854313339182, 19.70039097935498, 7.507150978766658, 35.81496421502917, -14.715972962167172, -15.112624156836025]
 [1.8422189507319665e-5, -7.293784313068863e-6, -3.6152132445151746e-6, -0.01277034295586152, 0.006546792017930737, 0.0031142894218963196, 0.00083616482052726, -0.0050499312036854285, -0.0021241589762742074, -0.0037706736461482937  …  -0.6435541355611202, -4.371541288303119, 15.889096388256299, 7.004081515381987, 21.56708652511459, 19.839950113028596, 7.568289795865951, 35.934018796395165, -14.56975994173848, -15.102747159638314]
 [1.9528386162006377e-5, 3.3262939514657635e-6, 9.096610556033493e-7, -0.014001168100877622, -0.0045892138624837545, -0.001627061546858441, 0.0010778010581308781, -0.00499476579022963, -0.002111821183346825, -0.0037507284954909864  …  -0.8074427869775495, -4.662501148922091, 15.788567292840264, 6.964142414016437, 21.387896914587838, 19.992413987860015, 7.635151793254869, 36.06429336395901, -14.40830086790954, -15.091480560674173]
 [1.2729354369534523e-5, 1.0623232358334157e-5, 4.2020995520387704e-6, -0.006942686061080501, -0.01224547384816105, -0.005078357334930667, 0.0012869121964130942, -0.004938439642149392, -0.0020975939409277737, -0.0037323393617425837  …  -0.9505015180693351, -4.916789336755138, 15.696648826990723, 6.927459569145205, 21.229373931097484, 20.12521273901333, 7.693450676653517, 36.17795589243498, -14.26615755047737, -15.081252479171315]
 [2.33239946080059e-6, 1.1631350388915576e-5, 4.887271482087873e-6, 0.0038863392921218965, -0.01331685726214166, -0.005799741699513141, 0.0014884671932481634, -0.0048763603846162795, -0.002080660738285056, -0.0037135545994000404  …  -1.0894880821004658, -5.164446745686518, 15.603424030605593, 6.890111835349553, 21.07321224696216, 20.25416191782829, 7.750115212924507, 36.288504097517475, -14.126752292520793, -15.070942830214259]
 [-4.998642334648639e-6, 7.266766578127084e-6, 3.1960207606871922e-6, 0.011515699408678384, -0.008759822883215505, -0.004032137381178998, 0.0016516585770397461, -0.0048203105659747345, -0.0020645580826820493, -0.0036975312604936453  …  -1.2028309419284258, -5.367074823899943, 15.524387464479677, 6.85834537275782, 20.94410683799006, 20.359392049171102, 7.796398351643997, 36.378856724134025, -14.01195473852491, -15.062247693208693]
 [-8.196754870727107e-6, -2.6748389504926338e-6, -9.844384223416686e-7, 0.014801734862602478, 0.0016317523766605242, 0.00034047591507541987, 0.0018636163158711409, -0.004739466781581778, -0.002040317084686035, -0.0036755520924993228  …  -1.3511678877199182, -5.633483652750338, 15.416624935179177, 6.814895102491169, 20.772478417684624, 20.49738562821115, 7.8571496785103365, 36.4975397492519, -13.859969063297443, -15.050452992615021]
 [-6.075512849690596e-6, -8.147820596725328e-6, -3.3798588857354807e-6, 0.01254481662737608, 0.007351192486272225, 0.002845623313972247, 0.0019809916131378175, -0.004690633663013654, -0.002025212599964893, -0.0036627763874156707  …  -1.4338718456440076, -5.782759883840572, 15.354302296423441, 6.789699152830791, 20.67535024070022, 20.574535098139414, 7.891143772620107, 36.5639959717329, -13.77426559411742, -15.043661040103219]
plt = plot()
for i in 1:N
    plot!(plt, sol, idxs = (u[:, i]...,), lab = planets[i])
end
plot!(plt; xlab = "x", ylab = "y", zlab = "z", title = "Outer solar system")
Example block output