[MATH-338] Wrong parameter for first step size guess for Embedded Runge Kutta methods

• Type: Bug
• Status: Closed
• Priority: Major
• Resolution: Fixed
• Affects Version/s: 2.0
• Fix Version/s: 2.1

Environment

Eclipse sous Red Hat 5

• Assignee: Unassigned
• Reporter: Morand Vincent
• Created: 28/Jan/10 10:59
• Updated: 03/Apr/10 03:54
• Resolved: 28/Jan/10 15:17

Description

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