Моделирование внешней Солнечной системы
Данные
В качестве единицы измерения выбрана масса относительно Солнца, то есть Солнце имеет массу . Мы принимаем для учета массы планет земной группы. Расстояния выражаются в астрономических единицах, время — в земных сутках, а постоянная силы тяжести, таким образом, имеет значение .
| Планета | Масса | Начальное положение | Начальная скорость |
|---|---|---|---|
Юпитер |
|
<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> |
Сатурн |
|
<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> |
Уран |
|
<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> |
Нептун |
|
<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> |
Плутон |
|
<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> |
Данные взяты из книги Geometric Numerical Integration (Геометрическое численное интегрирование) Э. Хайрера, К. Любича и Г. Ваннера.
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)Гамильтониан задачи N тел имеет следующий вид:
В данном случае мы хотим решить задачу движения пяти внешних планет, а именно Юпитера, Сатурна, Урана, Нептуна и Плутона относительно Солнца.
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(7.692353667846155e-9 / sqrt((-(u(t))[1, 1] + (u(t))[1, 6])^2 + (-(u(t))[2, 1] + (u(t))[2, 6])^2 + (-(u(t))[3, 1] + (u(t))[3, 6])^2) + 2.1967979473153846e-12 / sqrt((-(u(t))[1, 3] + (u(t))[1, 6])^2 + (-(u(t))[2, 3] + (u(t))[2, 6])^2 + (-(u(t))[3, 3] + (u(t))[3, 6])^2) + 0.00028558544003356794 / sqrt((-(u(t))[1, 1] + (u(t))[1, 3])^2 + (-(u(t))[2, 1] + (u(t))[2, 3])^2 + (-(u(t))[3, 1] + (u(t))[3, 3])^2) + 4.372757780489954e-5 / sqrt((-(u(t))[1, 1] + (u(t))[1, 4])^2 + (-(u(t))[2, 1] + (u(t))[2, 4])^2 + (-(u(t))[3, 1] + (u(t))[3, 4])^2) + 0.0009547918106276825 / sqrt((-(u(t))[1, 1] + (u(t))[1, 2])^2 + (-(u(t))[2, 1] + (u(t))[2, 2])^2 + (-(u(t))[3, 1] + (u(t))[3, 2])^2) + 1.4786358763131087e-8 / sqrt((-(u(t))[1, 3] + (u(t))[1, 5])^2 + (-(u(t))[2, 3] + (u(t))[2, 5])^2 + (-(u(t))[3, 3] + (u(t))[3, 5])^2) + 7.344508492638461e-12 / sqrt((-(u(t))[1, 2] + (u(t))[1, 6])^2 + (-(u(t))[2, 2] + (u(t))[2, 6])^2 + (-(u(t))[3, 2] + (u(t))[3, 6])^2) + 1.2487810273779817e-8 / sqrt((-(u(t))[1, 3] + (u(t))[1, 4])^2 + (-(u(t))[2, 3] + (u(t))[2, 4])^2 + (-(u(t))[3, 3] + (u(t))[3, 4])^2) + 2.7267137995329905e-7 / sqrt((-(u(t))[1, 2] + (u(t))[1, 3])^2 + (-(u(t))[2, 2] + (u(t))[2, 3])^2 + (-(u(t))[3, 2] + (u(t))[3, 3])^2) + 5.177622330021739e-5 / sqrt((-(u(t))[1, 1] + (u(t))[1, 5])^2 + (-(u(t))[2, 1] + (u(t))[2, 5])^2 + (-(u(t))[3, 1] + (u(t))[3, 5])^2) + 2.2640217694220477e-9 / sqrt((-(u(t))[1, 4] + (u(t))[1, 5])^2 + (-(u(t))[2, 4] + (u(t))[2, 5])^2 + (-(u(t))[3, 4] + (u(t))[3, 5])^2) + 4.175023411794291e-8 / sqrt((-(u(t))[1, 2] + (u(t))[1, 4])^2 + (-(u(t))[2, 2] + (u(t))[2, 4])^2 + (-(u(t))[3, 2] + (u(t))[3, 4])^2) + 3.363639727276923e-13 / sqrt((-(u(t))[1, 4] + (u(t))[1, 6])^2 + (-(u(t))[2, 4] + (u(t))[2, 6])^2 + (-(u(t))[3, 4] + (u(t))[3, 6])^2) + 4.9434923063238094e-8 / sqrt((-(u(t))[1, 2] + (u(t))[1, 5])^2 + (-(u(t))[2, 2] + (u(t))[2, 5])^2 + (-(u(t))[3, 2] + (u(t))[3, 5])^2) + 3.982762603453846e-13 / sqrt((-(u(t))[1, 5] + (u(t))[1, 6])^2 + (-(u(t))[2, 5] + (u(t))[2, 6])^2 + (-(u(t))[3, 5] + (u(t))[3, 6])^2))Гамильтонова система
NBodyProblem, так сказать «за кадром», формирует задачу ОДУ второго порядка. Известно, что гамильтонова система имеет следующую форму:
Для системы N тел ее можно упростить так:
Таким образом, определяется массами. Достаточно определить , и это делается автоматически путем получения градиента . Поэтому нам нужно лишь передать гармоническую функцию, а остальное будет сделано за нас.
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.90807728347
199550.3018430433
199620.52140424005
199697.88811276815
199765.85051860654
199832.37128055916
199887.05311094294
199959.3168813624
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.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.562666124366707e-7, -1.5291383285282026e-6, -6.369559468830294e-7, 0.006730897179308314, -0.00258065612604747, -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.2700600085295772e-6, -3.3960857323731747e-6, -1.4300961196408633e-6, 0.007365006196212581, -0.0006664595881567761, -0.0004651305514942438, 0.00044221410614463444, 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.0075134582329184e-6, -5.6313841026625595e-6, -2.3985120866280563e-6, 0.007324487954155873, 0.0016548597238463265, 0.0005308765916475192, -0.0003215351780469663, 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.55466763421876e-7, -8.08224049085999e-6, -3.4873974695025767e-6, 0.006174492455065487, 0.004236423620030293, 0.0016654664602168783, -0.0012242830831770996, 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.49018116256438e-6, -1.0242904678740652e-5, -4.497045822643825e-6, 0.003230726642949305, 0.006575790349577204, 0.0027399381613809654, -0.002325513928515103, 0.004752067449942419, 0.0020626776193825294, 0.0030372322298507144 … 0.980818022484472, 12.472577563460696, -14.085404632753704, -6.345502797745267, 15.000995668497582, -24.071220490411154, -10.2261812806247, -11.927011166605329, -27.14391250145533, -4.876236637628356]
[7.821334963992697e-6, -1.0696453179221612e-5, -4.805135026272915e-6, -0.000967826079200862, 0.0072014961995934944, 0.0031104154878608863, -0.003416481639860906, 0.00420282263218379, 0.0018827832856910382, 0.0028597245490350033 … 1.6967093574323593, 13.538407580658337, -13.289128533479984, -6.0118582932866556, 15.969038388790834, -23.51787256963099, -10.023818114603046, -10.849436012073225, -27.609190660804593, -5.345590705748484]
⋮
[-1.3756970627695715e-6, -1.1856984706630836e-5, -5.093597970681905e-6, 0.008108473909196294, 0.011342873757654758, 0.004663311500494419, 0.0003563597623574362, -0.005129121778814678, -0.0021361091055617973, -0.003806475284902511 … -0.3225260944445467, -3.8008830539358556, 16.07219051306757, 7.076252034834613, 21.91200257299775, 19.539243441348937, 7.436628618226829, 35.67770829208628, -14.882978228964651, -15.12352282822378]
[1.0290958757390644e-5, -1.3299598628325022e-5, -5.992767770954497e-6, -0.004187602530511477, 0.012846272613164754, 0.005605336866295312, 0.000613878310091792, -0.005091556922404325, -0.0021317369299466322, -0.0037878630536486288 … -0.49411123588117895, -4.1061594197485975, 15.976536986338717, 7.038648290737175, 21.72854209570609, 19.700391883051672, 7.507151374468558, 35.81496498537105, -14.715972020184203, -15.112624094207108]
[1.8422184495472906e-5, -7.293791089728327e-6, -3.6152160171201716e-6, -0.012770338150859593, 0.006546799045321322, 0.003114292316343306, 0.0008361663069990827, -0.005049930896357213, -0.002124158913622823, -0.0037706735275822843 … -0.6435551391128408, -4.371543069589128, 15.88909578781262, 7.004081277440546, 21.567085435117082, 19.839951048153193, 7.568290205733158, 35.934019594726436, -14.56975895702035, -15.102747092071935]
[1.95283884974606e-5, 3.3262854666349616e-6, 9.096573732368673e-7, -0.014001170983887508, -0.004589205063449963, -0.0016270577070957868, 0.0010778025191502305, -0.00499476542486096, -0.0021118210955948465, -0.003750728370773701 … -0.8074437824348492, -4.66250291676283, 15.788566667059692, 6.964142164792934, 21.387895818799922, 19.99241491252846, 7.635152198986754, 36.06429415474363, -14.408299883118902, -15.091480490811962]
[1.2729360798610836e-5, 1.0623228723530922e-5, 4.202097848860816e-6, -0.0069426932390552415, -0.0122454701471734, -0.005078355575343604, 0.0012869136813447037, -0.004938439213339738, -0.0020975938279811884, -0.003732339227278664 … -0.9505025380195785, -4.916791151439979, 15.696648157362487, 6.927459301379585, 21.22937279330445, 20.12521368525648, 7.693451092261363, 36.177956702978136, -14.266156532549475, -15.081252404894752]
[2.3324000737610397e-6, 1.1631350391530173e-5, 4.887271479261856e-6, 0.003886338171611383, -0.013316857396985259, -0.005799741730122335, 0.001488468791020256, -0.004876359861545542, -0.0020806605912545074, -0.0037135544461730507 … -1.089489188270397, -5.164448719814184, 15.60342327276436, 6.890111531194441, 21.073210995053113, 20.254162944228852, 7.750115664180507, 36.28850497817888, -14.126751177376857, -15.070942746647455]
[-4.998646920089926e-6, 7.266761699357857e-6, 3.1960187934719827e-6, 0.01151570369556528, -0.008759817931977743, -0.004032135363684724, 0.0016516602980104218, -0.004820309946896218, -0.0020645579012944304, -0.003697531087510426 … -1.2028321411625005, -5.3670769715147655, 15.52438661347566, 6.858345030242764, 20.944105463129624, 20.359393163214396, 7.796398841831731, 36.37885768134841, -14.011953518193495, -15.06224759979627]
[-8.196754138526984e-6, -2.6748476798018883e-6, -9.84442167632637e-7, 0.014801733571921667, 0.001631761336106387, 0.0003404797849933768, 0.0018636180623509, -0.004739466076900978, -0.0020403168690100244, -0.003675551905719826 … -1.3511691152659502, -5.633485863969452, 15.416624022391302, 6.8148947338155175, 20.772476984089963, 20.497386771923665, 7.857150182303481, 36.49754073387676, -13.859967796702009, -15.050452892986796]
[-6.075558781521503e-6, -8.14775531349734e-6, -3.3798298026366263e-6, 0.012544864734528162, 0.007351124111132237, 0.0028455928534587244, 0.0019809916133655347, -0.004690633662960368, -0.002025212599956213, -0.0036627763874126835 … -1.4338718457448323, -5.782759883865783, 15.354302296415245, 6.789699152827176, 20.675350240679148, 20.574535098153937, 7.891143772626507, 36.56399597173869, -13.774265594102516, -15.043661040101771]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")