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\left( t \right)_{1,1} + u\left( t \right)_{1,6} \right)^{2} + \left( - u\left( t \right)_{2,1} + u\left( t \right)_{2,6} \right)^{2} + \left( - u\left( t \right)_{3,1} + u\left( t \right)_{3,6} \right)^{2}}} + \frac{2.1968 \cdot 10^{-12}}{\sqrt{\left( - u\left( t \right)_{1,3} + u\left( t \right)_{1,6} \right)^{2} + \left( - u\left( t \right)_{2,3} + u\left( t \right)_{2,6} \right)^{2} + \left( - u\left( t \right)_{3,3} + u\left( t \right)_{3,6} \right)^{2}}} + \frac{0.00028559}{\sqrt{\left( - u\left( t \right)_{1,1} + u\left( t \right)_{1,3} \right)^{2} + \left( - u\left( t \right)_{2,1} + u\left( t \right)_{2,3} \right)^{2} + \left( - u\left( t \right)_{3,1} + u\left( t \right)_{3,3} \right)^{2}}} + \frac{4.3728 \cdot 10^{-5}}{\sqrt{\left( - u\left( t \right)_{1,1} + u\left( t \right)_{1,4} \right)^{2} + \left( - u\left( t \right)_{2,1} + u\left( t \right)_{2,4} \right)^{2} + \left( - u\left( t \right)_{3,1} + u\left( t \right)_{3,4} \right)^{2}}} + \frac{0.00095479}{\sqrt{\left( - u\left( t \right)_{1,1} + u\left( t \right)_{1,2} \right)^{2} + \left( - u\left( t \right)_{2,1} + u\left( t \right)_{2,2} \right)^{2} + \left( - u\left( t \right)_{3,1} + u\left( t \right)_{3,2} \right)^{2}}} + \frac{1.4786 \cdot 10^{-8}}{\sqrt{\left( - u\left( t \right)_{1,3} + u\left( t \right)_{1,5} \right)^{2} + \left( - u\left( t \right)_{2,3} + u\left( t \right)_{2,5} \right)^{2} + \left( - u\left( t \right)_{3,3} + u\left( t \right)_{3,5} \right)^{2}}} + \frac{7.3445 \cdot 10^{-12}}{\sqrt{\left( - u\left( t \right)_{1,2} + u\left( t \right)_{1,6} \right)^{2} + \left( - u\left( t \right)_{2,2} + u\left( t \right)_{2,6} \right)^{2} + \left( - u\left( t \right)_{3,2} + u\left( t \right)_{3,6} \right)^{2}}} + \frac{1.2488 \cdot 10^{-8}}{\sqrt{\left( - u\left( t \right)_{1,3} + u\left( t \right)_{1,4} \right)^{2} + \left( - u\left( t \right)_{2,3} + u\left( t \right)_{2,4} \right)^{2} + \left( - u\left( t \right)_{3,3} + u\left( t \right)_{3,4} \right)^{2}}} + \frac{2.7267 \cdot 10^{-7}}{\sqrt{\left( - u\left( t \right)_{1,2} + u\left( t \right)_{1,3} \right)^{2} + \left( - u\left( t \right)_{2,2} + u\left( t \right)_{2,3} \right)^{2} + \left( - u\left( t \right)_{3,2} + u\left( t \right)_{3,3} \right)^{2}}} + \frac{5.1776 \cdot 10^{-5}}{\sqrt{\left( - u\left( t \right)_{1,1} + u\left( t \right)_{1,5} \right)^{2} + \left( - u\left( t \right)_{2,1} + u\left( t \right)_{2,5} \right)^{2} + \left( - u\left( t \right)_{3,1} + u\left( t \right)_{3,5} \right)^{2}}} + \frac{2.264 \cdot 10^{-9}}{\sqrt{\left( - u\left( t \right)_{1,4} + u\left( t \right)_{1,5} \right)^{2} + \left( - u\left( t \right)_{2,4} + u\left( t \right)_{2,5} \right)^{2} + \left( - u\left( t \right)_{3,4} + u\left( t \right)_{3,5} \right)^{2}}} + \frac{4.175 \cdot 10^{-8}}{\sqrt{\left( - u\left( t \right)_{1,2} + u\left( t \right)_{1,4} \right)^{2} + \left( - u\left( t \right)_{2,2} + u\left( t \right)_{2,4} \right)^{2} + \left( - u\left( t \right)_{3,2} + u\left( t \right)_{3,4} \right)^{2}}} + \frac{3.3636 \cdot 10^{-13}}{\sqrt{\left( - u\left( t \right)_{1,4} + u\left( t \right)_{1,6} \right)^{2} + \left( - u\left( t \right)_{2,4} + u\left( t \right)_{2,6} \right)^{2} + \left( - u\left( t \right)_{3,4} + u\left( t \right)_{3,6} \right)^{2}}} + \frac{4.9435 \cdot 10^{-8}}{\sqrt{\left( - u\left( t \right)_{1,2} + u\left( t \right)_{1,5} \right)^{2} + \left( - u\left( t \right)_{2,2} + u\left( t \right)_{2,5} \right)^{2} + \left( - u\left( t \right)_{3,2} + u\left( t \right)_{3,5} \right)^{2}}} + \frac{3.9828 \cdot 10^{-13}}{\sqrt{\left( - u\left( t \right)_{1,5} + u\left( t \right)_{1,6} \right)^{2} + \left( - u\left( t \right)_{2,5} + u\left( t \right)_{2,6} \right)^{2} + \left( - u\left( t \right)_{3,5} + u\left( t \right)_{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.27337773985937
404.37880043365203
641.2461951827266
917.4024243418721
1259.375074434556
1620.723519074113
⋮
199469.8999485139
199550.29373994426
199620.51324402454
199697.8799456865
199765.84231811052
199832.3629902858
199887.04472779683
199959.30845589447
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.562666292142268e-7, -1.5291383715369453e-6, -6.369559649765034e-7, 0.006730897201673277, -0.0025806560822637585, -0.001270190055102589, 0.001096641355501883, 0.004998334418352889, 0.002017203304813606, 0.0034761487549033287 … -1.2631303536467062, 8.995429169166421, -16.010383695759355, -7.139325013872073, 12.032571252078975, -25.500564258550423, -10.737247074610087, -14.995086973641925, -25.55048864881513, -3.4560950438194866]
[-1.2700600202222442e-6, -3.3960859156204874e-6, -1.4300961982037678e-6, 0.007365006227690945, -0.0006664593990847896, -0.000465130471216582, 0.0004422140434641477, 0.005105982376182712, 0.00208981000745197, 0.0034005069415658402 … -0.8332345535094318, 9.714497391497536, -15.68463752177678, -7.006844890940794, 12.628673504082018, -25.24379557801365, -10.647004315296812, -14.405852122492028, -25.88957839011003, -3.7391503253370364]
[-1.0075133892324704e-6, -5.6313843294952595e-6, -2.3985121859603363e-6, 0.007324487906032101, 0.0016548599609000753, 0.0005308766944316987, -0.00032153525665023434, 0.005139085076549321, 0.002136329564543771, 0.0033083362456712614 … -0.33194277091171415, 10.509184152894242, -15.283739843913043, -6.842521975761539, 13.296699061272498, -24.938543951977763, -10.53870984206983, -13.729932034469996, -26.25846463383041, -4.057582856734454]
[3.554671854327567e-7, -8.082240979433938e-6, -3.4873976903991466e-6, 0.006174492073878075, 0.0042364241395163585, 0.0016654666921783507, -0.0012242832807843097, 0.005052449138437194, 0.0021393731032716677, 0.0031925159080346227 … 0.25974242939723763, 11.40701225360826, -14.77489011229215, -6.632380583628376, 14.06532395815657, -24.563583836576818, -10.404390682372131, -12.931149248426918, -26.667923865275117, -4.425639265202051]
[3.490182139991148e-6, -1.0242905047747281e-5, -4.497046006563054e-6, 0.0032307257053528857, 0.006575790764247339, 0.0027399383619676563, -0.0023255142074334366, 0.004752067344284884, 0.002062677587743199, 0.0030372321879999563 … 0.9808182051610246, 12.472577832446467, -14.085404445957762, -6.345502719745113, 15.000995908408392, -24.071220358929434, -10.226181232787367, -11.927010904521099, -27.143912620248173, -4.876236753538832]
[7.821335996287877e-6, -1.0696453026522622e-5, -4.805134987913176e-6, -0.0009678270857685672, 0.007201496082832918, 0.0031104154623282946, -0.003416481880143624, 0.004202822477539218, 0.001882783232168296, 0.0028597245062313105 … 1.6967095156096763, 13.538407820910853, -13.289128340584687, -6.011858212208256, 15.969038611287116, -23.517872437082897, -10.023818065895448, -10.849435759690564, -27.60919076445552, -5.345590814017202]
⋮
[-1.3756744229968758e-6, -1.1856993777287261e-5, -5.093602561207465e-6, 0.008108457980597006, 0.011342883850856458, 0.004663316210511856, 0.00035633367830419674, -0.0051291250113228055, -0.0021361093115094673, -0.003806477109536711 … -0.32250879505040436, -3.8008521073327577, 16.072199944917653, 7.07625573090995, 21.912021048390184, 19.539227073137685, 7.436621459039603, 35.677694362158256, -14.882995094372687, -15.123523908053537]
[1.029098928373866e-5, -1.329959190658e-5, -5.99276578936649e-6, -0.004187626827607001, 0.012846266491634149, 0.005605334832474342, 0.0006138526092754943, -0.005091561239020404, -0.0021317376007227357, -0.003787864987417807 … -0.4940940265906956, -4.106128721704839, 15.97654686887402, 7.038652187625974, 21.728560662928654, 19.700375713235523, 7.50714429419396, 35.81495120188578, -14.715988874803228, -15.112625214771747]
[1.8422201038654874e-5, -7.293776761652931e-6, -3.615210441504003e-6, -0.012770347834418713, 0.006546785245580066, 0.003114286639018903, 0.0008361407306580451, -0.005049936185066453, -0.0021241599917019606, -0.0037706755759733947 … -0.6435378696730796, -4.371512295487581, 15.889106161162271, 7.0040853881646274, 21.567104264190675, 19.839934894373435, 7.568283125502167, 35.93400580419888, -14.56977596724869, -15.102748259192193]
[1.952839074389161e-5, 3.326302749363971e-6, 9.096645589644748e-7, -0.014001165795204574, -0.004589221653123887, -0.0016270649402344186, 0.0010777773025554135, -0.004994771732040524, -0.0021118226103461753, -0.0037507305321528685 … -0.807426598801957, -4.662472279686814, 15.788577511855342, 6.964146483858048, 21.387914806941232, 19.992398889551627, 7.635145168330223, 36.06428045197091, -14.40831694766038, -15.09148170135516]
[1.2729355534881039e-5, 1.0623237014941348e-5, 4.202101362366456e-6, -0.006942680259770478, -0.012245477052168007, -0.005078358848319472, 0.0012868887339745352, -0.004938446418789154, -0.0020975957257873433, -0.003732341495518283 … -0.9504854000366415, -4.916760539860718, 15.696659453069875, 6.927463818232227, 21.229391984401875, 20.125197724977966, 7.693444082212711, 36.177943031771214, -14.266173701679623, -15.081253657663263]
[2.3324058961651186e-6, 1.1631352420154572e-5, 4.887272038195817e-6, 0.0038863395071219568, -0.013316857470596516, -0.005799741794611661, 0.001488443970222344, -0.004876367988605522, -0.002080662875585741, -0.0037135568360866576 … -1.0894720018046966, -5.16441792891375, 15.603435092846192, 6.890116275130058, 21.07323051922373, 20.254146936940117, 7.750108626598513, 36.28849124400135, -14.12676856841137, -15.070944049864142]
[-4.998631870479622e-6, 7.2667732266931395e-6, 3.1960231959236945e-6, 0.011515695348352117, -0.008759827726036576, -0.004032139358511234, 0.0016516355539531228, -0.004820318849497776, -0.0020645605095986146, -0.0036975335845717537 … -1.2028148959554052, -5.367045970073011, 15.52439889782596, 6.858349974506756, 20.94412530750435, 20.359377083368916, 7.79639176658115, 36.37884386537011, -14.011971131918841, -15.062248948035348]
[-8.19674967899049e-6, -2.6748285239534262e-6, -9.844342421184133e-7, 0.014801736209212563, 0.001631743829273253, 0.0003404722196000947, 0.001863593695957123, -0.004739475910084152, -0.0020403198784278223, -0.0036755545219508625 … -1.351151986197135, -5.63345489125107, 15.416636807738904, 6.8148998978425785, 20.77249706241604, 20.49737075352989, 7.857143126378928, 36.49752694383905, -13.85998553588676, -15.050454288286229]
[-6.074855057339981e-6, -8.148753589265148e-6, -3.3802745576123506e-6, 0.012544127654405973, 0.007352169655363238, 0.002846058669027302, 0.0019809917026013344, -0.004690633625871623, -0.002025212588104413, -0.0036627763878147894 … -1.433871906851642, -5.782759879418573, 15.354302298141548, 6.789699153547345, 20.675350241507825, 20.574535097506825, 7.891143772347629, 36.563995971508824, -13.774265594556638, -15.043661040109278]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")