You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2011/01/20 21:55:01 UTC
svn commit: r1061507 [1/2] - in /commons/proper/math/branches/MATH_2_X/src:
main/java/org/apache/commons/math/ode/
main/java/org/apache/commons/math/ode/events/
main/java/org/apache/commons/math/ode/nonstiff/
main/java/org/apache/commons/math/ode/sampl...
Author: luc
Date: Thu Jan 20 20:55:00 2011
New Revision: 1061507
URL: http://svn.apache.org/viewvc?rev=1061507&view=rev
Log:
separate discrete event detection from adaptive step size handling in ODE solvers,
thus improving robustness, maintainability and speed
JIRA: MATH-484
Modified:
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/EventState.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblem4.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/sampling/DummyStepInterpolatorTest.java
commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java Thu Jan 20 20:55:00 2011
@@ -20,15 +20,22 @@ package org.apache.commons.math.ode;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
+import java.util.List;
+import java.util.SortedSet;
+import java.util.TreeSet;
+import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MaxEvaluationsExceededException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.events.EventException;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.events.EventState;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
+import org.apache.commons.math.util.MathUtils;
/**
* Base class managing common boilerplate for all integrators.
@@ -46,8 +53,17 @@ public abstract class AbstractIntegrator
/** Current stepsize. */
protected double stepSize;
- /** Events handlers manager. */
- protected CombinedEventsManager eventsHandlersManager;
+ /** Indicator for last step. */
+ protected boolean isLastStep;
+
+ /** Indicator that a state or derivative reset was triggered by some event. */
+ protected boolean resetOccurred;
+
+ /** Events states. */
+ protected Collection<EventState> eventsStates;
+
+ /** Initialization indicator of events states. */
+ protected boolean statesInitialized;
/** Name of the method. */
private final String name;
@@ -69,7 +85,8 @@ public abstract class AbstractIntegrator
stepHandlers = new ArrayList<StepHandler>();
stepStart = Double.NaN;
stepSize = Double.NaN;
- eventsHandlersManager = new CombinedEventsManager();
+ eventsStates = new ArrayList<EventState>();
+ statesInitialized = false;
setMaxEvaluations(-1);
resetEvaluations();
}
@@ -101,22 +118,25 @@ public abstract class AbstractIntegrator
}
/** {@inheritDoc} */
- public void addEventHandler(final EventHandler function,
+ public void addEventHandler(final EventHandler handler,
final double maxCheckInterval,
final double convergence,
final int maxIterationCount) {
- eventsHandlersManager.addEventHandler(function, maxCheckInterval,
- convergence, maxIterationCount);
+ eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount));
}
/** {@inheritDoc} */
public Collection<EventHandler> getEventHandlers() {
- return eventsHandlersManager.getEventsHandlers();
+ final List<EventHandler> list = new ArrayList<EventHandler>();
+ for (EventState state : eventsStates) {
+ list.add(state.getEventHandler());
+ }
+ return Collections.unmodifiableCollection(list);
}
/** {@inheritDoc} */
public void clearEventHandlers() {
- eventsHandlersManager.clearEventsHandlers();
+ eventsStates.clear();
}
/** Check if one of the step handlers requires dense output.
@@ -185,6 +205,109 @@ public abstract class AbstractIntegrator
equations.computeDerivatives(t, y, yDot);
}
+ /** Accept a step, triggering events and step handlers.
+ * @param interpolator step interpolator
+ * @param handlers step handlers
+ * @param y state vector at step end time, must be reset if an event
+ * asks for resetting or if an events stops integration during the step
+ * @param yDot placeholder array where to put the time derivative of the state vector
+ * @param tEnd final integration time
+ * @return time at end of step
+ * @exception IntegratorException if the value of one event state cannot be evaluated
+ */
+ protected double acceptStep(final AbstractStepInterpolator interpolator,
+ final Collection<StepHandler> handlers,
+ final double[] y,
+ final double[] yDot, final double tEnd)
+ throws IntegratorException {
+
+ try {
+ double previousT = interpolator.getGlobalPreviousTime();
+ final double currentT = interpolator.getGlobalCurrentTime();
+ resetOccurred = false;
+
+ // initialize the events states if needed
+ if (! statesInitialized) {
+ for (EventState state : eventsStates) {
+ state.reinitializeBegin(interpolator);
+ }
+ statesInitialized = true;
+ }
+
+ // find all events that occur during the step
+ SortedSet<EventState> occuringEvents = new TreeSet<EventState>();
+ for (final EventState state : eventsStates) {
+ if (state.evaluateStep(interpolator)) {
+ // the event occurs during the current step
+ occuringEvents.add(state);
+ }
+ }
+
+ // handle the events chronologically
+ for (final EventState state : occuringEvents) {
+
+ // restrict the interpolator to the first part of the step, up to the event
+ final double eventT = state.getEventTime();
+ interpolator.setSoftBounds(previousT, eventT);
+
+ // trigger the event
+ interpolator.setInterpolatedTime(eventT);
+ final double[] eventY = interpolator.getInterpolatedState();
+ state.stepAccepted(eventT, eventY);
+ isLastStep = state.stop();
+
+ // handle the first part of the step, up to the event
+ for (final StepHandler handler : stepHandlers) {
+ handler.handleStep(interpolator, isLastStep);
+ }
+
+ if (isLastStep) {
+ // the event asked to stop integration
+ System.arraycopy(eventY, 0, y, 0, y.length);
+ return eventT;
+ }
+
+ if (state.reset(eventT, eventY)) {
+ // some event handler has triggered changes that
+ // invalidate the derivatives, we need to recompute them
+ System.arraycopy(eventY, 0, y, 0, y.length);
+ computeDerivatives(eventT, y, yDot);
+ resetOccurred = true;
+ return eventT;
+ }
+
+ // prepare handling of the remaining part of the step
+ previousT = eventT;
+ interpolator.setSoftBounds(eventT, currentT);
+
+ }
+
+ interpolator.setInterpolatedTime(currentT);
+ final double[] currentY = interpolator.getInterpolatedState();
+ for (final EventState state : eventsStates) {
+ state.stepAccepted(currentT, currentY);
+ isLastStep = isLastStep || state.stop();
+ }
+ isLastStep = isLastStep || MathUtils.equals(currentT, tEnd, 1);
+
+ // handle the remaining part of the step, after all events if any
+ for (StepHandler handler : stepHandlers) {
+ handler.handleStep(interpolator, isLastStep);
+ }
+
+ return currentT;
+ } catch (EventException se) {
+ final Throwable cause = se.getCause();
+ if ((cause != null) && (cause instanceof MathUserException)) {
+ throw (MathUserException) cause;
+ }
+ throw new IntegratorException(se);
+ } catch (ConvergenceException ce) {
+ throw new IntegratorException(ce);
+ }
+
+ }
+
/** Perform some sanity checks on the integration parameters.
* @param ode differential equations set
* @param t0 start time
@@ -225,7 +348,9 @@ public abstract class AbstractIntegrator
* @param manager manager containing the user-defined handlers
* @return a new manager containing all the user-defined handlers plus a
* dedicated manager triggering a stop event at entTime
+ * @deprecated as of 2.2, this method is not used any more
*/
+ @Deprecated
protected CombinedEventsManager addEndTimeChecker(final double startTime,
final double endTime,
final CombinedEventsManager manager) {
@@ -243,7 +368,10 @@ public abstract class AbstractIntegrator
return newManager;
}
- /** Specialized event handler to stop integration. */
+ /** Specialized event handler to stop integration.
+ * @deprecated as of 2.2, this class is not used anymore
+ */
+ @Deprecated
private static class EndTimeChecker implements EventHandler {
/** Desired end time. */
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java Thu Jan 20 20:55:00 2011
@@ -73,10 +73,8 @@ public interface ODEIntegrator {
* @see #getEventHandlers()
* @see #clearEventHandlers()
*/
- void addEventHandler(EventHandler handler,
- double maxCheckInterval,
- double convergence,
- int maxIterationCount);
+ void addEventHandler(EventHandler handler, double maxCheckInterval,
+ double convergence, int maxIterationCount);
/** Get all the event handlers that have been added to the integrator.
* @return an unmodifiable collection of the added events handlers
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java Thu Jan 20 20:55:00 2011
@@ -33,8 +33,9 @@ import org.apache.commons.math.ode.sampl
* @see EventState
* @version $Revision$ $Date$
* @since 1.2
+ * @deprecated as of 2.2, this class is not used anymore
*/
-
+@Deprecated
public class CombinedEventsManager {
/** Events states. */
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/EventState.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/EventState.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/EventState.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/events/EventState.java Thu Jan 20 20:55:00 2011
@@ -18,9 +18,9 @@
package org.apache.commons.math.ode.events;
import org.apache.commons.math.ConvergenceException;
-import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.solvers.BrentSolver;
+import org.apache.commons.math.exception.MathInternalError;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.ode.sampling.StepInterpolator;
import org.apache.commons.math.util.FastMath;
@@ -39,7 +39,7 @@ import org.apache.commons.math.util.Fast
* @version $Revision$ $Date$
* @since 1.2
*/
-public class EventState {
+public class EventState implements Comparable<EventState> {
/** Event handler. */
private final EventHandler handler;
@@ -187,8 +187,7 @@ public class EventState {
/** Evaluate the impact of the proposed step on the event handler.
* @param interpolator step interpolator for the proposed step
* @return true if the event handler triggers an event before
- * the end of the proposed step (this implies the step should be
- * rejected)
+ * the end of the proposed step
* @exception MathUserException if the interpolator fails to
* compute the switching function somewhere within the step
* @exception EventException if the switching function
@@ -202,16 +201,21 @@ public class EventState {
forward = interpolator.isForward();
final double t1 = interpolator.getCurrentTime();
- final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(t1 - t0) / maxCheckInterval));
- final double h = (t1 - t0) / n;
+ if (FastMath.abs(t1 - t0) < convergence) {
+ // we cannot do anything on such a small step, don't trigger any events
+ return false;
+ }
+ final double start = forward ? (t0 + convergence) : t0 - convergence;
+ final double dt = t1 - start;
+ final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
+ final double h = dt / n;
double ta = t0;
double ga = g0;
- double tb = t0 + (interpolator.isForward() ? convergence : -convergence);
for (int i = 0; i < n; ++i) {
// evaluate handler value at the end of the substep
- tb += h;
+ final double tb = start + (i + 1) * h;
interpolator.setInterpolatedTime(tb);
final double gb = handler.g(tb, interpolator.getInterpolatedState());
@@ -219,7 +223,22 @@ public class EventState {
if (g0Positive ^ (gb >= 0)) {
// there is a sign change: an event is expected during this step
- if (ga * gb > 0) {
+ // variation direction, with respect to the integration direction
+ increasing = gb >= ga;
+
+ final UnivariateRealFunction f = new UnivariateRealFunction() {
+ public double value(final double t) throws MathUserException {
+ try {
+ interpolator.setInterpolatedTime(t);
+ return handler.g(t, interpolator.getInterpolatedState());
+ } catch (EventException e) {
+ throw new MathUserException(e);
+ }
+ }
+ };
+ final BrentSolver solver = new BrentSolver(convergence);
+
+ if (ga * gb >= 0) {
// this is a corner case:
// - there was an event near ta,
// - there is another event between ta and tb
@@ -230,49 +249,33 @@ public class EventState {
final double epsilon = (forward ? 0.25 : -0.25) * convergence;
for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
ta += epsilon;
- interpolator.setInterpolatedTime(ta);
- ga = handler.g(ta, interpolator.getInterpolatedState());
+ ga = f.value(ta);
}
if (ga * gb > 0) {
// this should never happen
- throw MathRuntimeException.createInternalError(null);
+ throw new MathInternalError();
}
}
- // variation direction, with respect to the integration direction
- increasing = gb >= ga;
-
- final UnivariateRealFunction f = new UnivariateRealFunction() {
- public double value(final double t) throws MathUserException {
- try {
- interpolator.setInterpolatedTime(t);
- return handler.g(t, interpolator.getInterpolatedState());
- } catch (EventException e) {
- throw new MathUserException(e);
- }
- }
- };
- final BrentSolver solver = new BrentSolver(maxIterationCount, convergence);
- final double root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
- if ((FastMath.abs(root - ta) <= convergence) &&
- (FastMath.abs(root - previousEventTime) <= convergence)) {
+ final double root = (ta <= tb) ?
+ solver.solve(maxIterationCount, f, ta, tb) :
+ solver.solve(maxIterationCount, f, tb, ta);
+
+ if ((!Double.isNaN(previousEventTime)) &&
+ (FastMath.abs(root - ta) <= convergence) &&
+ (FastMath.abs(root - previousEventTime) <= convergence)) {
// we have either found nothing or found (again ?) a past event, we simply ignore it
ta = tb;
ga = gb;
} else if (Double.isNaN(previousEventTime) ||
(FastMath.abs(previousEventTime - root) > convergence)) {
pendingEventTime = root;
- if (pendingEvent && (FastMath.abs(t1 - pendingEventTime) <= convergence)) {
- // we were already waiting for this event which was
- // found during a previous call for a step that was
- // rejected, this step must now be accepted since it
- // properly ends exactly at the event occurrence
- return false;
- }
- // either we were not waiting for the event or it has
- // moved in such a way the step cannot be accepted
pendingEvent = true;
return true;
+ } else {
+ // no sign change: there is no event for now
+ ta = tb;
+ ga = gb;
}
} else {
@@ -321,7 +324,7 @@ public class EventState {
t0 = t;
g0 = handler.g(t, y);
- if (pendingEvent) {
+ if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
// force the sign to its value "just after the event"
previousEventTime = t;
g0Positive = increasing;
@@ -352,7 +355,7 @@ public class EventState {
public boolean reset(final double t, final double[] y)
throws EventException {
- if (! pendingEvent) {
+ if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
return false;
}
@@ -367,4 +370,21 @@ public class EventState {
}
+ /** Compare the instance with another event state.
+ * <p>
+ * Event state ordering is based on occurrence time within the last
+ * evaluated step. If no event occurs during the step, a time arbitrarily
+ * set to positive infinity is used.
+ * </p>
+ * @param state other event state to compare the instance to
+ * @return a negative integer, zero, or a positive integer as the event
+ * occurs before, simultaneous, or after the specified event of the
+ * specified state.
+ */
+ public int compareTo(final EventState state) {
+ final double instanceTime = pendingEvent ? pendingEventTime : Double.POSITIVE_INFINITY;
+ final double otherTime = state.pendingEvent ? state.pendingEventTime : Double.POSITIVE_INFINITY;
+ return Double.compare(instanceTime, otherTime);
+ }
+
}
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Thu Jan 20 20:55:00 2011
@@ -21,7 +21,6 @@ import org.apache.commons.math.exception
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
@@ -202,19 +201,16 @@ public class AdamsBashforthIntegrator ex
System.arraycopy(y0, 0, y, 0, n);
}
final double[] yDot = new double[n];
- final double[] yTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
interpolator.reinitialize(y, forward);
- final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();
- interpolatorTmp.reinitialize(yTmp, forward);
// set up integration control objects
for (StepHandler handler : stepHandlers) {
handler.reset();
}
- CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
+ statesInitialized = false;
// compute the initial Nordsieck vector using the configured starter integrator
start(t0, y, t);
@@ -226,14 +222,12 @@ public class AdamsBashforthIntegrator ex
double hNew = stepSize;
interpolator.rescale(hNew);
- boolean lastStep = false;
- while (!lastStep) {
+ // main integration loop
+ isLastStep = false;
+ do {
- // shift all data
- interpolator.shift();
-
- double error = 0;
- for (boolean loop = true; loop;) {
+ double error = 10;
+ while (error >= 1.0) {
stepSize = hNew;
@@ -249,92 +243,51 @@ public class AdamsBashforthIntegrator ex
}
error = FastMath.sqrt(error / mainSetDimension);
- if (error <= 1.0) {
-
- // predict a first estimate of the state at step end
- final double stepEnd = stepStart + stepSize;
- interpolator.setInterpolatedTime(stepEnd);
- System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
-
- // evaluate the derivative
- computeDerivatives(stepEnd, yTmp, yDot);
-
- // update Nordsieck vector
- final double[] predictedScaled = new double[y0.length];
- for (int j = 0; j < y0.length; ++j) {
- predictedScaled[j] = stepSize * yDot[j];
- }
- final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
- updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
-
- // discrete events handling
- interpolatorTmp.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
- interpolatorTmp.storeTime(stepStart);
- interpolatorTmp.shift();
- interpolatorTmp.storeTime(stepEnd);
- if (manager.evaluateStep(interpolatorTmp)) {
- final double dt = manager.getEventTime() - stepStart;
- if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
- // we cannot simply truncate the step, reject the current computation
- // and let the loop compute another state with the truncated step.
- // it is so small (much probably exactly 0 due to limited accuracy)
- // that the code above would fail handling it.
- // So we set up an artificial 0 size step by copying states
- interpolator.storeTime(stepStart);
- System.arraycopy(y, 0, yTmp, 0, y0.length);
- hNew = 0;
- stepSize = 0;
- loop = false;
- } else {
- // reject the step to match exactly the next switch time
- hNew = dt;
- interpolator.rescale(hNew);
- }
- } else {
- // accept the step
- scaled = predictedScaled;
- nordsieck = nordsieckTmp;
- interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
- loop = false;
- }
-
- } else {
+ if (error >= 1.0) {
// reject the step and attempt to reduce error by stepsize control
final double factor = computeStepGrowShrinkFactor(error);
hNew = filterStep(stepSize * factor, forward, false);
interpolator.rescale(hNew);
- }
-
- }
- // the step has been accepted (may have been truncated)
- final double nextStep = stepStart + stepSize;
- System.arraycopy(yTmp, 0, y, 0, n);
- interpolator.storeTime(nextStep);
- manager.stepAccepted(nextStep, y);
- lastStep = manager.stop();
-
- // provide the step data to the step handler
- for (StepHandler handler : stepHandlers) {
- interpolator.setInterpolatedTime(nextStep);
- handler.handleStep(interpolator, lastStep);
+ }
}
- stepStart = nextStep;
- if (!lastStep && manager.reset(stepStart, y)) {
+ // predict a first estimate of the state at step end
+ final double stepEnd = stepStart + stepSize;
+ interpolator.shift();
+ interpolator.setInterpolatedTime(stepEnd);
+ System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
- // some events handler has triggered changes that
- // invalidate the derivatives, we need to restart from scratch
- start(stepStart, y, t);
- interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
+ // evaluate the derivative
+ computeDerivatives(stepEnd, y, yDot);
+ // update Nordsieck vector
+ final double[] predictedScaled = new double[y0.length];
+ for (int j = 0; j < y0.length; ++j) {
+ predictedScaled[j] = stepSize * yDot[j];
}
-
- if (! lastStep) {
- // in some rare cases we may get here with stepSize = 0, for example
- // when an event occurs at integration start, reducing the first step
- // to zero; we have to reset the step to some safe non zero value
- stepSize = filterStep(stepSize, forward, true);
+ final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
+ updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
+ interpolator.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
+
+ // discrete events handling
+ interpolator.storeTime(stepEnd);
+ stepStart = acceptStep(interpolator, stepHandlers, y, yDot, t);
+ scaled = predictedScaled;
+ nordsieck = nordsieckTmp;
+ interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
+
+ if (!isLastStep) {
+
+ // prepare next step
+ interpolator.storeTime(stepStart);
+
+ if (resetOccurred) {
+ // some events handler has triggered changes that
+ // invalidate the derivatives, we need to restart from scratch
+ start(stepStart, y, t);
+ interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
+ }
// stepsize control for next step
final double factor = computeStepGrowShrinkFactor(error);
@@ -342,14 +295,21 @@ public class AdamsBashforthIntegrator ex
final double nextT = stepStart + scaledH;
final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
hNew = filterStep(scaledH, forward, nextIsLast);
+
+ final double filteredNextT = stepStart + hNew;
+ final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
+ if (filteredNextIsLast) {
+ hNew = t - stepStart;
+ }
+
interpolator.rescale(hNew);
+
}
- }
+ } while (!isLastStep);
- final double stopTime = stepStart;
- stepStart = Double.NaN;
- stepSize = Double.NaN;
+ final double stopTime = stepStart;
+ resetInternalState();
return stopTime;
}
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Thu Jan 20 20:55:00 2011
@@ -24,7 +24,6 @@ import org.apache.commons.math.linear.Ar
import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
@@ -220,19 +219,18 @@ public class AdamsMoultonIntegrator exte
}
final double[] yDot = new double[y0.length];
final double[] yTmp = new double[y0.length];
+ final double[] predictedScaled = new double[y0.length];
+ Array2DRowRealMatrix nordsieckTmp = null;
// set up two interpolators sharing the integrator arrays
final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
interpolator.reinitialize(y, forward);
- final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();
- interpolatorTmp.reinitialize(yTmp, forward);
// set up integration control objects
for (StepHandler handler : stepHandlers) {
handler.reset();
}
- CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
-
+ statesInitialized = false;
// compute the initial Nordsieck vector using the configured starter integrator
start(t0, y, t);
@@ -242,14 +240,11 @@ public class AdamsMoultonIntegrator exte
double hNew = stepSize;
interpolator.rescale(hNew);
- boolean lastStep = false;
- while (!lastStep) {
+ isLastStep = false;
+ do {
- // shift all data
- interpolator.shift();
-
- double error = 0;
- for (boolean loop = true; loop;) {
+ double error = 10;
+ while (error >= 1.0) {
stepSize = hNew;
@@ -262,96 +257,56 @@ public class AdamsMoultonIntegrator exte
computeDerivatives(stepEnd, yTmp, yDot);
// update Nordsieck vector
- final double[] predictedScaled = new double[y0.length];
for (int j = 0; j < y0.length; ++j) {
predictedScaled[j] = stepSize * yDot[j];
}
- final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
+ nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
// apply correction (C in the PECE sequence)
error = nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
- if (error <= 1.0) {
-
- // evaluate a final estimate of the derivative (second E in the PECE sequence)
- computeDerivatives(stepEnd, yTmp, yDot);
-
- // update Nordsieck vector
- final double[] correctedScaled = new double[y0.length];
- for (int j = 0; j < y0.length; ++j) {
- correctedScaled[j] = stepSize * yDot[j];
- }
- updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
-
- // discrete events handling
- interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
- interpolatorTmp.storeTime(stepStart);
- interpolatorTmp.shift();
- interpolatorTmp.storeTime(stepEnd);
- if (manager.evaluateStep(interpolatorTmp)) {
- final double dt = manager.getEventTime() - stepStart;
- if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
- // we cannot simply truncate the step, reject the current computation
- // and let the loop compute another state with the truncated step.
- // it is so small (much probably exactly 0 due to limited accuracy)
- // that the code above would fail handling it.
- // So we set up an artificial 0 size step by copying states
- interpolator.storeTime(stepStart);
- System.arraycopy(y, 0, yTmp, 0, y0.length);
- hNew = 0;
- stepSize = 0;
- loop = false;
- } else {
- // reject the step to match exactly the next switch time
- hNew = dt;
- interpolator.rescale(hNew);
- }
- } else {
- // accept the step
- scaled = correctedScaled;
- nordsieck = nordsieckTmp;
- interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
- loop = false;
- }
-
- } else {
+ if (error >= 1.0) {
// reject the step and attempt to reduce error by stepsize control
final double factor = computeStepGrowShrinkFactor(error);
hNew = filterStep(stepSize * factor, forward, false);
interpolator.rescale(hNew);
}
-
}
- // the step has been accepted (may have been truncated)
- final double nextStep = stepStart + stepSize;
- System.arraycopy(yTmp, 0, y, 0, n);
- interpolator.storeTime(nextStep);
- manager.stepAccepted(nextStep, y);
- lastStep = manager.stop();
-
- // provide the step data to the step handler
- for (StepHandler handler : stepHandlers) {
- interpolator.setInterpolatedTime(nextStep);
- handler.handleStep(interpolator, lastStep);
+ // evaluate a final estimate of the derivative (second E in the PECE sequence)
+ final double stepEnd = stepStart + stepSize;
+ computeDerivatives(stepEnd, yTmp, yDot);
+
+ // update Nordsieck vector
+ final double[] correctedScaled = new double[y0.length];
+ for (int j = 0; j < y0.length; ++j) {
+ correctedScaled[j] = stepSize * yDot[j];
}
- stepStart = nextStep;
-
- if (!lastStep && manager.reset(stepStart, y)) {
+ updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
- // some events handler has triggered changes that
- // invalidate the derivatives, we need to restart from scratch
- start(stepStart, y, t);
- interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
-
- }
+ // discrete events handling
+ System.arraycopy(yTmp, 0, y, 0, n);
+ interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
+ interpolator.storeTime(stepStart);
+ interpolator.shift();
+ interpolator.storeTime(stepEnd);
+ stepStart = acceptStep(interpolator, stepHandlers, y, yDot, t);
+ scaled = correctedScaled;
+ nordsieck = nordsieckTmp;
+
+ if (!isLastStep) {
+
+ // prepare next step
+ interpolator.storeTime(stepStart);
+
+ if (resetOccurred) {
+ // some events handler has triggered changes that
+ // invalidate the derivatives, we need to restart from scratch
+ start(stepStart, y, t);
+ interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
- if (! lastStep) {
- // in some rare cases we may get here with stepSize = 0, for example
- // when an event occurs at integration start, reducing the first step
- // to zero; we have to reset the step to some safe non zero value
- stepSize = filterStep(stepSize, forward, true);
+ }
// stepsize control for next step
final double factor = computeStepGrowShrinkFactor(error);
@@ -359,10 +314,17 @@ public class AdamsMoultonIntegrator exte
final double nextT = stepStart + scaledH;
final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
hNew = filterStep(scaledH, forward, nextIsLast);
+
+ final double filteredNextT = stepStart + hNew;
+ final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
+ if (filteredNextIsLast) {
+ hNew = t - stepStart;
+ }
+
interpolator.rescale(hNew);
}
- }
+ } while (!isLastStep);
final double stopTime = stepStart;
stepStart = Double.NaN;
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java Thu Jan 20 20:55:00 2011
@@ -21,7 +21,6 @@ import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
-import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.sampling.StepInterpolator;
@@ -396,6 +395,7 @@ class DormandPrince853StepInterpolator
double s;
final double[] yTmp = new double[currentState.length];
+ final double pT = getGlobalPreviousTime();
// k14
for (int j = 0; j < currentState.length; ++j) {
@@ -404,7 +404,7 @@ class DormandPrince853StepInterpolator
K14_11 * yDotK[10][j] + K14_12 * yDotK[11][j] + K14_13 * yDotK[12][j];
yTmp[j] = currentState[j] + h * s;
}
- integrator.computeDerivatives(previousTime + C14 * h, yTmp, yDotKLast[0]);
+ integrator.computeDerivatives(pT + C14 * h, yTmp, yDotKLast[0]);
// k15
for (int j = 0; j < currentState.length; ++j) {
@@ -414,7 +414,7 @@ class DormandPrince853StepInterpolator
K15_14 * yDotKLast[0][j];
yTmp[j] = currentState[j] + h * s;
}
- integrator.computeDerivatives(previousTime + C15 * h, yTmp, yDotKLast[1]);
+ integrator.computeDerivatives(pT + C15 * h, yTmp, yDotKLast[1]);
// k16
for (int j = 0; j < currentState.length; ++j) {
@@ -424,7 +424,7 @@ class DormandPrince853StepInterpolator
K16_14 * yDotKLast[0][j] + K16_15 * yDotKLast[1][j];
yTmp[j] = currentState[j] + h * s;
}
- integrator.computeDerivatives(previousTime + C16 * h, yTmp, yDotKLast[2]);
+ integrator.computeDerivatives(pT + C16 * h, yTmp, yDotKLast[2]);
}
@@ -437,7 +437,9 @@ class DormandPrince853StepInterpolator
// save the local attributes
finalizeStep();
} catch (MathUserException e) {
- throw MathRuntimeException.createIOException(e);
+ IOException ioe = new IOException(e.getLocalizedMessage());
+ ioe.initCause(e);
+ throw ioe;
}
final int dimension = (currentState == null) ? -1 : currentState.length;
out.writeInt(dimension);
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java Thu Jan 20 20:55:00 2011
@@ -20,7 +20,6 @@ package org.apache.commons.math.ode.nons
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
@@ -206,11 +205,12 @@ public abstract class EmbeddedRungeKutta
System.arraycopy(y0, 0, y, 0, y0.length);
}
final double[][] yDotK = new double[stages][y0.length];
- final double[] yTmp = new double[y0.length];
+ final double[] yTmp = new double[y0.length];
+ final double[] yDotTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
AbstractStepInterpolator interpolator;
- if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
+ if (requiresDenseOutput() || (! eventsStates.isEmpty())) {
final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
rki.reinitialize(this, yTmp, yDotK, forward);
interpolator = rki;
@@ -226,16 +226,17 @@ public abstract class EmbeddedRungeKutta
for (StepHandler handler : stepHandlers) {
handler.reset();
}
- CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
- boolean lastStep = false;
+ statesInitialized = false;
// main integration loop
- while (!lastStep) {
+ isLastStep = false;
+ do {
interpolator.shift();
- double error = 0;
- for (boolean loop = true; loop;) {
+ // iterate over step size, ensuring local normalized error is smaller than 1
+ double error = 10;
+ while (error >= 1.0) {
if (firstTime || !fsal) {
// first stage
@@ -248,11 +249,11 @@ public abstract class EmbeddedRungeKutta
for (int i = 0; i < scale.length; ++i) {
scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]);
}
- } else {
+ } else {
for (int i = 0; i < scale.length; ++i) {
scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
}
- }
+ }
hNew = initializeStep(equations, forward, getOrder(), scale,
stepStart, y, yDotK[0], yTmp, yDotK[1]);
firstTime = false;
@@ -286,83 +287,49 @@ public abstract class EmbeddedRungeKutta
// estimate the error at the end of the step
error = estimateError(yDotK, y, yTmp, stepSize);
- if (error <= 1.0) {
-
- // discrete events handling
- interpolator.storeTime(stepStart + stepSize);
- if (manager.evaluateStep(interpolator)) {
- final double dt = manager.getEventTime() - stepStart;
- if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
- // we cannot simply truncate the step, reject the current computation
- // and let the loop compute another state with the truncated step.
- // it is so small (much probably exactly 0 due to limited accuracy)
- // that the code above would fail handling it.
- // So we set up an artificial 0 size step by copying states
- interpolator.storeTime(stepStart);
- System.arraycopy(y, 0, yTmp, 0, y0.length);
- hNew = 0;
- stepSize = 0;
- loop = false;
- } else {
- // reject the step to match exactly the next switch time
- hNew = dt;
- }
- } else {
- // accept the step
- loop = false;
- }
-
- } else {
+ if (error >= 1.0) {
// reject the step and attempt to reduce error by stepsize control
final double factor =
FastMath.min(maxGrowth,
- FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
+ FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
hNew = filterStep(stepSize * factor, forward, false);
}
}
- // the step has been accepted
- final double nextStep = stepStart + stepSize;
+ // local error is small enough: accept the step, trigger events and step handlers
+ interpolator.storeTime(stepStart + stepSize);
System.arraycopy(yTmp, 0, y, 0, y0.length);
- manager.stepAccepted(nextStep, y);
- lastStep = manager.stop();
+ System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
+ stepStart = acceptStep(interpolator, stepHandlers, y, yDotTmp, t);
- // provide the step data to the step handler
- interpolator.storeTime(nextStep);
- for (StepHandler handler : stepHandlers) {
- handler.handleStep(interpolator, lastStep);
- }
- stepStart = nextStep;
+ if (!isLastStep) {
- if (fsal) {
- // save the last evaluation for the next step
- System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
- }
+ // prepare next step
+ interpolator.storeTime(stepStart);
- if (manager.reset(stepStart, y) && ! lastStep) {
- // some event handler has triggered changes that
- // invalidate the derivatives, we need to recompute them
- computeDerivatives(stepStart, y, yDotK[0]);
- }
+ if (fsal) {
+ // save the last evaluation for the next step
+ System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length);
+ }
+
+ // stepsize control for next step
+ final double factor =
+ FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
+ final double scaledH = stepSize * factor;
+ final double nextT = stepStart + scaledH;
+ final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+ hNew = filterStep(scaledH, forward, nextIsLast);
+
+ final double filteredNextT = stepStart + hNew;
+ final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
+ if (filteredNextIsLast) {
+ hNew = t - stepStart;
+ }
- if (! lastStep) {
- // in some rare cases we may get here with stepSize = 0, for example
- // when an event occurs at integration start, reducing the first step
- // to zero; we have to reset the step to some safe non zero value
- stepSize = filterStep(stepSize, forward, true);
-
- // stepsize control for next step
- final double factor = FastMath.min(maxGrowth,
- FastMath.max(minReduction,
- safety * FastMath.pow(error, exp)));
- final double scaledH = stepSize * factor;
- final double nextT = stepStart + scaledH;
- final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
- hNew = filterStep(scaledH, forward, nextIsLast);
}
- }
+ } while (!isLastStep);
final double stopTime = stepStart;
resetInternalState();
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java Thu Jan 20 20:55:00 2011
@@ -170,7 +170,7 @@ public class GraggBulirschStoerIntegrato
final double scalRelativeTolerance) {
super(METHOD_NAME, minStep, maxStep,
scalAbsoluteTolerance, scalRelativeTolerance);
- denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+ denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
setStabilityCheck(true, -1, -1, -1);
setStepsizeControl(-1, -1, -1, -1);
setOrderControl(-1, -1, -1);
@@ -193,7 +193,7 @@ public class GraggBulirschStoerIntegrato
final double[] vecRelativeTolerance) {
super(METHOD_NAME, minStep, maxStep,
vecAbsoluteTolerance, vecRelativeTolerance);
- denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+ denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
setStabilityCheck(true, -1, -1, -1);
setStepsizeControl(-1, -1, -1, -1);
setOrderControl(-1, -1, -1);
@@ -339,7 +339,7 @@ public class GraggBulirschStoerIntegrato
public void addStepHandler (final StepHandler handler) {
super.addStepHandler(handler);
- denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+ denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
// reinitialize the arrays
initializeArrays();
@@ -353,7 +353,7 @@ public class GraggBulirschStoerIntegrato
final double convergence,
final int maxIterationCount) {
super.addEventHandler(function, maxCheckInterval, convergence, maxIterationCount);
- denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+ denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
// reinitialize the arrays
initializeArrays();
@@ -556,7 +556,7 @@ public class GraggBulirschStoerIntegrato
@Override
public double integrate(final FirstOrderDifferentialEquations equations,
final double t0, final double[] y0, final double t, final double[] y)
- throws MathUserException, IntegratorException {
+ throws MathUserException, IntegratorException {
sanityChecks(equations, t0, y0, t, y);
setEquations(equations);
@@ -594,10 +594,9 @@ public class GraggBulirschStoerIntegrato
System.arraycopy(y0, 0, y, 0, y0.length);
}
- double[] yDot1 = null;
+ double[] yDot1 = new double[y0.length];
double[][] yMidDots = null;
if (denseOutput) {
- yDot1 = new double[y0.length];
yMidDots = new double[1 + 2 * sequence.length][];
for (int j = 0; j < yMidDots.length; ++j) {
yMidDots[j] = new double[y0.length];
@@ -614,13 +613,13 @@ public class GraggBulirschStoerIntegrato
// initial order selection
final double tol =
(vecRelativeTolerance == null) ? scalRelativeTolerance : vecRelativeTolerance[0];
- final double log10R = FastMath.log(FastMath.max(1.0e-10, tol)) / FastMath.log(10.0);
+ final double log10R = FastMath.log10(FastMath.max(1.0e-10, tol));
int targetIter = FastMath.max(1,
FastMath.min(sequence.length - 2,
(int) FastMath.floor(0.5 - 0.6 * log10R)));
// set up an interpolator sharing the integrator arrays
AbstractStepInterpolator interpolator = null;
- if (denseOutput || (! eventsHandlersManager.isEmpty())) {
+ if (denseOutput || (! eventsStates.isEmpty())) {
interpolator = new GraggBulirschStoerStepInterpolator(y, yDot0,
y1, yDot1,
yMidDots, forward);
@@ -635,13 +634,14 @@ public class GraggBulirschStoerIntegrato
boolean previousRejected = false;
boolean firstTime = true;
boolean newStep = true;
- boolean lastStep = false;
boolean firstStepAlreadyComputed = false;
for (StepHandler handler : stepHandlers) {
handler.reset();
}
+ statesInitialized = false;
costPerTimeUnit[0] = 0;
- while (! lastStep) {
+ isLastStep = false;
+ do {
double error;
boolean reject = false;
@@ -656,15 +656,9 @@ public class GraggBulirschStoerIntegrato
}
if (firstTime) {
-
hNew = initializeStep(equations, forward,
2 * targetIter + 1, scale,
stepStart, y, yDot0, yTmp, yTmpDot);
-
- if (! forward) {
- hNew = -hNew;
- }
-
}
newStep = false;
@@ -679,7 +673,7 @@ public class GraggBulirschStoerIntegrato
stepSize = t - stepStart;
}
final double nextT = stepStart + stepSize;
- lastStep = forward ? (nextT >= t) : (nextT <= t);
+ isLastStep = forward ? (nextT >= t) : (nextT <= t);
// iterate over several substep sizes
int k = -1;
@@ -804,7 +798,7 @@ public class GraggBulirschStoerIntegrato
break;
default :
- if ((firstTime || lastStep) && (error <= 1.0)) {
+ if ((firstTime || isLastStep) && (error <= 1.0)) {
loop = false;
}
break;
@@ -816,6 +810,11 @@ public class GraggBulirschStoerIntegrato
}
}
+ if (! reject) {
+ // derivatives at end of step
+ computeDerivatives(stepStart + stepSize, y1, yDot1);
+ }
+
// dense output handling
double hInt = getMaxStep();
if (denseOutput && ! reject) {
@@ -825,9 +824,6 @@ public class GraggBulirschStoerIntegrato
extrapolate(0, j, diagonal, yMidDots[0]);
}
- // derivative at end of step
- computeDerivatives(stepStart + stepSize, y1, yDot1);
-
final int mu = 2 * k - mudif + 3;
for (int l = 0; l < mu; ++l) {
@@ -880,52 +876,21 @@ public class GraggBulirschStoerIntegrato
}
}
- // Discrete events handling
- if (!reject) {
- interpolator.storeTime(stepStart + stepSize);
- if (eventsHandlersManager.evaluateStep(interpolator)) {
- final double dt = eventsHandlersManager.getEventTime() - stepStart;
- if (FastMath.abs(dt) > FastMath.ulp(stepStart)) {
- // reject the step to match exactly the next switch time
- hNew = FastMath.abs(dt);
- reject = true;
- }
- }
- }
-
- }
-
- if (!reject) {
- // we will reuse the slope for the beginning of next step
- firstStepAlreadyComputed = true;
- System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
}
}
if (! reject) {
- // store end of step state
- final double nextStep = stepStart + stepSize;
- System.arraycopy(y1, 0, y, 0, y0.length);
+ // Discrete events handling
+ interpolator.storeTime(stepStart + stepSize);
+ stepStart = acceptStep(interpolator, stepHandlers, y1, yDot1, t);
- eventsHandlersManager.stepAccepted(nextStep, y);
- if (eventsHandlersManager.stop()) {
- lastStep = true;
- }
-
- // provide the step data to the step handler
- interpolator.storeTime(nextStep);
- for (StepHandler handler : stepHandlers) {
- handler.handleStep(interpolator, lastStep);
- }
- stepStart = nextStep;
-
- if (eventsHandlersManager.reset(stepStart, y) && ! lastStep) {
- // some switching function has triggered changes that
- // invalidate the derivatives, we need to recompute them
- firstStepAlreadyComputed = false;
- }
+ // prepare next step
+ interpolator.storeTime(stepStart);
+ System.arraycopy(y1, 0, y, 0, y0.length);
+ System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
+ firstStepAlreadyComputed = true;
int optimalIter;
if (k == 1) {
@@ -987,15 +952,17 @@ public class GraggBulirschStoerIntegrato
firstTime = false;
if (reject) {
- lastStep = false;
+ isLastStep = false;
previousRejected = true;
} else {
previousRejected = false;
}
- }
+ } while (!isLastStep);
- return stepStart;
+ final double stopTime = stepStart;
+ resetInternalState();
+ return stopTime;
}
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java Thu Jan 20 20:55:00 2011
@@ -21,7 +21,6 @@ import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
-import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepInterpolator;
import org.apache.commons.math.util.FastMath;
@@ -310,8 +309,7 @@ class GraggBulirschStoerStepInterpolator
/** {@inheritDoc} */
@Override
protected void computeInterpolatedStateAndDerivatives(final double theta,
- final double oneMinusThetaH)
- throws MathUserException {
+ final double oneMinusThetaH) {
final int dimension = currentState.length;
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java Thu Jan 20 20:55:00 2011
@@ -22,7 +22,6 @@ import org.apache.commons.math.exception
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
@@ -112,11 +111,12 @@ public abstract class RungeKuttaIntegrat
for (int i = 0; i < stages; ++i) {
yDotK [i] = new double[y0.length];
}
- final double[] yTmp = new double[y0.length];
+ final double[] yTmp = new double[y0.length];
+ final double[] yDotTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
AbstractStepInterpolator interpolator;
- if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
+ if (requiresDenseOutput() || (! eventsStates.isEmpty())) {
final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
rki.reinitialize(this, yTmp, yDotK, forward);
interpolator = rki;
@@ -131,90 +131,61 @@ public abstract class RungeKuttaIntegrat
for (StepHandler handler : stepHandlers) {
handler.reset();
}
- CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
- boolean lastStep = false;
+ statesInitialized = false;
// main integration loop
- while (!lastStep) {
+ isLastStep = false;
+ do {
interpolator.shift();
- for (boolean loop = true; loop;) {
+ // first stage
+ computeDerivatives(stepStart, y, yDotK[0]);
- // first stage
- computeDerivatives(stepStart, y, yDotK[0]);
-
- // next stages
- for (int k = 1; k < stages; ++k) {
+ // next stages
+ for (int k = 1; k < stages; ++k) {
for (int j = 0; j < y0.length; ++j) {
- double sum = a[k-1][0] * yDotK[0][j];
- for (int l = 1; l < k; ++l) {
- sum += a[k-1][l] * yDotK[l][j];
- }
- yTmp[j] = y[j] + stepSize * sum;
+ double sum = a[k-1][0] * yDotK[0][j];
+ for (int l = 1; l < k; ++l) {
+ sum += a[k-1][l] * yDotK[l][j];
+ }
+ yTmp[j] = y[j] + stepSize * sum;
}
computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
- }
+ }
- // estimate the state at the end of the step
- for (int j = 0; j < y0.length; ++j) {
+ // estimate the state at the end of the step
+ for (int j = 0; j < y0.length; ++j) {
double sum = b[0] * yDotK[0][j];
for (int l = 1; l < stages; ++l) {
- sum += b[l] * yDotK[l][j];
+ sum += b[l] * yDotK[l][j];
}
yTmp[j] = y[j] + stepSize * sum;
- }
-
- // discrete events handling
- interpolator.storeTime(stepStart + stepSize);
- if (manager.evaluateStep(interpolator)) {
- final double dt = manager.getEventTime() - stepStart;
- if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
- // we cannot simply truncate the step, reject the current computation
- // and let the loop compute another state with the truncated step.
- // it is so small (much probably exactly 0 due to limited accuracy)
- // that the code above would fail handling it.
- // So we set up an artificial 0 size step by copying states
- interpolator.storeTime(stepStart);
- System.arraycopy(y, 0, yTmp, 0, y0.length);
- stepSize = 0;
- loop = false;
- } else {
- // reject the step to match exactly the next switch time
- stepSize = dt;
- }
- } else {
- loop = false;
- }
-
}
- // the step has been accepted
- final double nextStep = stepStart + stepSize;
+ // discrete events handling
+ interpolator.storeTime(stepStart + stepSize);
System.arraycopy(yTmp, 0, y, 0, y0.length);
- manager.stepAccepted(nextStep, y);
- lastStep = manager.stop();
+ System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
+ stepStart = acceptStep(interpolator, stepHandlers, y, yDotTmp, t);
- // provide the step data to the step handler
- interpolator.storeTime(nextStep);
- for (StepHandler handler : stepHandlers) {
- handler.handleStep(interpolator, lastStep);
- }
- stepStart = nextStep;
+ if (!isLastStep) {
- if (manager.reset(stepStart, y) && ! lastStep) {
- // some events handler has triggered changes that
- // invalidate the derivatives, we need to recompute them
- computeDerivatives(stepStart, y, yDotK[0]);
- }
+ // prepare next step
+ interpolator.storeTime(stepStart);
- // make sure step size is set to default before next step
- stepSize = forward ? step : -step;
+ // stepsize control for next step
+ final double nextT = stepStart + stepSize;
+ final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+ if (nextIsLast) {
+ stepSize = t - stepStart;
+ }
+ }
- }
+ } while (!isLastStep);
final double stopTime = stepStart;
stepStart = Double.NaN;
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java Thu Jan 20 20:55:00 2011
@@ -21,7 +21,6 @@ import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
-import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.MathUserException;
/** This abstract class represents an interpolator over the last step
@@ -44,11 +43,17 @@ import org.apache.commons.math.exception
public abstract class AbstractStepInterpolator
implements StepInterpolator {
- /** previous time */
- protected double previousTime;
+ /** global previous time */
+ private double globalPreviousTime;
- /** current time */
- protected double currentTime;
+ /** global current time */
+ private double globalCurrentTime;
+
+ /** soft previous time */
+ private double softPreviousTime;
+
+ /** soft current time */
+ private double softCurrentTime;
/** current time step */
protected double h;
@@ -87,8 +92,10 @@ public abstract class AbstractStepInterp
* initializing the copy.
*/
protected AbstractStepInterpolator() {
- previousTime = Double.NaN;
- currentTime = Double.NaN;
+ globalPreviousTime = Double.NaN;
+ globalCurrentTime = Double.NaN;
+ softPreviousTime = Double.NaN;
+ softCurrentTime = Double.NaN;
h = Double.NaN;
interpolatedTime = Double.NaN;
currentState = null;
@@ -106,10 +113,12 @@ public abstract class AbstractStepInterp
*/
protected AbstractStepInterpolator(final double[] y, final boolean forward) {
- previousTime = Double.NaN;
- currentTime = Double.NaN;
- h = Double.NaN;
- interpolatedTime = Double.NaN;
+ globalPreviousTime = Double.NaN;
+ globalCurrentTime = Double.NaN;
+ softPreviousTime = Double.NaN;
+ softCurrentTime = Double.NaN;
+ h = Double.NaN;
+ interpolatedTime = Double.NaN;
currentState = y;
interpolatedState = new double[y.length];
@@ -140,10 +149,12 @@ public abstract class AbstractStepInterp
*/
protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
- previousTime = interpolator.previousTime;
- currentTime = interpolator.currentTime;
- h = interpolator.h;
- interpolatedTime = interpolator.interpolatedTime;
+ globalPreviousTime = interpolator.globalPreviousTime;
+ globalCurrentTime = interpolator.globalCurrentTime;
+ softPreviousTime = interpolator.softPreviousTime;
+ softCurrentTime = interpolator.softCurrentTime;
+ h = interpolator.h;
+ interpolatedTime = interpolator.interpolatedTime;
if (interpolator.currentState != null) {
currentState = interpolator.currentState.clone();
@@ -168,10 +179,12 @@ public abstract class AbstractStepInterp
*/
protected void reinitialize(final double[] y, final boolean isForward) {
- previousTime = Double.NaN;
- currentTime = Double.NaN;
- h = Double.NaN;
- interpolatedTime = Double.NaN;
+ globalPreviousTime = Double.NaN;
+ globalCurrentTime = Double.NaN;
+ softPreviousTime = Double.NaN;
+ softCurrentTime = Double.NaN;
+ h = Double.NaN;
+ interpolatedTime = Double.NaN;
currentState = y;
interpolatedState = new double[y.length];
@@ -208,7 +221,9 @@ public abstract class AbstractStepInterp
* interpolator for future calls to {@link #storeTime storeTime}
*/
public void shift() {
- previousTime = currentTime;
+ globalPreviousTime = globalCurrentTime;
+ softPreviousTime = globalPreviousTime;
+ softCurrentTime = globalCurrentTime;
}
/** Store the current step time.
@@ -216,8 +231,9 @@ public abstract class AbstractStepInterp
*/
public void storeTime(final double t) {
- currentTime = t;
- h = currentTime - previousTime;
+ globalCurrentTime = t;
+ softCurrentTime = globalCurrentTime;
+ h = globalCurrentTime - globalPreviousTime;
setInterpolatedTime(t);
// the step is not finalized anymore
@@ -225,14 +241,53 @@ public abstract class AbstractStepInterp
}
- /** {@inheritDoc} */
+ /** Restrict step range to a limited part of the global step.
+ * <p>
+ * This method can be used to restrict a step and make it appear
+ * as if the original step was smaller. Calling this method
+ * <em>only</em> changes the value returned by {@link #getPreviousTime()}
+ * and {@link #getCurrentTime()}, it does not change any
+ * </p>
+ * @param softPreviousTime start of the restricted step
+ * @param softCurrentTime end of the restricted step
+ */
+ public void setSoftBounds(final double softPreviousTime, final double softCurrentTime) {
+ this.softPreviousTime = softPreviousTime;
+ this.softCurrentTime = softCurrentTime;
+ }
+
+ /**
+ * Get the previous global grid point time.
+ * @return previous global grid point time
+ */
+ public double getGlobalPreviousTime() {
+ return globalPreviousTime;
+ }
+
+ /**
+ * Get the current global grid point time.
+ * @return current global grid point time
+ */
+ public double getGlobalCurrentTime() {
+ return globalCurrentTime;
+ }
+
+ /**
+ * Get the previous soft grid point time.
+ * @return previous soft grid point time
+ * @see #setSoftBounds(double, double)
+ */
public double getPreviousTime() {
- return previousTime;
+ return softPreviousTime;
}
- /** {@inheritDoc} */
+ /**
+ * Get the current soft grid point time.
+ * @return current soft grid point time
+ * @see #setSoftBounds(double, double)
+ */
public double getCurrentTime() {
- return currentTime;
+ return softCurrentTime;
}
/** {@inheritDoc} */
@@ -270,7 +325,7 @@ public abstract class AbstractStepInterp
// lazy evaluation of the state
if (dirtyState) {
- final double oneMinusThetaH = currentTime - interpolatedTime;
+ final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
dirtyState = false;
@@ -285,7 +340,7 @@ public abstract class AbstractStepInterp
// lazy evaluation of the state
if (dirtyState) {
- final double oneMinusThetaH = currentTime - interpolatedTime;
+ final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
dirtyState = false;
@@ -376,8 +431,10 @@ public abstract class AbstractStepInterp
} else {
out.writeInt(currentState.length);
}
- out.writeDouble(previousTime);
- out.writeDouble(currentTime);
+ out.writeDouble(globalPreviousTime);
+ out.writeDouble(globalCurrentTime);
+ out.writeDouble(softPreviousTime);
+ out.writeDouble(softCurrentTime);
out.writeDouble(h);
out.writeBoolean(forward);
@@ -396,7 +453,9 @@ public abstract class AbstractStepInterp
try {
finalizeStep();
} catch (MathUserException e) {
- throw MathRuntimeException.createIOException(e);
+ IOException ioe = new IOException(e.getLocalizedMessage());
+ ioe.initCause(e);
+ throw ioe;
}
}
@@ -414,11 +473,13 @@ public abstract class AbstractStepInterp
throws IOException {
final int dimension = in.readInt();
- previousTime = in.readDouble();
- currentTime = in.readDouble();
- h = in.readDouble();
- forward = in.readBoolean();
- dirtyState = true;
+ globalPreviousTime = in.readDouble();
+ globalCurrentTime = in.readDouble();
+ softPreviousTime = in.readDouble();
+ softCurrentTime = in.readDouble();
+ h = in.readDouble();
+ forward = in.readBoolean();
+ dirtyState = true;
if (dimension < 0) {
currentState = null;
Modified: commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml Thu Jan 20 20:55:00 2011
@@ -52,6 +52,10 @@ The <action> type attribute can be add,u
If the output is not quite correct, check for invisible trailing spaces!
-->
<release version="2.2" date="TBD" description="TBD">
+ <action dev="luc" type="fix" issue="MATH-484">
+ separate discrete event detection from adaptive step size handling in ODE solvers,
+ thus improving robustness, maintainability and speed
+ </action>
<action dev="sebb" type="fix" issue="MATH-486">
FastMath toRadian and toDegree don't handle large double numbers well
</action>
Modified: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblem4.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblem4.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblem4.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblem4.java Thu Jan 20 20:55:00 2011
@@ -80,6 +80,20 @@ public TestProblem4 copy() {
return new EventHandler[] { new Bounce(), new Stop() };
}
+ /**
+ * Get the theoretical events times.
+ * @return theoretical events times
+ */
+ public double[] getTheoreticalEventsTimes() {
+ return new double[] {
+ 1 * FastMath.PI - a,
+ 2 * FastMath.PI - a,
+ 3 * FastMath.PI - a,
+ 4 * FastMath.PI - a,
+ 12.0
+ };
+ }
+
@Override
public void doComputeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = y[1];
Modified: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java Thu Jan 20 20:55:00 2011
@@ -159,6 +159,14 @@ public abstract class TestProblemAbstrac
}
/**
+ * Get the theoretical events times.
+ * @return theoretical events times
+ */
+ public double[] getTheoreticalEventsTimes() {
+ return new double[0];
+ }
+
+ /**
* Get the number of calls.
* @return nuber of calls
*/
Modified: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java Thu Jan 20 20:55:00 2011
@@ -80,7 +80,13 @@ public class TestProblemHandler
// multistep integrators do not handle the first steps themselves
// so we have to make sure the integrator we look at has really started its work
if (!Double.isNaN(expectedStepStart)) {
- maxTimeError = FastMath.max(maxTimeError, FastMath.abs(start - expectedStepStart));
+ // the step should either start at the end of the integrator step
+ // or at an event if the step is split into several substeps
+ double stepError = FastMath.max(maxTimeError, FastMath.abs(start - expectedStepStart));
+ for (double eventTime : problem.getTheoreticalEventsTimes()) {
+ stepError = FastMath.min(stepError, FastMath.abs(start - eventTime));
+ }
+ maxTimeError = FastMath.max(maxTimeError, stepError);
}
expectedStepStart = start + integrator.getCurrentSignedStepsize();
}
Modified: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java?rev=1061507&r1=1061506&r2=1061507&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java Thu Jan 20 20:55:00 2011
@@ -111,8 +111,8 @@ public class FirstOrderIntegratorWithJac
extInt.setMaxEvaluations(5000);
extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
Assert.assertEquals(5000, extInt.getMaxEvaluations());
- Assert.assertTrue(extInt.getEvaluations() > 1500);
- Assert.assertTrue(extInt.getEvaluations() < 2100);
+ Assert.assertTrue(extInt.getEvaluations() > 1400);
+ Assert.assertTrue(extInt.getEvaluations() < 2000);
Assert.assertEquals(4 * integ.getEvaluations(), extInt.getEvaluations());
residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());