Differences between numerical propagators

Question

I am a trainee who is working on a numeric orbital propagator developed in the company. I can't show you the code but I can tell you that the propagator was developed to work in Simulink. My job was to extract the propagator in order to use it as a Matlab script. To evaluate the functioning of my script, I am comparing the results obtained with those of another numerical propagator (GMAT). The initial orbital elements are from a SSO Repeat Ground Track orbit. At the moment I am considering the disturbances due only to J2 and I am using a fourth order runge kutta type integrator. The analyzes with the two propagators over a period of 52 days with a fixed step size of 60 sec have the different trends shown in the figure (the left shows the orbital elements while the right shows the state vector).

I cannot explain this difference, particularly in the state vector from which the orbital elements are derived. At first I thought it was due to the reference systems, as the GMAT data is in the ICRF while my propagator data is in the GCRS. So I used the same conversion functions on the GMAT data and comparing (these new results) with the original ones there were no errors, except those due to the numeric. So I thought it was the geopotential function that was wrong. But even in this case, passing the same state vector from GMAT, the accelerations calculated by the function were the same as those provided by the software. So I think it may be the integrator's fault, as there is the periodic drift of the error shown in the next figure.

The Runge-Kutta, however, also comparing with what I have found on the web seems to be set up well. Can you think of a reason that explains these discrepancies?

I add more plots to show you the trend of orbital elements and state vector in 24 hours. Over time, the differences tend to increase. In particular, my propagator seems accelerated compared to that of GMAT.

Answer 1

You should check to see what propagator GMAT is using. A constant step size RK4 method is going to drift from the true solution of the ODE over time, so if GMAT is using an adaptive step size method that at some points goes from the order of milliseconds to tens of seconds, its likely going to be closer to the solution than RK4.

In essence, this is what RK4 constant step size is doing: Its taking multiple evaluations of the derivative between now and the next time step in order to get a better estimate of how the function is changing (as opposed to say a 1st order euler's method). It then takes a weighted mean of those estimates in order to determine what derivative it will use to take this step. So as it can be seen from the plot, at some points it underestimates the area under the curve and at others it overestimates, which will build error over time.

Here is a zoomed out version (and 500 second time step) to show the effect.

Since you're using 60 seconds the error will grow less quickly but nevertheless will drift over time.

An adaptive step size solver will change its step size based on how fast the derivative is changing (the second derivative):

So if the derivative is closer to constant, it means that the solution is close to linear, so the solver can take bigger time steps. But if the derivative is changing, the solution is nonlinear, the solver needs to take smaller time steps.

The most common of these methods is RK4-5, which compares the solution from the RK4 and RK5 methods, and based on how different they are will change the time step accordingly.

So overall there may not be anything wrong with your model / RK4 implementation, it may be depending on what propagator GMAT is using that you're comparing your solution to

Also if you're curious these screenshots come from a video I made on ODE solvers:

Answer 2

The previous answers contain some good general information, but TBH I don't find it very helpful to the question.

The most fundamental sanity check with a numerical solver is: how does the result change when you vary the step size? Keep it as simple as possible. Simulate over a single day, and try different step sizes. Like, try a step of 6 seconds, 30 seconds, 200 seconds, 600 seconds. How do they differ?

The 600 seconds probably will probably behave very weird, but the others – I don't know, you should try. If the h = 6 seconds solution is essentially the same as your h=60 s one, in particular if it has consistently the same deviations from the GMAT – then you evidently have a problem that is more fundamental and doesn't have to do with the integrator at all. This could be lots of different things.

Only if you've established that you've definitely set up everything as it should be, should you focus on getting the best out of your numerical solver.

Answer 3

Your numerical integrator, RK4, is inappropriate for the task. The order is far too low for the accuracy you want. It's also fixed step size. You haven't done your research on this. Error propagation rates for fixed-step RK integrators is about as bad as it gets. Research high-order predictor-corrector integrators...order 8 with (at least) 64-bit floating point operations (even if they are fixed step).

You might look at embedded RK integrators that use a lower-order RK integrator as a predictor within a higher order RK integrator. Look at RK78 as a replacement for your RK4. In the past, that integrator has been a workhorse for orbital mechanicians. This one will give you much better data agreements and good accuracy over longer time intervals.

Adams-Moulton, or better yet, Adams-Bashforth-Moulton, are predictor-correctors that can be coded easily to any order you want. Their error propagation rates are less than any RK method of the same order. Gauss-Radau integrators are also very accurate and well-suited to orbital mechanics integrations. Don't code them yourself. Download.

There are many variable-step predictor-corrector methods which will give the best accuracy and smallest error propagation rates over longer time intervals. With these methods your numerical integrations will produce the best accuracy you can expect...but they have to be high order. Order 8 is good.

Before jumping into your project with just any integrator mentioned here or elsewhere, test your integrator selections against the solution to the 2-body problem expressed in Cartesian coordinates. That's the worst way to write the 2-body equations, and that form of the equations will stress your integrators best. Since the 2-body problem solution is closed form in any coordinate system, you will see how selected integrators fair with respect to error propagation and accuracy over short and long time intervals.

Your project isn't trivial by any means, and you will need an integrator that is up to the task. You will get a good idea of the integrator to use by numerical integrations of the 2-body problem. Make sure that your initial conditions will give you a 2-body orbit of decently high eccentricity: e >= 0.3. Don't go crazy with very high e (0.8). Fixed-step integrators will die quickly with high-eccentricity orbits (e = 0.8). Their rate of death is inversely proportional to the integrator's order.

Good luck.