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

4 messages
Open this post in threaded view
|

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

 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.
Open this post in threaded view
|

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

 [ 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.