Error control of the Runge Kutta variable stepsize integrator

[Reneh107]

The way error control is handled by the current implementation of the runge kutta variable stepsize integrator is not correct.

The implementation uses the relative error and the higher order estimate to compute an absolute error tolerance and adds the absolute error tolerance to obtain a total tolerance.

The implementation thus computes a total absolute tolerance based on the relative and absolute tolerance that are specified. The way this is implemented is incorrect. Some examples illustrate this:

state = 100 ;
absTol = 1E-5 , relTol = 1E-12 ; total abs tolerance = 1E-5 + 1E-12 * 100 approx equal to 1E-5 ;
absTol = 1E-12 , relTol = 1E-5 ; total abs tolerance = 1E-12 + 1E-5 * 100 approx equal to 1E-3 ;

This means that in case 2, the absolute tolerance is never satisfied and in case 1 the relative tolerance is never satisfied. Effectively, only the least constraining tolerance is used.

This code is used to calculate the total error tolerance:
// Compute error tolerance based on relative and absolute error tolerances. const StateType errorTolerance_ = ( higherOrderEstimate.array( ).abs( ) * relativeErrorTolerance.array( ) ).matrix( ) + absoluteErrorTolerance;

Furthermore, maybe it is better to put the error controler in a different class, such that different methods can be implement.

[Reneh107]

It would be better to use:

totalAbsoluteErrorTolerance = min( relTol * state , absTol ) ;

In this case the most constraining tolerance is used.

[DominicDirkx]

Hi Rene,

I've thought about it a bit more, and I'm pretty sure that the way it is programmed is the way it is intended to be. However, it should be more clearly documented what is happening (and why!).

The reason that it uses the sum of the two is so that the absolute tolerance is indeed only used where the state apporaches 0 (for instance also what is used in the Matlab numerical integrators). For all other cases, the relative tolerance is used, which means that the integrator looks at the number of digits in the state for which it is 'confident' that it works.

The problem with using an absolute tolerance that should be active everywhere is that this can cause severe issues for problems ranging over many orders of mangitude. For instance, having a 1 mm absolute tolerance may be fine for a state on the order of 100-1000 km, but will may a nightmare on the order of many AUs.

Nevertheless, I understand your point on why the minimum of the two should be used in certain cases. As you say, having a more flexible stepsize control would be a lot better than what we have now.

[Reneh107]

Hi Dominic,

Thank you for the response. I also thought it through and came with the following examples that illustrates what you are saying:

state = 1E4 ;
absTol = 1E-10 , relTol = 1E-12 ; total abs tolerance = 1E-10 + 1E-12 * 1E4 approx equal to 1E-8 ;
So tolerance determined using relative tolerance for large states

state = 1E-5 ;
absTol = 1E-10 , relTol = 1E-12 ; total abs tolerance = 1E-10 + 1E-12 * 1E-5 approx equal to 1E-10 ;
So tolerance determined by absolute tolerance for small states

This means that adding the tolerances ensures that there is a minimum absolute tolerance. Otherwise the absolute tolerance that is used can get very small if the state becomes small.

[DominicDirkx]

That's the point exactly! I've changed the label of this issue to 'Upcoming feature', as this is hopefully something we can make more general in the near future.

[Reneh107]

Since we are now talking about an upcoming feature:

Maybe it is an idea to:

• write a base class rkErrorControler that contains all error control parameters as private members.
• Add functions that compute the step size and validate the step size etc. And use get functions that only retrieve one variable instead of std::pairs (since now everything is saved in the class)
• Add a pointer to the base class to the RK variable numerical integrator class, such that we can remove sending all the control parameters and just use get functions. I prefer it that way, because you don't have long lists with input in your functions (don't say it is the best)
• write a derived class of standardRKErrorControler that contains the implementation that we now have.

In this way we can add several different error controlers and I think it structures the code.

However, I don't think I have the time to implement this.. Just some ideas

[mvandenbroeck]

Hi,

I would like to know if these lines of code describe the default variable step-size control for the RK variable step-size integrator?
This step-size control is given in Eqs. 4.19 to 4.20 in Satellite Orbits of Montenbruck and Gill (2001).

.

Thanks,

Michael

[DominicDirkx]

Hi Michael,

Indeed this code describes the stepsize control that is used.

[mvandenbroeck]

Thanks!

[mvandenbroeck]

It looks like the variable step-size control described in Montenbruck and Gill (2001) (see Eq. (4.21) in the above post) differs from the equation that is implemented in Tudat at this line. In Montenbruck and Gill (2001), p is the higher order, so p+1 becomes the higher order + 1, whereas in the Tudat file, the higherOrder is used in the power.

The following definition is used for p in Montenbruck and Gill (2001):

Maybe I'm missing something, but at least it is good to check.

[DominicDirkx]

Hi Michael,

I'd have to check the precise implementation (any comments, notes etc.?) to see which model is used exactly. Perhaps it is the model from numerical recipes:

https://www2.units.it/ipl/students_area/imm2/files/Numerical_Recipes.pdf

Let's keep this issue open, so I'm reminded to check,

Dominic

[DominicDirkx]

Tracked in #437