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

R0=x02+y02+z02
r0=3661801.7212+2616815.3682+6388793.0132=7814933.824 м
Вычисляем правые части дифференциальных уравнений
f0=-μ∙x0r03+32∙Jz∙μ∙a02r03∙5∙x0∙z02r02-x0
f0=-3.9860044∙1014∙3661801.7217814933.8243+32·0.001082636∙3.9860044∙1014∙637813627814933.8243∙5∙3661801.721∙6388793.01327814933.8242-3661801.721=-3.050389104
φ0=-μ∙y0r03+32∙Jz∙μ∙a02r03∙5∙y0∙z02r02-y0
φ0=-3.9860044∙1014∙2616815.3687814933.8243+32·0.001082636∙3.9860044∙1014∙637813627814933.8243∙5∙2616815.368∙6388793.01327814933.8242-2616815.368=-2.1798845
ψ0=-μ∙z0r03+32∙Jz∙μ∙a02r03∙5∙z0∙z02r02-3∙z0
ψ0=-3.9860044∙1014∙6388793.0137814933.8243+32·0.001082636∙3.9860044∙1014∙637813627814933.8243∙5∙6388793.013∙6388793.01327814933.8242-3∙6388793.013=-5.33359655
Вычисляем в первом приближении координаты ИСЗ в подшагах
xk(1)=x0+αk∙H∙x0+12∙αk2∙H2∙f0
yk(1)=y0+αk∙H∙y0+12∙αk2∙H2∙φ0
zk(1)=z0+αk∙H∙z0+12∙αk2∙H2∙ψ0
Где
α1 = 0,212340538;
α2 = 0,590533136;
α3 = 0,911412040;
H = 60s
x1(1)=3661801.721+0.212340538∙60∙(1236.484)+12∙0.2123405382∙602∙(-3.050389104)
x1(1)=3677122.191
y1(1)=2616815.368+0.212340538∙60∙6225.923+12∙0.2123405382∙602∙(-2.1798845)
y1(1)=2695959.400
z1(1)=6388793.013+0.212340538∙60∙-3269.095+12∙0.2123405382∙602·(-5.33359655)
z1(1)=6346710.459
Вычисляем в первом приближении r
r1(1)=x1(1)2+y1(1)2+z1(1)2
r1(1)=3677122.1912+(2695959.400)2+6346710.4592=7814739.813 м
x2(1)=3661801.721+0,590533136∙60∙(1236.484)+12∙0,5905331362∙602∙(-3.050389104)
x2(1)=3702264.840
y2(1)=2616815.368+0,590533136∙60∙6225.923+12∙0,5905331362∙602∙(-2.1798845)
y2(1)=2836043.856
z2(1)=6388793.013+0,590533136∙60∙-3269.095+12∙0,5905331362∙602∙(-5.33359655)
z2(1)=6269614.510
r2=x22+y22+z22
r2=3702264.8402+(2836043.856)2+6269614.5102=7813960.315 м
x3(1)=3661801.721+0,911412040∙60∙(1236.484)+12∙0,9114120402∙602∙(-3.050389104)
x3(1)=3721443.661
y3(1)=2616815.368+0,911412040∙60∙6225.923+12∙0,9114120402∙602∙(-2.1798845)
y3(1)=2954018.855
z3(1)=6388793.013+0,911412040∙60∙-3269.095+12∙0,9114120402∙602∙(-5.33359655)
z3(1)=6202048.617
r3=x32+y32+z32
r3=3721443.6612+(2954018.855)2+6202048.6172=7812859.743 м
Вычисляем в первом приближении значения правых частей в подшагах
f11=-μ∙x1(1)r13+32∙Jz∙μ∙aе2r15∙5∙x1(1)∙z1(1)2r12-x1(1)
f11=-3.063524495
φ11=-μ∙y1(1)r13+32∙Jz∙μ∙aе2r15∙5∙y1(1)∙z1(1)2r12-y1(1)
φ11=-2.246087357
ψ11=-μ∙z1(1)r13+32∙Jz∙μ∙aе2r13∙5∙z1(1)∙z1(1)2r12-3∙z1(1)
ψ11=-5.287641244
f21=-μ∙x2(1)r23+32∙Jz∙μ∙aе2r25∙5∙x2(1)∙z2(1)2r22-x2(1)
f21=-3.08565763
φ21=-μ∙y2(1)r23+32∙Jz∙μ∙aе2r25∙5∙y2(1)∙z2(1)2r22-y2(1)
φ21=-2.363704581
ψ21=-μ∙z2(1)r23+32∙Jz∙μ∙aе2r25∙5∙z2(1)∙z2(1)2r22-3∙z2(1)
ψ21=-5.225418679
f31=-μ∙x3(1)r33+32∙Jz∙μ∙aе2r35∙5∙x3(1)∙z3(1)2r32-x3(1)
f31=-3.103180387
φ31=-μ∙y3(1)r33+32∙Jz∙μ∙aе2r35∙5∙y3(1)∙z3(1)2r32-y3(1)
φ31=-2.463251955
ψ31=-μ∙z3(1)r33+32∙Jz∙μ∙aе2r35∙5∙z3(1)∙z3(1)2r32-3∙z3(1)
ψ31=-5.171669217
Вычисляем в первом приближении p, q, r
p1(1)q1(1)r1(1)p2(1)q2(1)r2(1)p2(1)q3(1)r3(1)=-0.063066671-0.3120871120.2477608270.0056942510.001515124-0.169052314-0.000054967-0.000245680.101101201
Вычисляем во втором приближении координаты ИСЗ в подшагах
xk(2)=x0+αk∙H∙x0+12∙αk2∙H2∙f0+αk2∙H2∙16∙αk∙p11+112∙αk∙p21+120∙αk∙p31
yk(1)=y0+αk∙H∙y0+12∙αk2∙H2∙φ0+αk2∙H2∙16∙αk∙q11+112∙αk∙q21+120∙αk∙q31
zk(1)=z0+αk∙H∙z0+12∙αk2∙H2∙ψ0+αk2∙H2∙16∙αk∙r11+112∙αk∙r21+120∙αk∙r31
x1(2)=3677306.663
y1(2)=2695959.146
z1(2)=6346710.632
r1=x12+y12+z12
r1=3677306.6632+(2695959.146)2+6346710.6322=7814826.668 м
x2(2)=3703680.149
y2(2)=2836038.387
z2(2)=6269618.236
r2(2)=x22+y22+z22
r2(2)=3703680.1492+(2836038.387)2+6269618.2362=7814631.994 м
x3(2)=3724791.767
y3(2)=2953998.748
z3(2)=6202062.313
r3(2)=x32+y32+z32
r3(2)=3724791.7672+(2953998.748)2+6202062.3132=7814458.347м
Вычисляем во втором приближении значения правых частей в подшагах
f12=-μ∙x1(2)r13+32∙Jz∙μ∙aе2r15∙5∙x1(2)∙z1(2)2r12-x1(2)
f12=-3.063576448
φ12=-μ∙y1(2r13+32∙Jz∙μ∙aе2r15∙5∙y1(2)∙z1(2)2r12-y1(2)
φ12=-2.246012558
ψ12=-μ∙z1(2)r13+32∙Jz∙μ∙aе2r13∙5∙z1(2)∙z1(2)2r12-3∙z1(1)
ψ12=-5.287465799
f22=-μ∙x2(1)r23+32∙Jz∙μ∙aе2r25∙5∙x2(1)∙z2(1)2r22-x2(1)
f22=-3.08604445
φ21=-μ∙y2(1)r23+32∙Jz∙μ∙aе2r25∙5∙y2(1)∙z2(1)2r22-y2(1)
φ22=-2.36309297
ψ22=-μ∙z2(1)r23+32∙Jz∙μ∙aе2r25∙5∙z2(1)∙z2(1)2r22-3∙z2(1)
ψ22=-5.224079776
f32=-μ∙x3(1)r33+32∙Jz∙μ∙aе2r35∙5∙x3(1)∙z3(1)2r32-x3(1)
f32=-3.104073733
φ32=-μ∙y3(1)r33+32∙Jz∙μ∙aе2r35∙5∙y3(1)∙z3(1)2r32-y3(1)
φ32=-2.461729539
ψ32=-μ∙z3(1)r33+32∙Jz∙μ∙aе2r35∙5∙z3(1)∙z3(1)2r32-3∙z3(1)
ψ32=-5.168519455
Вычисляем во втором приближении p, q, r
p1(2)q1(2)r1(2)p2(2)q2(2r2(2)p2(2)q3(2)r3(2)=-0.063068047-0.3120879210.2477589710.0045266940.003125955-0.1651202490.000047749-0.0000012140.100949565
Вычисляем окончательные значения скорости и координаты
х=x0+H∙f0+∙H12∙p1+13∙p2+14∙p3
х=1049.043 м/с
y=y0+H∙φ0+∙H12∙q1+13∙q2+14∙q3
y=6092.851 м/с
z=z0+H∙Ψ0+∙H12∙r1+13∙r2+14∙r3
z=-3587.595141 м/с
х=х0+H∙х0+12∙H2∙f0+ H16∙p1+112∙p2+120∙p3
х= 3730413.190 м
y=y0+H∙y0+12∙H2∙φ0+ H16∙q1+112∙q2+120∙q3
y=2986420.888 м
z=z0+H∙z0+12∙H2∙Ψ0+ H16∙r1+112∙r2+120∙r3
z= 6183065.038 м