Step size computation in AdaptiveStepsizeFieldIntegrator

For DerivativeStructure instances (and perhaps for the other classes that implement the Derivative interface) the step size computation in AdaptiveStepsizeFieldIntegrator.initializeStep() and AdaptiveStepsizeFieldIntegrator.filterStep() should not consider the derivatives of the instance but only its real value returned by getReal().

A simple (but not efficient) solution would be to replace return h with return h.getField().getZero().add(h.getReal()) (line 316) and return filteredH with return filteredH.getField().getZero().add(filteredH.getReal()) (line 348).

This change do not influence the real part of neither the step nor the output state, but does improve a lot the accuracy of the higher-order terms, especially in the intermediate steps. In fact, the coefficients of the highest-order derivatives as computed now can be orders of magnitude bigger in absolute value than the expected ones.

Hi @afossa

I have looked at the issue and am puzzled. I find the rationale for filterStep very sound so will most certainly include it. Concerning initializeStep, I would like to have more arguments for it. Fields are not used only for DerivativeStructure (or Gradient, UnivariateDerivative1 and UnivariateDerivative2) but also for other types like Dfp, so using a primitive double instead of a Field instance always looks risky.

Could you provide a unit test showing a case where derivatives are wrong (and not in singular cases like square roots near origin or similar corner cases)?

Hi @maisonobe

I perfectly agree with you that for Dfp objects and similar the remark do not apply, so that why I was suggesting that perhaps it should be applied to classes that implement the Derivative interface.

However, I dug a bit more into the code and for EmbeddedRungeKuttaFieldIntegrator integrators the problem is due to the error estimation routines, since an estimation in double using only the constant coefficient is more coherent with the fact that usually the time is a DerivativeStructure object with only the real part different from zero and also the if condition in integrate() that checks if the step should be rejected calls the getReal() method.

So the modification should be limited to initializeStep() (computation of ratio, ratioDot and ratioDotDot in double) and estimateError() which should return a constant T object. I haven't look at them yet, but I assume that similar modifications will be needed for Adams integrators too.

I'll attach the script I was testing when I discovered the error and the results I obtained with both Hipparchus, a modified version to truncate the estimated error and a third-party reference solution. They are the expected states (array of 4 DerivativeStructure object) at the second to last step (step 58) since the error is much more visible here. I'll attach also the time at each step for reference (you still need to propagate until the end since the last step is computed differently and that helps a lot to reduce the error).

Here is the error w.r.t. the reference I obtain truncating/non truncating the estimated error during integration grouped by component (rows) and derivative order (cols)

Truncated error

                             0                              1                              2                              3                              4                              5                              6 
        3.2003288907844762e-12         4.3848653874523080e-11         4.9865986448578700e-11         3.9119596451087090e-11         4.6877498410813345e-11         2.4347704408178572e-11         3.4548266157619080e-11 
        7.1677899726729780e-11         6.8767019856252890e-11         4.9570667015608194e-11         4.4327368164953640e-11         1.7067357943600925e-11         2.8495284257190612e-11         1.9217690806760945e-11 
        5.7284732513096515e-11         5.6854894707122905e-11         8.6436247048737870e-11         8.3968450248494970e-11         6.4252554665689130e-11         1.2195921876567706e-10         8.5858055907939160e-11 
        1.6100454303114020e-12         6.3324789856267220e-11         6.0740981688844900e-11         7.1678597085567120e-11         1.0325487687090629e-10         4.9174383315309100e-11         8.8443736709353790e-11

Polynomial error

                             0                              1                              2                              3                              4                              5                              6 
        3.2003288907844762e-12         1.1994260208779417e-03         3.6084408905483840e-01         2.9466679697181486e-01         6.2871584658579090e-02         9.4522026172128140e-02         5.0265904269655055e-02 
        7.1677899726729780e-11         9.4956500533512830e-01         4.4142033681465714e-01         1.5801316588742179e-01         2.4104978704899568e-01         4.5431531056133023e-01         6.5553490130984000e-01 
        5.7284732513096515e-11         7.7531518858084290e-01         3.6043318430899050e-01         2.5629676808285210e-01         3.8524207979496680e-01         3.8668124640913165e-01         5.2902208918578900e-01 
        1.6100454303114020e-12         1.4689891084217511e-03         4.4192860910630320e-01         3.6508354293464523e-01         1.4778146609717357e-01         2.4735758458354146e-01         1.4040344821248948e-01

hope this can help test_stepsize.zip

I will give a try at this in the next few days. I am a little overloaded right now.

Thank you.

I also found today that this problem affects issue #587 in Orekit and posted a comment there. You could also use that test case to reproduce the error.

I have investigated the problem further and think I have found the reason for the wrong high order derivatives. In order to do that, I simplified your test case (which was nice, with very complete output that helped me a lot, thanks for that). The simplified problem is a simple ellipse using the eccentric anomaly as the free parameter so all derivatives can be computed analytically without relying on solving Kepler equation. The code is here:

    private static class Ellipse<T extends RealFieldElement<T>> implements FieldOrdinaryDifferentialEquation<T> {

        private final T a;
        private final T b;
        private final T omega;

        Ellipse(final T a, final T b, final T omega) {
            this.a     = a;
            this.b     = b;
            this.omega = omega;
        }

        @Override
        public int getDimension() {
            return 4;
        }

        @Override
        public T[] computeDerivatives(T t, T[] y) {
            final T[] yDot = MathArrays.buildArray(t.getField(), getDimension());
            yDot[0] = y[2];
            yDot[1] = y[3];
            yDot[2] = y[0].multiply(omega.multiply(omega).negate());
            yDot[3] = y[1].multiply(omega.multiply(omega).negate());
            return yDot;
        }

        public T[] computeTheoreticalState(T t) {
            final T[] y = MathArrays.buildArray(t.getField(), getDimension());
            final FieldSinCos<T> sc = FastMath.sinCos(t.multiply(omega));
            y[0] = a.multiply(sc.cos());
            y[1] = b.multiply(sc.sin());
            y[2] = a.multiply(omega).multiply(sc.sin()).negate();
            y[3] = b.multiply(omega).multiply(sc.cos());
            return y;
        }

    }

With this problem, I added a few columns to the state output from your test case. In addition to displaying the integrated value ∂ⁿsᵢ/∂pᵤ∂pᵥ (tₖ) for state parameter sᵢ at time tₖ corresponding to step k with respect to derivation parameters pᵤ, pᵥ… I also print the corresponding theoretical value (let's call it thᵢ (tₖ)) and the theoretical value thᵢ (trunc(tₖ)). These three values are computed as follows in the step handler:

            DerivativeStructure   tK       = interpolator.getCurrentState().getTime();
            DerivativeStructure[] state    = interpolator.getCurrentState().getPrimaryState();
            DerivativeStructure   tKtrunc  = tK.getFactory().constant(tK.getReal());
            DerivativeStructure[] thK      = ellipse.computeTheoreticalState(tK);
            DerivativeStructure[] tkKtrunc = ellipse.computeTheoreticalState(tKtrunc);

At initial time and also after one period, all values are equal, but as you noticed, there are problems with high order derivatives at intermediate steps. Here are the results I get for step 3, computing derivatives up to order 3 with respect to three parameters: the length of the X ellipse axis (set to 2.0), the length of the Y ellipse axis (set to 1.0), the angular pulsation (set to 0.5).

Step number:   3
Time:             0.6063921404363682

     I          INTEGRATED              THEORETICAL(tk)      THEORETICAL(trunc(tk))  ORDER EXPONENTS
     1    1.9087742161796480e+00   1.9087742161796320e+00   1.9087742161796320e+00    0    0  0  0
     2    9.5828540375391780e-01   9.5828540375391250e-01   9.5438710808981600e-01    1    1  0  0
     3    1.0901351604216203e-03   1.0901351604217770e-03   0.0000000000000000e+00    2    2  0  0
     4   -1.8441046723209948e-04  -1.8441046723216214e-04   0.0000000000000000e+00    3    3  0  0
     5    5.7807906805598410e-03   5.7807906805627510e-03   0.0000000000000000e+00    1    0  1  0
     6    1.9173795499964910e-04   1.9173795499973250e-04   0.0000000000000000e+00    2    1  1  0
     7   -4.3950687130086990e-04  -4.3950687130078997e-04   0.0000000000000000e+00    3    2  1  0
     8   -4.9647587729954130e-04  -4.9647587730015730e-04   0.0000000000000000e+00    2    0  2  0
     9    1.0681579079090601e-03   1.0681579079093086e-03   0.0000000000000000e+00    3    1  2  0
    10   -7.8637783078656550e-04  -7.8637783078677220e-04   0.0000000000000000e+00    3    0  3  0
    11   -9.5260360373726900e-02  -9.5260360373733170e-02  -3.6210344603499467e-01    1    0  0  1
    12   -3.9941289723244490e-02  -3.9941289723247710e-02  -1.8105172301749733e-01    2    1  0  1
    13    2.7778792890761914e-03   2.7778792890765947e-03   0.0000000000000000e+00    3    2  0  1
    14    4.3700010087671440e-03   4.3700010087715510e-03   0.0000000000000000e+00    2    0  1  1
    15   -3.0971844118563257e-03  -3.0971844118559100e-03   0.0000000000000000e+00    3    1  1  1
    16    1.6716814762223766e-03   1.6716814762193048e-03   0.0000000000000000e+00    3    0  2  1
    17    5.6144481749214810e-04   5.6144481749330340e-04  -3.5093904636427214e-01    2    0  0  2
    18    4.6704222639700746e-03   4.6704222639724850e-03  -1.7546952318213607e-01    3    1  0  2
    19   -1.0659985629716543e-02  -1.0659985629720833e-02   0.0000000000000000e+00    3    0  1  2
    20   -5.8893424329735470e-03  -5.8893424329820750e-03   2.2191595869848854e-02    3    0  0  3
------------------------------------------------

     I          INTEGRATED              THEORETICAL(tk)      THEORETICAL(trunc(tk))  ORDER EXPONENTS
     1    2.9857201461613586e-01   2.9857201461616900e-01   2.9857201461616900e-01    0    0  0  0
     2   -6.2304619040000200e-03  -6.2304619040049886e-03   0.0000000000000000e+00    1    1  0  0
     3    1.3015499588523640e-03   1.3015499588546647e-03   0.0000000000000000e+00    2    2  0  0
     4   -3.6176109797284730e-04  -3.6176109797388936e-04   0.0000000000000000e+00    3    3  0  0
     5    2.8933284980830130e-01   2.8933284980832960e-01   2.9857201461616900e-01    1    0  1  0
     6   -2.1289929895928275e-03  -2.1289929895957870e-03   0.0000000000000000e+00    2    1  1  0
     7   -1.3745261875320810e-05  -1.3745261874219695e-05   0.0000000000000000e+00    3    2  1  0
     8   -8.6026121394751790e-03  -8.6026121394790260e-03   0.0000000000000000e+00    2    0  2  0
     9    2.1524198614097344e-03   2.1524198614110597e-03   0.0000000000000000e+00    3    1  2  0
    10    1.9154841622539536e-03   1.9154841622552954e-03   0.0000000000000000e+00    3    0  3  0
    11    1.5225013631897788e-01   1.5225013631896740e-01   5.7873284127945900e-01    1    0  0  1
    12   -8.8007646453718840e-03  -8.8007646453674080e-03   0.0000000000000000e+00    2    1  0  1
    13    7.6368024432561130e-04   7.6368024432305930e-04   0.0000000000000000e+00    3    2  0  1
    14    1.5043817091788933e-01   1.5043817091787250e-01   5.7873284127945900e-01    2    0  1  1
    15   -3.0268321036552200e-03  -3.0268321036494590e-03   0.0000000000000000e+00    3    1  1  1
    16   -4.9251541458121190e-03  -4.9251541458118400e-03   0.0000000000000000e+00    3    0  2  1
    17   -4.3514751044075414e-02  -4.3514751044068210e-02  -5.4894170925136340e-02    2    0  0  2
    18   -2.8246833790752292e-03  -2.8246833790789740e-03   0.0000000000000000e+00    3    1  0  2
    19   -2.6554122543562125e-02  -2.6554122543543873e-02  -5.4894170925136340e-02    3    0  1  2
    20    3.1646816919467975e-02   3.1646816919481346e-02  -3.5467779914588140e-02    3    0  0  3
------------------------------------------------

     I          INTEGRATED              THEORETICAL(tk)      THEORETICAL(trunc(tk))  ORDER EXPONENTS
     1   -2.9857201461613586e-01  -2.9857201461616900e-01  -2.9857201461616900e-01    0    0  0  0
     2   -1.4305554540406798e-01  -1.4305554540407950e-01  -1.4928600730808450e-01    1    1  0  0
     3    1.8136809931476480e-03   1.8136809931478296e-03  -0.0000000000000000e+00    2    2  0  0
     4   -2.8901388145333504e-04  -2.8901388145344290e-04  -0.0000000000000000e+00    3    3  0  0
     5    9.2391648078345170e-03   9.2391648078394060e-03  -0.0000000000000000e+00    1    0  1  0
     6    5.1811348951006680e-04   5.1811348951050130e-04  -0.0000000000000000e+00    2    1  1  0
     7   -7.3543923647591360e-04  -7.3543923647571650e-04  -0.0000000000000000e+00    3    2  1  0
     8   -6.3655266835934140e-04  -6.3655266836037930e-04  -0.0000000000000000e+00    2    0  2  0
     9    1.6307727188177904e-03   1.6307727188179526e-03  -0.0000000000000000e+00    3    1  2  0
    10   -1.2789314938946121e-03  -1.2789314938949161e-03  -0.0000000000000000e+00    3    0  3  0
    11   -7.4939416555124900e-01  -7.4939416555130540e-01  -1.1758768705117970e+00    1    0  0  1
    12   -3.5343539432225260e-01  -3.5343539432227533e-01  -5.8793843525589850e-01    2    1  0  1
    13    7.2640640646556230e-03   7.2640640646563035e-03  -0.0000000000000000e+00    3    2  0  1
    14    2.0290295016758010e-02   2.0290295016773690e-02  -0.0000000000000000e+00    2    0  1  1
    15   -3.8317228621520590e-03  -3.8317228621495072e-03  -0.0000000000000000e+00    3    1  1  1
    16    1.8400834080044464e-03   1.8400834079962038e-03  -0.0000000000000000e+00    3    0  2  1
    17   -2.6098552159388044e-01  -2.6098552159386657e-01  -1.1025715116337818e+00    2    0  0  2
    18   -1.1006654812712208e-01  -1.1006654812711950e-01  -5.5128575581689090e-01    3    1  0  2
    19   -1.3336697698335792e-02  -1.3336697698334583e-02  -0.0000000000000000e+00    3    0  1  2
    20    5.5382685168683520e-02   5.5382685168655080e-02   1.4525612176486083e-01    3    0  0  3
------------------------------------------------

     I          INTEGRATED              THEORETICAL(tk)      THEORETICAL(trunc(tk))  ORDER EXPONENTS
     1    4.7719355404491200e-01   4.7719355404490800e-01   4.7719355404490800e-01    0    0  0  0
     2    9.7457391602344250e-04   9.7457391602412390e-04   0.0000000000000000e+00    1    1  0  0
     3   -2.1475316790631540e-04  -2.1475316790661768e-04   0.0000000000000000e+00    2    2  0  0
     4    6.1273967145132620e-05   6.1273967145268300e-05   0.0000000000000000e+00    3    3  0  0
     5    4.7863875171505194e-01   4.7863875171504866e-01   4.7719355404490800e-01    1    0  1  0
     6    2.9990956970337500e-04   2.9990956970371314e-04   0.0000000000000000e+00    2    1  1  0
     7    1.2702287428500056e-05   1.2702287428390213e-05   0.0000000000000000e+00    3    2  1  0
     8    1.3210787008150725e-03   1.3210787008156484e-03   0.0000000000000000e+00    2    0  2  0
     9   -3.4556538468035840e-04  -3.4556538468056390e-04   0.0000000000000000e+00    3    1  2  0
    10   -3.2071342702152630e-04  -3.2071342702173240e-04   0.0000000000000000e+00    3    0  3  0
    11    9.3057201799639230e-01   9.3057201799638270e-01   8.6386124658106730e-01    1    0  0  1
    12    3.8713704479515438e-03   3.8713704479529670e-03   0.0000000000000000e+00    2    1  0  1
    13   -6.9614782149591090e-04  -6.9614782149644630e-04   0.0000000000000000e+00    3    2  0  1
    14    9.3455491358886400e-01   9.3455491358885690e-01   8.6386124658106730e-01    2    0  1  1
    15    1.2014955262514386e-03   1.2014955262517240e-03   0.0000000000000000e+00    3    1  1  1
    16    4.1525780228775390e-03   4.1525780228790110e-03   0.0000000000000000e+00    3    0  2  1
    17   -4.7489818982490580e-02  -4.7489818982493257e-02  -2.6878648460856536e-01    2    0  0  2
    18    4.9418701956157670e-03   4.9418701956158970e-03   0.0000000000000000e+00    3    1  0  2
    19   -4.7969814885536144e-02  -4.7969814885537690e-02  -2.6878648460856536e-01    3    0  1  2
    20   -1.1916131994970103e-03  -1.1916131994988670e-03  -1.6992162421467386e-01    3    0  0  3
------------------------------------------------

There are three things to notice in this output: 1) the first two columns are equal 2) the third column is different with the first two columns and shows regular patterns of zeros 3) the first two columns do not exhibit the pattern of zeros

The fact the first two columns are equal means that the integration algorithm works perfectly and computes derivatives correctly, including high order ones that match the expected theoretical values that are in the simple implementation of the method computeTheoreticalState.

The fact the third column exhibits a pattern of zeros is expected and corresponds to mathematical equations: the X coordinate and its first time-derivative depend only on parameters a and omega and not on b, the Y coordinate and its first time-derivative depend only on parameters b and omega and not on a.

So the problem is due to fact 3: the first two columns do not exhibit the expected mathematical pattern of zeroes despite columns 2 and 3 are computed by exactly the same method computeTheoreticalState.

The reason for this discrepancy is that when we compute column 2, we use the time as extracted from the computed step, and as we are using an adaptive step handler, this time tk is not an independant variable anymore, it is a complex computation and depends on a, b and omega. Just removing this dependency by replacing time tk by tK.getFactory().constant(tK.getReal()) before calling computeTheoreticalState is sufficient to get a perfectly correct result.

So what can we do? Obviously, derivatives are computed correctly under the assumption that tₖ is itself the result of a computation and they are computed incorrectly under the assumption that tₖ should always remain an independent variable.

As I understand your patch, you already did the same analysis and decided to select the second option: forcing tₖ to be an independent variable by dropping the derivatives in the computations that make tₖ evolve, i.e. in time step computation.

I have been playing with this idea for a few hours and tend to agree with you.

The problem is that we cannot have a specific implementation for any algorithm parameterized by Field due to a design choice in the Java language itself: generics are implemented by type erasure, which means that after compilation, the generic parameter is lost. The top level interface for "fieldified" integrators is FieldODEIntegrator<T extends RealFieldElement<T>>, so I think that the compiled integration classes only refer to RealFieldElement and cannot separate DerivativeStructure (or the sibling classes Gradient, UnivariateDerivative1 and UnivariateDerivative2) from Dfp, Decimal64 or Tuple for example.

One solution could be to just assume step sizes in adaptive stepsize integrators are never fields and must always be treated as simple primitive double. I think it does make sense, for example in the Tuple case, you already use the status of the first tuple component to decide how step is computed, the other components are ignored and will be propagated at the same rate as the first one. This move is an incompatible move (but anyway we have other incompatible move waiting to be published,also concerning fields and Complex) so publishing Hipparchus 2.0 is already something we have in mind.

I would however like to have the opinions of other Hipparchus developers. @wardev, @psteitz, @BryanCazabonne, @oertl , @HankG, do you have any thought about this?

Hi,

I didn't think that dealing with the step size in field or real would have such a big impact on the calculation of derivatives. In any case, I am particularly interested in this problem because the issue #587 of Orekit has a significant impact on the accuracy of the orbit determination. I recently studied the difference in accuracy of a state transition matrix computation between a finite differences method and automatic differentiation method using the Grandient class. (See Table 4 and Table 5 of the following reference)

I noticed that by adding the atmospheric drag, the difference increases strongly (4 orders of maginute). This behavior is certainly related to the present issure and therefore, correcting it would allow to improve considerably the results in Orekit.

So I agree with the proposal to assume step sizes in adaptive step size integrators are never fields and must always be treated as simple primitive double.

Best, Bryan

A few more thoughts on this: step size by itself should remain a field.

For all steps, current time at step k is computed using expression tₖ = tₖ₋₁ + hₖ₋₁ where the hₖ₋₁ is the current step size.

When we set up DerivativeStructure so that t₀ is one of the free parameters (i.e. we compute partial derivatives with respect to initial time), then all times tₖ have a partial derivative ∂tₖ/∂t₀ = 1, even if all hᵢ are simple primitive doubles, which is fine.

Last step is special. It is always truncated to match end time, so hₙ = tₑ - tₙ When we set up DerivativeStructure so that tₑ is one of the free parameters (i.e. we compute partial derivatives with respect to final time), then we must compute hₙ as a derivative so ∂hₙ/∂tₑ = 1 and ∂hₙ/∂t₀ = -1 and when we use the first expression to compute last time (i.e. tₙ₊₁ = tₙ + hₙ), we properly get back tₙ₊₁=tₑ and have ∂tₙ₊₁/∂tₑ = 1 and ∂tₙ₊₁/∂t₀ = 0 (i.e. final time depends only on itself).

If we force step size to be a primitive double, we lose the partial derivatives when we compute hₙ = tₑ - tₙ and we get ∂tₙ₊₁/∂tₑ = 0 and ∂tₙ₊₁/∂t₀ = 1, which is wrong.

Fortunately, this problem is caught by unit tests, this is how I noticed it.

I have implemented the following changes:

initializeStep now returns a double
errorEstimation (several implementations) now returns a double This was done on all field versions of adaptive step size integrators, i.e. embedded Runge-Kutta types and multistep (Adams) types. We don't have a field implementation of Gragg-Bulirsch-Stoer yet.

This has been done in a dedicated issue-118 branch. This branch will be merged on master only after more review.

Beware as since this change is an incompatible one (we change signature of public methods), I had to change pom version to 2.0-SNAPSHOT, so everyone wanting to test this must update the dependencies of downstream projects.

@afossa could you check this solve you issue? @BryanCazabonne could you look if this solve issue 587

Thank you @maisonobe and @BryanCazabonne for your interest.

Regarding your previous observation on the computation of the partials w.r.t. the initial/final time, I was expecting that a dedicated handling of the step-size computation was needed in this case, but I hadn't look into that yet.

I will test your changes on my problem and come back with the results soon.

I will check it today. I will also verify if it improves the calculation of the state transition matrix for my semi-analytical simulation using LEO satellite.

I tested your changes and I can confirm they solve my issue for the test case posted above and also for more complex cases.

Indeed, your modification to the initializeStep and errorEstimation methods are equivalent to what I implemented locally to obtain the results for the "truncated error" case.

I also run the test case for issue #587 in Orekit and I obtain errors between DerivativeStructure and PartialDerivatives based partials in the order of 1e-13 without drag and 1e-5 with drag instead of 1e-6 and 1e+14 before.

I tested your changes and I can confirm they solve my issue for the test case posted above and also for more complex cases.

So it seems the fix works, great!

Indeed, your modification to the initializeStep and errorEstimation methods are equivalent to what I implemented locally to obtain the results for the "truncated error" case.

Yes, of course. I started from your proposal, studied it and just generalized it. I have to add you name in the contributor section of the pom and in the changes file, sorry to have forgotten that.

I also run the test case for issue #587 in Orekit and I obtain errors between DerivativeStructure and PartialDerivatives based partials in the order of 1e-13 without drag and 1e-5 with drag instead of 1e-6 and 1e+14 before.

Nice!

I confirm, it fixes issue #587 of Orekit.

However, it doesn't fix the issue I have in STM calculation when drag is added. Probably, we should not forget the first proposed fix about the calculation of the derivatives of tabulated functions.

Hipparchus-Math / hipparchus

Step size computation in AdaptiveStepsizeFieldIntegrator #118