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.
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\_{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{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{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{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{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{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{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{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.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{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.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{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{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{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{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}}} \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.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}}:
[0.0, 0.00565429, 0.00168318, 0.00354178, 0.0028893, 0.00276725, 0.0, -0.0041249, 0.00483525, 0.00137102 … -1.6483708, 8.310142, -16.2901086, -7.2521278, 11.4707666, -25.7294829, -10.8169456, -15.5387357, -25.2225594, -3.1902382]
[-2.3461051313643834e-9, 0.005657137663985553, 0.0016819059109041002, 0.0035416394156575887, 0.0028892462117755886, 0.002767325672882192, -3.1048711440018414e-9, -0.004121794983806699, 0.004835677090231295, 0.0013712954942851254 … -1.6475343823037847, 8.31168109382013, -16.289512746418094, -7.251888638015288, 11.472022169413352, -25.728985182525822, -10.816773167663484, -15.537533140796315, -25.223301179998256, -3.190831391666314]
[-2.5706941518695228e-8, 0.005685511868486027, 0.0016691568418040644, 0.003540232218112257, 0.0028887080028757346, 0.002768082053772159, -3.4233554324164556e-8, -0.004090656421676488, 0.004839932607591813, 0.0013740497263540709 … -1.6391646887336664, 8.327068669865051, -16.283547627110032, -7.249494087175589, 11.48457657746139, -25.724005124769256, -10.815047632292558, -15.52550574063629, -25.230716044387762, -3.196762937128868]
[-2.4909025937027336e-7, 0.0059587932915218115, 0.0015408613020808024, 0.003526025124018189, 0.002883293267559013, 0.0027756110286434707, -3.53353735173139e-7, -0.003770599600091882, 0.004880940100833716, 0.0014015206514237767 … -1.554925794138472, 8.480606040821046, -16.22323923167503, -7.225255950862118, 11.609991524392697, -25.673916441625973, -10.797671145111675, -15.405051486626984, -25.304570881838067, -3.2560410257920296]
[-8.562666292142268e-7, 0.006730897201673277, 0.001096641355501883, 0.0034761487549033287, 0.0028645330514698616, 0.002800613918439817, -1.529138371536945e-6, -0.0025806560822637585, 0.004998334418352889, 0.0014934905193989781 … -1.2631303536467064, 8.995429169166421, -16.010383695759355, -7.139325013872073, 12.032571252078975, -25.500564258550423, -10.737247074610087, -14.995086973641925, -25.55048864881513, -3.4560950438194866]
[-1.2700600202222442e-6, 0.007365006227690945, 0.00044221404346414734, 0.0034005069415658402, 0.002836692350484805, 0.002834905191574055, -3.396085915620486e-6, -0.0006664593990847885, 0.005105982376182712, 0.0016215875225696298 … -0.833234553509432, 9.714497391497536, -15.68463752177678, -7.006844890940794, 12.628673504082018, -25.24379557801365, -10.647004315296812, -14.405852122492028, -25.88957839011003, -3.7391503253370364]
[-1.0075133892324712e-6, 0.007324487906032102, -0.00032153525665023434, 0.0033083362456712614, 0.0028035224672855984, 0.002871931193437516, -5.631384329495253e-6, 0.0016548599609000731, 0.005139085076549321, 0.0017626814231114039 … -0.33194277091171437, 10.509184152894242, -15.283739843913043, -6.842521975761539, 13.296699061272498, -24.938543951977763, -10.53870984206983, -13.729932034469996, -26.25846463383041, -4.057582856734454]
[3.554671854327595e-7, 0.006174492073878071, -0.00122428328078431, 0.0031925159080346227, 0.0027626862901216956, 0.0029126298323277314, -8.082240979433938e-6, 0.004236424139516364, 0.005052449138437194, 0.0019215037459887952 … 0.2597424293972374, 11.40701225360826, -14.77489011229215, -6.632380583628376, 14.06532395815657, -24.563583836576818, -10.404390682372131, -12.931149248426918, -26.667923865275117, -4.425639265202051]
[3.4901821399911367e-6, 0.003230725705352898, -0.002325514207433437, 0.0030372321879999563, 0.0027089345099226463, 0.002959304010544335, -1.0242905047747271e-5, 0.006575790764247332, 0.004752067344284884, 0.0021091920800701013 … 0.9808182051610244, 12.472577832446465, -14.085404445957762, -6.345502719745113, 15.000995908408392, -24.071220358929434, -10.226181232787367, -11.927010904521099, -27.143912620248173, -4.876236753538832]
[7.821335996287854e-6, -0.0009678270857685428, -0.0034164818801436228, 0.0028597245062313105, 0.002648372176935374, 0.0030041096081688755, -1.069645302652259e-5, 0.007201496082832887, 0.004202822477539219, 0.002296031801046592 … 1.6967095156096763, 13.538407820910852, -13.289128340584687, -6.011858212208256, 15.969038611287116, -23.517872437082897, -10.023818065895448, -10.849435759690564, -27.60919076445552, -5.345590814017202]
⋮
[-1.3756693948065007e-6, 0.008108454173305517, 0.0003563287887563882, -0.003806477450287144, -0.0022727417098904676, 0.0017136284822634343, -1.185699607645447e-5, 0.01134288637029172, -0.005129125617158283, -0.0011601455648747171 … -0.3225055526973783, -3.80084632786859, 16.07220170640706, 7.076256421183583, 21.912024499155155, 19.539224015936842, 7.436620121868759, 35.67769176029667, -14.882998244513848, -15.123524109751513]
[1.0290996628378757e-5, -0.004187633084377313, 0.0006138478015937808, -0.003787865347792654, -0.002291273610594168, 0.0017009850441485415, -1.329959030311957e-5, 0.012846264984274715, -0.0050915620463481874, -0.0012194307251699133 … -0.4940908078582244, -4.106123000682689, 15.976548710660907, 7.03865291387637, 21.728564123573904, 19.700372699422818, 7.507142974535007, 35.81494863277498, -14.715992016328638, -15.112625423639356]
[1.8422204096992478e-5, -0.012770349595343568, 0.0008361359044072283, -0.0037706759610563424, -0.0023073242169766445, 0.0016899424648723028, -7.293774303388143e-6, 0.006546782898839691, -0.0050499371828510616, -0.0012710345365739772 … -0.6435346114512911, -4.371506510012115, 15.88910811136718, 7.004086160982199, 21.56710780440406, 19.83993185715745, 7.568281794283585, 35.934003211244324, -14.569779165569413, -15.102748478644132]
[1.9528391601416645e-5, -0.014001165263383377, 0.001077772521183009, -0.0037507309404401132, -0.0023248604628924884, 0.0016777770342673783, 3.3263048234364494e-6, -0.004589223538620178, -0.004994772927719658, -0.0013276831227973659 … -0.8074233411188485, -4.662466492151749, 15.788579560538553, 6.964147299765075, 21.387918394305576, 19.99239586237846, 7.635143840047986, 36.06427786307255, -14.408320171689876, -15.091481930070932]
[1.2729356638888073e-5, -0.006942679994640394, 0.001286883984482558, -0.0037323419257314088, -0.002340135972218436, 0.0016670915913214366, 1.06232374929685e-5, -0.01224547721394641, -0.0049384477902885025, -0.0013772523961088111 … -0.9504821378662232, -4.916754733668444, 15.69666159559459, 6.927464674965965, 21.22939562482573, 20.125194697411878, 7.693442752448569, 36.177940438353644, -14.266176958634356, -15.081253895317323]
[2.332407338051759e-6, 0.0038863394000273107, 0.0014884392853463401, -0.003713557285506897, -0.0023549697738096983, 0.0016566341977526057, 1.1631352784062142e-5, -0.013316857464432426, -0.004876369522254775, -0.0014255822908936055 … -1.08946875848104, -5.164412138586617, 15.603437315682786, 6.890117167250405, 21.073234191194505, 20.254143926393027, 7.750107303014459, 36.28848866089309, -14.126771839281062, -15.070944294976355]
[-4.998629839283249e-6, 0.011515694611077484, 0.001651630915142788, -0.003697534050983003, -0.0023670759141878896, 0.0016480391398844246, 7.266774386690968e-6, -0.00875982851352661, -0.004820320518135741, -0.0014651640165654746 … -1.202811663571521, -5.367040179366535, 15.524401192430672, 6.858350898045243, 20.944129014588455, 20.359374079521185, 7.796390444864761, 36.37884128435904, -14.011974422384405, -15.062249199909228]
[-8.196748588569976e-6, 0.014801736432386369, 0.0018635891439395197, -0.0036755550089151356, -0.002382952482684855, 0.0016366823258962035, -2.674826944228626e-6, 0.0016317426524223222, -0.004739477746682838, -0.0015172592179612584 … -1.3511487868570342, -5.633449126112974, 15.416639187583709, 6.814900859061516, 20.7725008001034, 20.497367771613536, 7.857141812876675, 36.49752437667235, -13.859988838214676, -15.050454548041067]
[-6.0747229013532665e-6, 0.012543989239067051, 0.0019809917035306527, -0.003662776387803894, -0.002391831329951801, 0.0016302880834588448, -8.148941151253767e-6, 0.007352366100350984, -0.004690633625287779, -0.0015464832920667488 … -1.4338719077983646, -5.782759879477863, 15.354302298126019, 6.789699153540188, 20.675350241475257, 20.574535097530337, 7.891143772357861, 36.563995971518366, -13.77426559455976, -15.043661040110429]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")