commons-issues mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From "Morand Vincent (JIRA)" <j...@apache.org>
Subject [jira] Created: (MATH-338) Wrong parameter for first step size guess for Embedded Runge Kutta methods
Date Thu, 28 Jan 2010 10:59:34 GMT
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.


Mime
View raw message