J2 orbit propagation problem in MATLAB (solver instability)

Question

The problem is that when I try to solve the following equations of J2 perturbed gravity:

function xdot = Gravity( ~ , x )

Re  = 6378.137;
J2  = 1.08262668e-3;
mu  = 398600.4418;

r   = norm(x(1:3));

% xdot = [x(4:6)
%         -mu/r^3*x(1:3)];

xdot    = [ x(4);
            x(5);
            x(6);
            -mu*x(1)/r^3*(1+1.5*J2*(Re/r)^2*(1-5*(x(3)/r)^2));
            -mu*x(2)/r^3*(1+1.5*J2*(Re/r)^2*(1-5*(x(3)/r)^2));
            -mu*x(3)/r^3*(1+1.5*J2*(Re/r)^2*(3-5*(x(3)/r)^2))];
end

using any of the solvers (ode45, ode113, etc.), the propagated position and velocity has a considerable error compared to a trusted commercial software (STK), and the error keeps growing. I've tried various settings of tolerances, initial steps, and max steps, but the problem persists!

It might be surprising, but when I wrote a fixed step 4th order Runge-Kutta code as below, it didn't get much better.

function [t,x] = RK4( fn,tspan,x0 )
% An implementation of fixed step 4th order Runge-Kutta method

f = str2func(fn);
t = tspan;
h = tspan(2)-tspan(1);
x(:,1) = x0;

for i = 1:size(t,2)-1
    k1 = f(t(i),x(:,i));
    k2 = f(t(i)+h/2,x(:,i)+k1*h/2);
    k3 = f(t(i)+h/2,x(:,i)+k2*h/2);
    k4 = f(t(i)+h,x(:,i)+k3*h);
    x(:,i+1) = x(:,i) + (k1 + 2*k2 + 2*k3 + k4)*h/6;
end

end

any idea of what the problem might be? I also have to mention that to a great surprise, I did get a result using RK4, but the same code with the same conditions misbehaved later!!

A sample result:

Answer 1

As far as I can see, you are trying to integrate a nonlinear ODE of a conservative system. It is very likely that your system exhibits chaotic behavior. RK4 is a fourth order numerical scheme. These numerical integrators tend to change the total energy of the system, which contradicts the conservativeness of the system but explains the differences between your solution and the solution of the program you are talking about (btw I am not sure if the program does integrate the equations correctly, but let's assume it does. You could look up the user manual or contact the support in order to find out if the software is suited for your problem). In order to get better numerical integration, you have different possibilities.

1. Use higher order approximation scheme (e.g. 8th order numerical scheme): This system will still show a drift in the total energy, but the drift will be less rapidly than for RK4. So your simulation will follow the "real" trajectory for a longer period. But as you can imagine at some point your trajectory will not follow the "real" trajectory anymore.

2. Use custom integrators like symplectic integrators (for Hamiltonian systems: keeps Hamiltonian constant) or variational integrators (for Lagrangian systems: keeps the Lagrangian constant). The hard part about these integrators is that you need to create them problem dependent, but they are usually much better than just increasing the order of the numerical scheme.

Answer 2

Solved it, thanks to a colleague. The problem is, I used "J2 Propagator" in STK, which in fact uses secular rates of orbital precession. Comparing that result with the real integration that MATLAB does, led to this result. When I used HPOP propagator with degree 2 and order 0 (which is the same as just using J2 perturbation), the results got much better. HPOP does indeed integrate the differential equations. There was a growing error at a rate of 500 meters per day, and that is due to the fact that HPOP considers earth's axis precession, nutation, and polar motion, which those equations above do not contain these phenomena.

I also need to mention, to get a correct solution of J2 perturbed orbit, use a relative tolerance of 1e-8 or finer.

Answer 3

https://www.mathworks.com/matlabcentral/fileexchange/55167-high-precision-orbit-propagator?s_tid=srchtitle

https://www.researchgate.net/publication/340793133_High_Precision_Orbit_Propagator_C_code

The motion of a near-Earth satellite is affected by various forces. One of these forces is the Earth's central gravitation and the others are known as perturbations. These perturbations are classified into gravitational and non-gravitational forces. In this case, the equation of motion can be written as: r ̈=-(GM/r^3)*r+γ_p

γ_p is the vector of additional accelerations induced by the disturbing forces. γ_p=r ̈_E+r ̈_S+r ̈_M+r ̈_P+r ̈_e+r ̈_o+r ̈_D+r ̈_SP+r ̈_A+r ̈_emp r ̈_E = Accelerations due to the non-spherically and inhomogeneous mass distribution within Earth (central body) r ̈_S, r ̈_M, r ̈_P = Accelerations due to other celestial bodies (Sun, Moon, and planets) r ̈_e, r ̈_o = Accelerations due to Earth and oceanic tides r ̈_D = Accelerations due to atmospheric drag r ̈_SP, r ̈_A = Accelerations due to direct and Earth-reflected solar radiation pressure r ̈_emp = Accelerations due to unmodeled forces Here I have used the following integrator and force model to simulate the satellite's perturbed motion: Integrator: Variable-order Radau IIA integrator with step-size control Force Model:

• gravity field of the Earth (GGM03S model)
• gravity of the solar system planets (positions of the planets are computed by JPLDE436)
• drag effect using Jacchia-Bowman 2008, NRLMSISE-00, MSIS-86, Jacchia 70 or modified Harris-Priester atmospheric density model (in Accel.m you can uncomment your favorite model)
• solar radiation pressure using geometrical or cylindrical shadow model
• solid Earth tides (IERS Conventions 2010)
• ocean tides
• general relativity
• ECEF2ECI and ECI2ECEF transformations using IAU 2006 Resolution

The Simulation starts by running the test_HPOP.m. In the InitialState.txt set initial values for your favorite satellite; lines 2-7 are the state vector of satellite/spacecraft in the International Terrestrial Reference Frame (ITRF). Lines 8 to 12 are the satellite parameters related to the forces from atmospheric drag and solar radiation pressure. Lines 8-10 are in units m^2 and kg. Line 11: Cr is the radiation pressure coefficient (Cr = 1 + reflectivity of satellite). Line 12: Cd is the atmospheric drag coefficient of the satellite. In the test_HPOP.m you can consider different perturbations by setting them 1 as follows:

AuxParam.n = 70; % maximum degree of central body's gravitation field

AuxParam.m = 70; % maximum order of central body's gravitation field

AuxParam.sun = 1; % perturbation of the Sun

AuxParam.moon = 1; % perturbation of Moon

AuxParam.planets = 1; % perturbations of planets

AuxParam.sRad = 1; % solar radiation pressure

AuxParam.drag = 1; % atmospheric drag

AuxParam.SolidEarthTides = 1; % solid Earth tides

AuxParam.OceanTides = 1; % ocean tides

AuxParam.Relativity = 1; % general relativity

References:

Montenbruck O., Gill E.; Satellite Orbits: Models, Methods and Applications; Springer Verlag, Heidelberg; Corrected 3rd Printing (2005).

Montenbruck O., Pfleger T.; Astronomy on the Personal Computer; Springer Verlag, Heidelberg; 4th edition (2000).

Seeber G.; Satellite Geodesy; Walter de Gruyter, Berlin, New York; 2nd completely revised and extended edition (2003).

Vallado D. A; Fundamentals of Astrodynamics and Applications; McGraw-Hill, New York; 3rd edition(2007).

NIMA. 2000. Department of Defense World Geodetic System 1984. NIMA-TR 8350.2, 3rd ed, Amendment 1. Washington, DC: Headquarters, National Imagery, and Mapping Agency.

http://sol.spacenvironment.net/jb2008/ https://ssd.jpl.nasa.gov/planets/eph_export.html https://celestrak.org/SpaceData/

Answer 4

Try to use ECEF data in perturbation calculation then execute integration part with ECI data after making the proper conversion from ECEF to ECI reference frame.

Ref: D. A. Vallado, Fundamentals of Astrodynamics and Applications, 4th Editon, Microcosm Press, 2013, pp.591-595