Engee documentation

Simulating the Outer Solar System

Data

The chosen units are masses relative to the sun, meaning the sun has mass . We have taken to take account of the inner planets. Distances are in astronomical units, times in earth days, and the gravitational constant is thus .

planet mass initial position initial velocity

Jupiter

<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

<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

<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

<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

<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

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)
-0.000295912208286(0.0009547918106276825 / sqrt(((u(t))[1, 2] - (u(t))[1, 1])^2 + ((u(t))[2, 2] - (u(t))[2, 1])^2 + ((u(t))[3, 2] - (u(t))[3, 1])^2) + 0.00028558544003356794 / sqrt(((u(t))[1, 3] - (u(t))[1, 1])^2 + ((u(t))[2, 3] - (u(t))[2, 1])^2 + ((u(t))[3, 3] - (u(t))[3, 1])^2) + 2.7267137995329905e-7 / sqrt(((u(t))[1, 3] - (u(t))[1, 2])^2 + ((u(t))[2, 3] - (u(t))[2, 2])^2 + ((u(t))[3, 3] - (u(t))[3, 2])^2) + 4.372757780489954e-5 / sqrt(((u(t))[1, 4] - (u(t))[1, 1])^2 + ((u(t))[2, 4] - (u(t))[2, 1])^2 + ((u(t))[3, 4] - (u(t))[3, 1])^2) + 4.175023411794291e-8 / sqrt(((u(t))[1, 4] - (u(t))[1, 2])^2 + ((u(t))[2, 4] - (u(t))[2, 2])^2 + ((u(t))[3, 4] - (u(t))[3, 2])^2) + 1.2487810273779817e-8 / sqrt(((u(t))[1, 4] - (u(t))[1, 3])^2 + ((u(t))[2, 4] - (u(t))[2, 3])^2 + ((u(t))[3, 4] - (u(t))[3, 3])^2) + 5.177622330021739e-5 / sqrt(((u(t))[1, 5] - (u(t))[1, 1])^2 + ((u(t))[2, 5] - (u(t))[2, 1])^2 + ((u(t))[3, 5] - (u(t))[3, 1])^2) + 7.692353667846155e-9 / sqrt(((u(t))[1, 6] - (u(t))[1, 1])^2 + ((u(t))[2, 6] - (u(t))[2, 1])^2 + ((u(t))[3, 6] - (u(t))[3, 1])^2) + 4.9434923063238094e-8 / sqrt(((u(t))[1, 5] - (u(t))[1, 2])^2 + ((u(t))[2, 5] - (u(t))[2, 2])^2 + ((u(t))[3, 5] - (u(t))[3, 2])^2) + 7.344508492638461e-12 / sqrt(((u(t))[1, 6] - (u(t))[1, 2])^2 + ((u(t))[2, 6] - (u(t))[2, 2])^2 + ((u(t))[3, 6] - (u(t))[3, 2])^2) + 1.4786358763131087e-8 / sqrt(((u(t))[1, 5] - (u(t))[1, 3])^2 + ((u(t))[2, 5] - (u(t))[2, 3])^2 + ((u(t))[3, 5] - (u(t))[3, 3])^2) + 2.1967979473153846e-12 / sqrt(((u(t))[1, 6] - (u(t))[1, 3])^2 + ((u(t))[2, 6] - (u(t))[2, 3])^2 + ((u(t))[3, 6] - (u(t))[3, 3])^2) + 2.2640217694220477e-9 / sqrt(((u(t))[1, 5] - (u(t))[1, 4])^2 + ((u(t))[2, 5] - (u(t))[2, 4])^2 + ((u(t))[3, 5] - (u(t))[3, 4])^2) + 3.363639727276923e-13 / sqrt(((u(t))[1, 6] - (u(t))[1, 4])^2 + ((u(t))[2, 6] - (u(t))[2, 4])^2 + ((u(t))[3, 6] - (u(t))[3, 4])^2) + 3.982762603453846e-13 / sqrt(((u(t))[1, 6] - (u(t))[1, 5])^2 + ((u(t))[2, 6] - (u(t))[2, 5])^2 + ((u(t))[3, 6] - (u(t))[3, 5])^2))

Hamiltonian System

NBodyProblem constructs a second order ODE problem under the hood. We know that a Hamiltonian system has the form of

For an N-body system, we can symplify this as:

Thus, is defined by the masses. We only need to define , and this is done internally by taking the gradient of . 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.90711457172
 199550.30091018166
 199620.5204162125
 199697.88712605814
 199765.84950142307
 199832.37014898597
 199887.05186152307
 199959.31559004786
 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.1048711440018423e-9, -1.278530308573887e-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.5706941518695234e-8, -3.423355432416458e-8, -1.4100666088910381e-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.4909025937027347e-7, -3.5335373517313927e-7, -1.4594000158085006e-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.562666124366708e-7, -1.5291383285282024e-6, -6.369559468830295e-7, 0.006730897179308314, -0.002580656126047471, -0.0012701900733253107, 0.0010966413711339348, 0.004998334414871358, 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.270060008529577e-6, -3.396085732373174e-6, -1.430096119640863e-6, 0.007365006196212581, -0.000666459588156777, -0.0004651305514942436, 0.0004422141061446342, 0.005105982369438331, 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.007513458232918e-6, -5.631384102662555e-6, -2.3985120866280567e-6, 0.0073244879541558726, 0.0016548597238463245, 0.0005308765916475187, -0.00032153517804696673, 0.005139085078150228, 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.554667634218723e-7, -8.082240490859979e-6, -3.487397469502574e-6, 0.00617449245506549, 0.004236423620030282, 0.0016654664602168748, -0.0012242830831771, 0.005052449172558457, 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.490181162564379e-6, -1.0242904678740642e-5, -4.497045822643826e-6, 0.0032307266429492964, 0.006575790349577194, 0.002739938161380965, -0.002325513928515104, 0.004752067449942418, 0.002062677619382529, 0.0030372322298507144  …  0.9808180224844723, 12.472577563460696, -14.085404632753704, -6.345502797745267, 15.000995668497582, -24.071220490411154, -10.2261812806247, -11.92701116660533, -27.14391250145533, -4.876236637628356]
 [7.821334963992685e-6, -1.0696453179221606e-5, -4.805135026272916e-6, -0.0009678260792008517, 0.00720149619959349, 0.0031104154878608867, -0.0034164816398609093, 0.004202822632183787, 0.0018827832856910373, 0.0028597245490350033  …  1.6967093574323595, 13.538407580658337, -13.289128533479984, -6.0118582932866556, 15.969038388790834, -23.51787256963099, -10.023818114603046, -10.849436012073227, -27.609190660804593, -5.345590705748484]
 ⋮
 [-1.37567814544341e-6, -1.1856996002307405e-5, -5.093603287107087e-6, 0.008108455021373986, 0.01134288565897322, 0.0046633170571464405, 0.0003563566609783263, -0.005129122163104611, -0.0021361091300507312, -0.003806475500967649  …  -0.3225240378071363, -3.8008793893153974, 16.07219162999071, 7.076252472523893, 21.912004761019208, 19.539241502893763, 7.436627770377224, 35.677706642356696, -14.882980226344653, -15.123522956113641]
 [1.029098651113033e-5, -1.3299591622612805e-5, -5.992765460085725e-6, -0.004187630711624032, 0.012846265381976127, 0.005605334450865502, 0.0006138753396056656, -0.005091557421254793, -0.0021317370074751436, -0.0037878632762399795  …  -0.4941092470831806, -4.106155886115259, 15.976538123937592, 7.038648739314681, 21.72854423317387, 19.700390021589417, 7.5071505593893315, 35.81496339858804, -14.715973960522717, -15.11262422321371]
 [1.842219652666053e-5, -7.293774520492941e-6, -3.615209229248645e-6, -0.012770349822796888, 0.00654678183821452, 0.0031142852268028222, 0.0008361631980679568, -0.005049931539148188, -0.002124159044666142, -0.0037706737755620347  …  -0.6435530402043296, -4.371539343978699, 15.889097043670548, 7.004081775107142, 21.56708771484374, 19.839949092352104, 7.568289348500879, 35.93401792501758, -14.569761016558825, -15.102747233387811]
 [1.952838260395013e-5, 3.3263057878745786e-6, 9.096662026701873e-7, -0.014001163896529643, -0.004589226163635941, -0.001627066919100134, 0.001077799460581168, -0.004994766189758649, -0.0021118212793072356, -0.0037507286318641353  …  -0.8074416985010336, -4.6624992158002625, 15.788567977139534, 6.964142686545758, 21.387898112789877, 19.992412976779896, 7.635151349606442, 36.06429249926517, -14.408301944743643, -15.091480637066821]
 [1.2729344544157419e-5, 1.0623237566824988e-5, 4.202102010825884e-6, -0.006942675285028866, -0.01224547918760745, -0.00507835988367803, 0.001286910574743401, -0.004938440110460052, -0.002097594064282013, -0.003732339508590371  …  -0.9505004042017027, -4.916787354885344, 15.69664955832317, 6.9274598615847305, 21.229375173679212, 20.12521170563226, 7.693450222772768, 36.17795500724116, -14.266158662157904, -15.081252560289625]
 [2.3323988349254586e-6, 1.163135033594507e-5, 4.887271462247707e-6, 0.0038863404877014147, -0.013316857057563481, -0.005799741640844737, 0.0014884653899861471, -0.00487636097497631, -0.002080660904233445, -0.003713554772336299  …  -1.0894868336681027, -5.164444517575353, 15.603424885959482, 6.890112178640866, 21.073213659901413, 20.254160759413455, 7.750114703628272, 36.28850310357726, -14.126753551106589, -15.070942924531085]
 [-4.99863577906425e-6, 7.266773623302826e-6, 3.196023606547223e-6, 0.011515693136271371, -0.008759830079494829, -0.004032140312448565, 0.0016516565960137486, -0.004820311278617864, -0.002064558291486837, -0.0036975314596195507  …  -1.2028295614813833, -5.367072351666058, 15.524388444127235, 6.858345767049813, 20.94410842062927, 20.35939076677482, 7.796397787379347, 36.37885562225705, -14.011956143285085, -15.062247800739291]
 [-8.196755983857613e-6, -2.67482757930912e-6, -9.844335377497345e-7, 0.014801736628967317, 0.001631740677065171, 0.00034047085962588664, 0.001863614313615383, -0.004739467589480904, -0.0020403173319562123, -0.0036755523066370503  …  -1.3511664803978896, -5.6334811175959, 15.416625981698106, 6.814895525180893, 20.77248006125926, 20.49738431698723, 7.857149100929919, 36.49753862040674, -13.859970515415805, -15.050453106836915]
 [-6.075461069455182e-6, -8.147894411939697e-6, -3.3798917660126695e-6, 0.012544762394824703, 0.007351269797483901, 0.002845657751503947, 0.001980991612933968, -0.0046906336630678095, -0.0020252125999756796, -0.0036627763874170966  …  -1.4338718455693498, -5.782759883786203, 15.35430229645479, 6.789699152843823, 20.6753502407047, 20.57453509814778, 7.891143772623499, 36.56399597173851, -13.774265594120587, -15.043661040104551]
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")