[jira] Created: (MATH-338) Wrong parameter for first step size guess for Embedded Runge Kutta methods

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[jira] Created: (MATH-338) Wrong parameter for first step size guess for Embedded Runge Kutta methods

ASF GitHub Bot (Jira)
Wrong parameter for first step size guess for Embedded Runge Kutta methods
--------------------------------------------------------------------------

                 Key: MATH-338
                 URL: https://issues.apache.org/jira/browse/MATH-338
             Project: Commons Math
          Issue Type: Bug
    Affects Versions: 2.0
         Environment: Eclipse sous Red Hat 5
            Reporter: Morand Vincent


In a space application using DOP853 i detected what seems to be a bad parameter in the call to the method  initializeStep of class AdaptiveStepsizeIntegrator.

Here, DormandPrince853Integrator is a subclass for EmbeddedRungeKuttaIntegrator which perform the call to initializeStep at the beginning of its method integrate(...)

The problem comes from the array "scale" that is used as a parameter in the call off initializeStep(..)

Following the theory described by Hairer in his book "Solving Ordinary Differential Equations 1 : Nonstiff Problems", the scaling should be :

sci = Atol i + |y0i| * Rtoli

Whereas EmbeddedRungeKuttaIntegrator uses :  sci = Atoli

Note that the Gragg-Bulirsch-Stoer integrator uses the good implementation "sci = Atol i + |y0i| * Rtoli  " when he performs the call to the same method initializeStep(..)

In the method initializeStep, the error leads to a wrong step size h used to perform an  Euler step. Most of the time it is unvisible for the user.
But in my space application the Euler step with this wrong step size h (much bigger than it should be)  makes an exception occur (my satellite hits the ground...)


To fix the bug, one should use the same algorithm as in the rescale method in GraggBulirschStoerIntegrator
For exemple :

 final double[] scale= new double[y0.length];;
         
          if (vecAbsoluteTolerance == null) {
              for (int i = 0; i < scale.length; ++i) {
                final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
                scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * yi;
              }
            } else {
              for (int i = 0; i < scale.length; ++i) {
                final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
                scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yi;
              }
            }
         
          hNew = initializeStep(equations, forward, getOrder(), scale,
                           stepStart, y, yDotK[0], yTmp, yDotK[1]);



Sorry for the length of this message, looking forward to hearing from you soon

Vincent Morand






--
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.

Reply | Threaded
Open this post in threaded view
|

[jira] Resolved: (MATH-338) Wrong parameter for first step size guess for Embedded Runge Kutta methods

ASF GitHub Bot (Jira)

     [ https://issues.apache.org/jira/browse/MATH-338?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Luc Maisonobe resolved MATH-338.
--------------------------------

    Resolution: Fixed

Fixed in subversion repository as of r904112.
Note that I have changed slightly the fix you proposed: the call to Math.max was not needed because both arguments were the same (in GBS integrator, they are different).
Note that instead of letting the integrator guess the first step by itself, you can provide it yourself by calling setInitialStepSize. This setting must be done before the call to integrate, which is called by Orekit propagate method if you happen to use Orekit for your application ;-) For such applications, an initial step of the order of magnitude of 1/100 of the keplerian period is a fair bet, it will be adjusted by the integrator if inconsistent with your accuracy settings.

Thanks for the report.

> Wrong parameter for first step size guess for Embedded Runge Kutta methods
> --------------------------------------------------------------------------
>
>                 Key: MATH-338
>                 URL: https://issues.apache.org/jira/browse/MATH-338
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 2.0
>         Environment: Eclipse sous Red Hat 5
>            Reporter: Morand Vincent
>
> In a space application using DOP853 i detected what seems to be a bad parameter in the call to the method  initializeStep of class AdaptiveStepsizeIntegrator.
> Here, DormandPrince853Integrator is a subclass for EmbeddedRungeKuttaIntegrator which perform the call to initializeStep at the beginning of its method integrate(...)
> The problem comes from the array "scale" that is used as a parameter in the call off initializeStep(..)
> Following the theory described by Hairer in his book "Solving Ordinary Differential Equations 1 : Nonstiff Problems", the scaling should be :
> sci = Atol i + |y0i| * Rtoli
> Whereas EmbeddedRungeKuttaIntegrator uses :  sci = Atoli
> Note that the Gragg-Bulirsch-Stoer integrator uses the good implementation "sci = Atol i + |y0i| * Rtoli  " when he performs the call to the same method initializeStep(..)
> In the method initializeStep, the error leads to a wrong step size h used to perform an  Euler step. Most of the time it is unvisible for the user.
> But in my space application the Euler step with this wrong step size h (much bigger than it should be)  makes an exception occur (my satellite hits the ground...)
> To fix the bug, one should use the same algorithm as in the rescale method in GraggBulirschStoerIntegrator
> For exemple :
>  final double[] scale= new double[y0.length];;
>          
>           if (vecAbsoluteTolerance == null) {
>               for (int i = 0; i < scale.length; ++i) {
>                 final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
>                 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * yi;
>               }
>             } else {
>               for (int i = 0; i < scale.length; ++i) {
>                 final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
>                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yi;
>               }
>             }
>          
>           hNew = initializeStep(equations, forward, getOrder(), scale,
>                            stepStart, y, yDotK[0], yTmp, yDotK[1]);
> Sorry for the length of this message, looking forward to hearing from you soon
> Vincent Morand

--
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.

Reply | Threaded
Open this post in threaded view
|

[jira] Closed: (MATH-338) Wrong parameter for first step size guess for Embedded Runge Kutta methods

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

     [ https://issues.apache.org/jira/browse/MATH-338?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Phil Steitz closed MATH-338.
----------------------------


> Wrong parameter for first step size guess for Embedded Runge Kutta methods
> --------------------------------------------------------------------------
>
>                 Key: MATH-338
>                 URL: https://issues.apache.org/jira/browse/MATH-338
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 2.0
>         Environment: Eclipse sous Red Hat 5
>            Reporter: Morand Vincent
>             Fix For: 2.1
>
>
> In a space application using DOP853 i detected what seems to be a bad parameter in the call to the method  initializeStep of class AdaptiveStepsizeIntegrator.
> Here, DormandPrince853Integrator is a subclass for EmbeddedRungeKuttaIntegrator which perform the call to initializeStep at the beginning of its method integrate(...)
> The problem comes from the array "scale" that is used as a parameter in the call off initializeStep(..)
> Following the theory described by Hairer in his book "Solving Ordinary Differential Equations 1 : Nonstiff Problems", the scaling should be :
> sci = Atol i + |y0i| * Rtoli
> Whereas EmbeddedRungeKuttaIntegrator uses :  sci = Atoli
> Note that the Gragg-Bulirsch-Stoer integrator uses the good implementation "sci = Atol i + |y0i| * Rtoli  " when he performs the call to the same method initializeStep(..)
> In the method initializeStep, the error leads to a wrong step size h used to perform an  Euler step. Most of the time it is unvisible for the user.
> But in my space application the Euler step with this wrong step size h (much bigger than it should be)  makes an exception occur (my satellite hits the ground...)
> To fix the bug, one should use the same algorithm as in the rescale method in GraggBulirschStoerIntegrator
> For exemple :
>  final double[] scale= new double[y0.length];;
>          
>           if (vecAbsoluteTolerance == null) {
>               for (int i = 0; i < scale.length; ++i) {
>                 final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
>                 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * yi;
>               }
>             } else {
>               for (int i = 0; i < scale.length; ++i) {
>                 final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
>                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yi;
>               }
>             }
>          
>           hNew = initializeStep(equations, forward, getOrder(), scale,
>                            stepStart, y, yDotK[0], yTmp, yDotK[1]);
> Sorry for the length of this message, looking forward to hearing from you soon
> Vincent Morand

--
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.

Reply | Threaded
Open this post in threaded view
|

[jira] Updated: (MATH-338) Wrong parameter for first step size guess for Embedded Runge Kutta methods

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

     [ https://issues.apache.org/jira/browse/MATH-338?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Phil Steitz updated MATH-338:
-----------------------------

    Fix Version/s: 2.1

> Wrong parameter for first step size guess for Embedded Runge Kutta methods
> --------------------------------------------------------------------------
>
>                 Key: MATH-338
>                 URL: https://issues.apache.org/jira/browse/MATH-338
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 2.0
>         Environment: Eclipse sous Red Hat 5
>            Reporter: Morand Vincent
>             Fix For: 2.1
>
>
> In a space application using DOP853 i detected what seems to be a bad parameter in the call to the method  initializeStep of class AdaptiveStepsizeIntegrator.
> Here, DormandPrince853Integrator is a subclass for EmbeddedRungeKuttaIntegrator which perform the call to initializeStep at the beginning of its method integrate(...)
> The problem comes from the array "scale" that is used as a parameter in the call off initializeStep(..)
> Following the theory described by Hairer in his book "Solving Ordinary Differential Equations 1 : Nonstiff Problems", the scaling should be :
> sci = Atol i + |y0i| * Rtoli
> Whereas EmbeddedRungeKuttaIntegrator uses :  sci = Atoli
> Note that the Gragg-Bulirsch-Stoer integrator uses the good implementation "sci = Atol i + |y0i| * Rtoli  " when he performs the call to the same method initializeStep(..)
> In the method initializeStep, the error leads to a wrong step size h used to perform an  Euler step. Most of the time it is unvisible for the user.
> But in my space application the Euler step with this wrong step size h (much bigger than it should be)  makes an exception occur (my satellite hits the ground...)
> To fix the bug, one should use the same algorithm as in the rescale method in GraggBulirschStoerIntegrator
> For exemple :
>  final double[] scale= new double[y0.length];;
>          
>           if (vecAbsoluteTolerance == null) {
>               for (int i = 0; i < scale.length; ++i) {
>                 final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
>                 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * yi;
>               }
>             } else {
>               for (int i = 0; i < scale.length; ++i) {
>                 final double yi = Math.max(Math.abs(y0[i]), Math.abs(y0[i]));
>                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yi;
>               }
>             }
>          
>           hNew = initializeStep(equations, forward, getOrder(), scale,
>                            stepStart, y, yDotK[0], yTmp, yDotK[1]);
> Sorry for the length of this message, looking forward to hearing from you soon
> Vincent Morand

--
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.