Численное интегрирование системы дифференциальных уравнений возмущённого движения ИСЗ

R0=x02+y02+z02
r0=3902901.7882+4677876.6612+4890253.5192=7812154.114 м
Вычисляем правые части дифференциальных уравнений
f0=-μ∙x0r03+32∙Jz∙μ∙a02r03∙5∙x0∙z02r02-x0
f0=-3.9860044∙1014∙3902901.7887812154.1143+32·0.001082636∙3.9860044∙1014∙637813627812154.1143∙5∙3902901.788∙4890253.51927812154.1142-3902901.788=-3.259581132
φ0=-μ∙y0r03+32∙Jz∙μ∙a02r03∙5∙y0∙z02r02-y0
φ0=-3.9860044∙1014∙3902901.7887812154.1143+32·0.001082636∙3.9860044∙1014∙637813627812154.1143∙5∙4677876.661∙4890253.51927812154.1142-4677876.661=-3.906815834
ψ0=-μ∙z0r03+32∙Jz∙μ∙a02r03∙5∙z0∙z02r02-3∙z0
ψ0=-3.9860044∙1014∙3902901.7887812154.1143+32·0.001082636∙3.9860044∙1014∙637813627812154.1143∙5∙4890253.519∙4890253.51927812154.1142-4890253.519=-4.093037609
Вычисляем в первом приближении координаты ИСЗ в подшагах
xk1=x0+αk∙H∙x0+12∙αk2∙H2∙f0
yk1=y0+αk∙H∙y0+12∙αk2∙H2∙φ0
zk1=z0+αk∙H∙z0+12∙αk2∙H2∙ψ0
Где
α1 = 0.212340538
α2 = 0.590533136
α3 = 0.91141204
H = 60s
x11=3902901.788+0.212340538∙60∙90.381+12∙0.21234053882∙602∙-3.259581132
x11=3903721.093 м
y11=4677876.661+0.212340538∙60∙5120.61+12∙0.21234053882∙602∙-3.906815834
y11=4742798.371 м
z11=4890253.519+0.212340538∙60∙-4981.564+12∙0.21234053882∙602∙-4.093037609
z11=4826454.052 м
Вычисляем в первом приближении r11
r11=x112+y112+z112
r11=3903721.0932+4742798.3712+4826454.0522=7812031.329 м
x21=3902901.788+0.590533136∙60∙90.381+12∙0.59053313682∙602∙-3.259581132
x21=3903534.914 м
y21=4677876.661+0.590533136∙60∙5120.61+12∙0.59053313682∙602∙-3.906815834
y21=4856857.695 м
z21=4890253.519+0.590533136∙60∙-4981.564+12∙0.59053313682∙602∙-4.093037609
z21=4711177.550 м
r21=x212+y212+z212
r21=3903534.9142+4856857.6952+4711177.5502=7811648.059 м
x31=3902901.788+0.91141204∙60∙90.381+12∙0.9114120482∙602∙-3.259581132
x31=3901724.299 м
y31=4677876.661+0.91141204∙60∙5120.61+12∙0.9114120482∙602∙-3.906815834
y31=4952054.289 м
z31=4890253.519+0.91141204∙60∙-4981.564+12∙0.9114120482∙602∙-4.093037609
z31=4611718.126 м
r31=x312+y312+z312
r21=3901724.2992+4952054.2892+4611718.1262=7811161.135 м
Вычисляем в первом приближении значения правых частей в подшагах
f11=-μ∙x11r13+32∙Jz∙μ∙aе2r15∙5∙x11∙z112r12-x11
f11=
-3.260598244
φ11=-μ∙y11r13+32∙Jz∙μ∙aе2r15∙5∙y11∙z112r12-y11
φ11=
-3.961440809
ψ11=-μ∙z11r13+32∙Jz∙μ∙aе2r13∙5∙z11∙z112r12-3∙z11
ψ11=
-4.031314542
f21=-μ∙x21r23+32∙Jz∙μ∙aе2r25∙5∙x21∙z212r22-x21
f21=
-3.261240048
φ21=-μ∙y21r23+32∙Jz∙μ∙aе2r25∙5∙y21∙z212r22-y21
φ21=
-4.057701332
ψ21=-μ∙z21r23+32∙Jz∙μ∙aе2r25∙5∙z21∙z212r22-3∙z21
ψ21=
-3.935991668
f31=-μ∙x31r33+32∙Jz∙μ∙aе2r35∙5∙x31∙z312r32-x31
f31=
-3.260604317
φ31=-μ∙y31r33+32∙Jz∙μ∙aе2r35∙5∙y31∙z312r32-y31
φ31=
-4.13834714
ψ31=-μ∙z31r33+32∙Jz∙μ∙aе2r35∙5∙z31∙z312r32-3∙z31
ψ31=
-3.853934025
Вычисляем в первом приближении p, q, r
p11q11r11p21q21r21p21q31r31= -0.005898859 -0.258236082 0.314294761
0.005216456 0.004643857 -0.127680634
0.000026341 -3.87719E-05 0.077550498
Вычисляем во втором приближении координаты ИСЗ в подшагах
xk2=x0+αk∙H∙x0+12∙αk2∙H2∙f0+αk2∙H2∙16∙αk∙p11+112∙αk∙p21+120∙αk∙p31
yk1=y0+αk∙H∙y0+12∙αk2∙H2∙φ0+αk2∙H2∙16∙αk∙q11+112∙αk∙q21+120∙αk∙q31
zk1=z0+αk∙H∙z0+12∙αk2∙H2∙ψ0+αk2∙H2∙16∙αk∙r11+112∙αk∙r21+120∙αk∙r31
x12=
3903788.717 м
y12=
4742796.901 м
z12=
4826455.624 м
r12=x122+y122+z122
r12=3903788.7172+4742796.9012+4826455.6242=7812065.2 м
x22=
3904057.680 м
y22=
4856826.073 м
z22=
4711211.371 м
r22=x222+y222+z222
r22=3904057.682+4856826.0732+4711211.3712=7811910.039 м
x32=
3902969.000 м
y32=
4951938.035 м
z32=
4611842.463 м
r32=x322+y322+z322
r32=39029692+4951938.0352+4611842.4632=7811782.654 м
Вычисляем во втором приближении значения правых частей в подшагах
f12=-μ∙x12r13+32∙Jz∙μ∙aе2r15∙5∙x12∙z122r12-x12
f12=
-3.260612397
φ12=-μ∙y12r13+32∙Jz∙μ∙aе2r15∙5∙y12∙z122r12-y12
φ12=
-3.961388153
ψ12=-μ∙z12r13+32∙Jz∙μ∙aе2r13∙5∙z12∙z122r12-3∙z11
ψ12=
-4.03126352
f22=-μ∙x22r23+32∙Jz∙μ∙aе2r25∙5∙x22∙z222r22-x22
f22=
-3.261349191
φ22=-μ∙y22r23+32∙Jz∙μ∙aе2r25∙5∙y22∙z222r22-y22
φ22=
-4.057267356
ψ22=-μ∙z22r23+32∙Jz∙μ∙aе2r25∙5∙z22∙z222r22-3∙z22
ψ22=
-3.935624587
f32=-μ∙x32r33+32∙Jz∙μ∙aе2r35∙5∙x32∙z322r32-x32
f32=
-3.260867114
φ32=-μ∙y32r33+32∙Jz∙μ∙aе2r35∙5∙y32∙z322r32-y32
φ32=
-4.137263679
ψ32=-μ∙z32r33+32∙Jz∙μ∙aе2r35∙5∙z32∙z322r32-3∙z32
ψ32=
-3.853119361
Вычисляем во втором приближении p, q, r
p12q12r12p22q22r22p22q32r32= -0.005900984 -0.2582387 0.314292265
0.004915656 0.005785429 -0.126488727
0.000011820 0.000142945 0.077321799
Вычисляем в третьем приближении координаты ИСЗ в подшагах
xk3=x0+αk∙H∙x0+12∙αk2∙H2∙f0+αk2∙H2∙16∙αk∙p11+112∙αk∙p21+120∙αk∙p31
yk3=y0+αk∙H∙y0+12∙αk2∙H2∙φ0+αk2∙H2∙16∙αk∙q11+112∙αk∙q21+120∙αk∙q31
zk3=z0+αk∙H∙z0+12∙αk2∙H2∙ψ0+αk2∙H2∙16∙αk∙r11+112∙αk∙r21+120∙αk∙r31
x13=
3903788.716 м
y13=
4742796.905 м
z13=
4826455.628 м
r13=x132+y132+z132
r13=3903788.7162+4742796.9052+4826455.6282=7812065.203 м
x23=
3904057.660 м
y23=
4856826.150 м
z23=
4711211.436 м
r23=x232+y232+z232
r23=3904057.662+4856826.152+4711211.4362=7811910.116 м
x33=
3902968.929 м
y33=
4951938.318 м
z33=
4611842.702 м
r33=x332+y332+z332
r33=3902968.9292+4951938.3182+4611842.7022=7811782.938 м
Вычисляем в третьем приближении значения правых частей в подшагах
f13=-μ∙x13r13+32∙Jz∙μ∙aе2r15∙5∙x13∙z132r12-x13
f13=
-3.260612392
φ13=-μ∙y13r13+32∙Jz∙μ∙aе2r15∙5∙y13∙z132r12-y13
φ13=
-3.961388151
ψ13=-μ∙z13r13+32∙Jz∙μ∙aе2r13∙5∙z13∙z132r12-3∙z13
ψ13=
-4.031263517
f23=-μ∙x23r23+32∙Jz∙μ∙aе2r25∙5∙x23∙z232r22-x23
f23=
-3.261349078
φ23=-μ∙y23r23+32∙Jz∙μ∙aе2r25∙5∙y23∙z232r22-y23
φ23=
-4.0572673
ψ23=-μ∙z23r23+32∙Jz∙μ∙aе2r25∙5∙z23∙z232r22-3∙z23
ψ23=
-3.935624524
f33=-μ∙x33r33+32∙Jz∙μ∙aе2r35∙5∙x33∙z332r32-x33
f33=
-3.260866698
φ33=-μ∙y33r33+32∙Jz∙μ∙aе2r35∙5∙y33∙z332r32-y33
φ33=
-4.137263464
ψ33=-μ∙z33r33+32∙Jz∙μ∙aе2r35∙5∙z33∙z332r32-3∙z33
ψ33=
Вычисляем в третьем приближении p, q, r
p13q13r13p23q23r23p23q33r33= -0.005900984 -0.258238696 0.314292261
0.004915655 0.005785396 -0.126488698
0.000012371 0.000143261 0.077322064
Вычисляем окончательные значения скорости и координаты
х=x0+H∙f0+∙H12∙p1+13∙p2+14∙p3
х=
-105.82 м/с
y=y0+H∙φ0+∙H12∙q1+13∙q2+14∙q3
y=
4884.57 м/с
z=z0+H∙Ψ0+∙H12∙r1+13∙r2+14∙r3
z=
-5225.98 м/с
х=х0+H∙х0+12∙H2∙f0+ H16∙p1+112∙p2+120∙p3
x=
3902432.962 м
y=y0+H∙y0+12∙H2∙φ0+ H16∙q1+112∙q2+120∙q3
y=
4978062.910 м
z=z0+H∙z0+12∙H2∙Ψ0+ H16∙r1+112∙r2+120∙r3
z=
4584006.180 м