⚙ Numerical accuracies in orbit propagation via RK4 for the two-body case.

[sammmlow]

🏹 Goal: Figure out the extent of numerical integration errors using various tested numerical propagators.

We attempt to investigate a numerical propagator that should be efficient (i.e. Runge-Kutta variants) but not overly complex to the point that it makes orbit propagation slow on MATLAB, as we might need to propagate orbits for days, and for N number of satellites in an arbitrary configuration.

Current position propagation errors for a satellite with initial conditions:

a1  = 7000000;     % Semi-major axis (m)
e1  = 0.001;       % Eccentricity (unitless)
i1  = 10.00;       % Inclination (degrees)
w1  = 90.00;       % Arg of Periapsis (degrees)
R1  = 45.00;       % Right Ascension (degrees)
M1  = -90.00;      % Mean Anomaly (degrees)
Cd1 = 2.2;         % Drag coefficient
Ar1 = 0.374;       % Drag area (m^2)
Ms1 = 17.90;       % Spacecraft mass (kg)

Propagator is a modified RK4 using the 3/8 rule, and the errors are compared with the analytical two-body solution (i.e. rotational solutions of an object about an elliptical orbit). Pretty significant errors.

[sammmlow]

Related case to using RK4 or ODE45 for the J2-inclusive numerical propagator (one full day duration).

RK4 propagator using J2 and two-body produces about 1-km error as compared the STK HPOP's J2 (STK uses an RKF78 propagator).

Important point to note: STK's two-body and J2 orbit propagator uses the analytical solution (the J2 or J4 propagator uses the analytical solution and applies the mean element change rate to the osculating elements at each time step, rather than a force model integration). Thus, using J2-mode and HPOP-with-J2-only will produce slightly different results.

Related topic on stack exchange: https://engineering.stackexchange.com/questions/14822/j2-orbit-propagation-problem-in-matlab-solver-instability

[ANDREWNGT]

Hi, would it be right to say that the main contributor to the errors we see above is because the STK's J2 propagator uses secular rates of orbital precession, as identified by AliRD in the stack exchange thread?

[sammmlow]

Hi, would it be right to say that the main contributor to the errors we see above is because the STK's J2 propagator uses secular rates of orbital precession, as identified by AliRD in the stack exchange thread?

Yup, so what it does is that it directly applies the rate of change of the RAAN, and any positions and velocities are converted from the osculating orbit elements at that point in time. On the other hand, in Matthew's code, the J2 effect is taken as an additional force term after the GMm-Over-R-Squared force, and then numerically integrated. I suppose just from a modelling perspective we should definitely see some differences.