# Differences between numerical propagators

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.

• Interesting question and squarely on-topic here in Space SE!
– uhoh
Nov 24, 2021 at 20:10
• Crucial to know would be: are those GMAT reference data computed with exactly the same setup, gravitational potential etc., or just with a setup you think should be the same? Nov 25, 2021 at 16:46
• Yeah, I've set for both the same degree and order for the Geopotential which Is the only perturbation that I'm considering, also the same Orbital elements, epoch and integrator (RK4 with 60 sec of fixed step size). The only difference Is that GMAT considers the EGM-96 model while my propagator considers the EGM-2008 and Is based on IERS Tech Note 36 conventions. Also GMAT results are in the ICRF while mine are in the GCRS. That was the reason why I expected different results but not a such a high error. Nov 25, 2021 at 23:54
• ...well that's already a lot of things that could go wrong. It's almost certainly these details that the discrepancies come from, and the ODE solver has nothing to do with them. Different index / complex etc. conventions for the geopotential coefficients? Something up in the coordinate transformations? And, you're only using quadrupole but GMAT uses the full thing? Nov 26, 2021 at 0:05

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:

• @ChrisR I didn't post the question? I also don't use Simulink, Matlab, or GMAT Nov 24, 2021 at 17:34
• Hi @AlfonsoGonzalez, thank you for your answer. In GMAT (to have consistent results) I've used a Runge Kutta 4 and forced the step size to remain constant setting 60 seconds as minimum and maximum step size. For this reason I expected to get the same results. In any case, I was already thinking about implementing an adaptive step size. I will try to solve it like this. Thanks again for the answer. Nov 24, 2021 at 19:25
• I withdraw my previous comment. I think my eye caught on the "Last modified" author and I thought Alfonso had written the question and answered it. My apologies. Nov 24, 2021 at 22:30
• @Frank we can replace "fuzzy" with "thick" or "wide". Well since there is disagreement we can know that one or both are pretty much just plain wrong. What I recommend is to isolate the integration from anything else. Get you initial state vector and start a simple program in matlab or better yet open-source Python and integrate for one and then a few orbits and compare directly. I'll add an answer about how to do that in a few hours.
– uhoh
Nov 25, 2021 at 10:29
• Adding to @uhoh 's comment, you could also try testing the difference between this GMAT propagator and Matlab's ode45 (as well as your propagator). I would also usually encourage Python but I think OP's company is looking for this in Matlab Nov 25, 2021 at 10:38

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.

• Thanks for your answer! You're right, in fact I've already noted that changing the step size doesn't affect the solution or the deviations with GMAT. The same happens with different integrators. The problem Is that I've already validated the function which computes the derivative of the state vector, because the results are the same of GMATS. Therefore, if the integrator works well and also does the function which computes the accelerations due to the perturbations.. what causes this problem? Also, with a lower altitude and inclination, the error seems lower even If non neglectible. Nov 25, 2021 at 23:42

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.

• Excellent answer! It's nice to see another fan of Adams methods. Do you have any advice on when to choose Adams-Bashforth-Moulton versus Gauss-Radau? I've heard of the latter, but never used them, while ABM I use nearly all the time. Nov 25, 2021 at 13:24
• This answer certainly contains some valid points, but they're phrased in way too much of a blanket fashion. Increasing the solver order because something weird is going on that I don't understand is actually a pretty horrible idea. Yes, higher order can of course give better accuracy, but it doesn't solve fundamental stability problems, and definitely not setup mistakes. RK4 is plenty good enough to explore these things, and the lower order actually makes it more suited to get a feeling for how numerical integrators behave with different ODEs and different constant or adaptive step sizes. Nov 25, 2021 at 16:22
• I agree with @leftaroundabout 1000% I'd say that for such a soft problems like a near circular orbit with only a 1 part per thousand quadrupole term you can get perfectly good results with RK4 or even lower as long as your step size is appropriate (hat tip to DavidHammen). For this kind of problem what a higher order integrator gets you is the same accuracy with a larger step size and so potentially a shorter computation time, which for this will be fractions of a second or less
– uhoh
Nov 25, 2021 at 20:53
• @uhoh yes. Though I wouldn't say “near circular” and “only 1‰ quadrupole” are reasons why a low-order solver is fine. They're rather reasons why there's nothing to be gained from adaptive step size. And this is actually also the scenario where high-order multi-step methods can really shine, allowing large step size and low computational effort to compute millions of orbits into the future. ...only, that is rather pointless if the results are already bogus after only a handful of orbits... To compute long into the future in less regular situations, a lower-order symplectic method may be better. Nov 25, 2021 at 21:27
• @leftaroundabout both can be true, they're not exclusive; but yes certainly adaptive step size doesn't help much in this case. If eccentricity were 0.9 or the magnitude of the quadrupole acceleration were similar to the monopole, and step size was 30 seconds fixed, RK4 would be lousy and an exotic higher order would do better, but yes you are right RK45 would dig in and lower that step size where necessary. I'm working on an answer post now...
– uhoh
Nov 25, 2021 at 22:44