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.8991880535
199550.29297594895
199620.512476788
199697.87917603034
199765.84155104912
199832.36223669595
199887.04398528096
199959.30772135654
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.3756745426261215e-6, -1.1856993370044154e-5, -5.093602398425872e-6, 0.008108458836404285, 0.011342883480072785, 0.004663316030889641, 0.0003563312302039047, -0.005129125314631053, -0.0021361093308668324, -0.003806477280192513 … -0.322507171720548, -3.8008492128113947, 16.072200827129425, 7.076256076621034, 21.91202277666881, 19.539225541961, 7.436620789327366, 35.67769305902247, -14.882996672105438, -15.123524009072867]
[1.0290988755740245e-5, -1.3299592240934033e-5, -5.992765934740819e-6, -0.004187625548608727, 0.012846266928800682, 0.005605334988727657, 0.0006138501782006211, -0.005091561647229093, -0.002131737664191976, -0.003787865169697781 … -0.49409239904785285, -4.10612582794665, 15.976547800478043, 7.038652554973895, 21.72856241339921, 19.700374188772376, 7.507143626676891, 35.81494990235186, -14.715990463874167, -15.112625320422667]
[1.842220086557274e-5, -7.293778041553316e-6, -3.6152110008688388e-6, -0.01277034693219267, 0.0065467867000238615, 0.003114287240245503, 0.0008361383181451783, -0.005049936683797766, -0.002124160093404663, -0.0037706757685206974 … -0.6435362410331006, -4.3715094026393295, 15.88910713630922, 7.004085774590366, 21.56710603440199, 19.839933375665424, 7.568282459849049, 35.93400450762391, -14.569777566522005, -15.102748368925784]
[1.952839204508782e-5, 3.3263010210562888e-6, 9.096637714818286e-7, -0.014001166445020441, -0.0045892196999754754, -0.001627064087709437, 0.00107777491846707, -0.004994772328192684, -0.002111822753561921, -0.0037507307357908326 … -0.8074249745143877, -4.66246939306433, 15.78857853367493, 6.964146890806444, 21.387916596231804, 19.992397379660687, 7.635144505810718, 36.064279160670075, -14.408318555746828, -15.091481815434232]
[1.2729357726315085e-5, 1.0623236365196787e-5, 4.202101015283348e-6, -0.0069426818543548915, -0.012245476204630935, -0.005078358446560554, 0.0012868863928879673, -0.004938447094777882, -0.002097595903877149, -0.0037323417076362622 … -0.950483792134452, -4.916757677074987, 15.696660509464628, 6.927464240653249, 21.229393779377002, 20.125196232173042, 7.693443426544644, 36.17794175302008, -14.266175307597267, -15.081253774844011]
[2.3324061746960045e-6, 1.1631352467832665e-5, 4.887272036387912e-6, 0.0038863398932770755, -0.013316857333436426, -0.005799741745169013, 0.0014884417067964942, -0.00487636872951648, -0.0020806630838849796, -0.003713557053280832 … -1.0894704349114115, -5.16441513056371, 15.603436167106192, 6.890116706276257, 21.07323229385188, 20.254145481961604, 7.750107986918667, 36.28848999558906, -14.126770149210293, -15.070944168325807]
[-4.998632662128737e-6, 7.26677185913009e-6, 3.1960226138296287e-6, 0.011515696836337373, -0.008759826091143542, -0.004032138694060933, 0.0016516333554088686, -0.004820319640290368, -0.0020645607413344604, -0.0036975338056919617 … -1.202813364052655, -5.367043224743483, 15.524399985687747, 6.85835041235144, 20.94412706504233, 20.359375659227055, 7.796391139947666, 36.378842641685274, -14.011972691956299, -15.062249067450741]
[-8.196748733281896e-6, -2.674831556281519e-6, -9.84435579463287e-7, 0.014801735857826768, 0.0016317472288112597, 0.00034047368464461465, 0.0018635915650100296, -0.00473947676979794, -0.002040320141590123, -0.003675554749983683 … -1.351150488561992, -5.633452191563266, 15.41663792217533, 6.814900347962489, 20.772498812731143, 20.497369357123507, 7.857142511276883, 36.49752574164415, -13.859987082347754, -15.050454409927772]
[-6.074791004805948e-6, -8.148844615174545e-6, -3.380315109041009e-6, 0.012544060567910648, 0.00735226499218046, 0.002846101140970487, 0.001980991703779547, -0.0046906336252165205, -0.0020252125879273476, -0.003662776387797641 … -1.4338719079024254, -5.782759879564289, 15.354302298091426, 6.789699153525852, 20.67535024145355, 20.57453509753865, 7.891143772361628, 36.563995971518985, -13.774265594543378, -15.043661040108478]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")