Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

separate discrete event detection from adaptive step size handling in ODE solvers, 

thus improving robustness, maintainability and speed
JIRA: MATH-484


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1061508 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2011-01-20 20:57:11 +00:00
parent 168139c161
commit 73a8e5d165
36 changed files with 639 additions and 887 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

196
src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
View File

@ -20,15 +20,21 @@ 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 +52,17 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
/** 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 +84,8 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
stepHandlers = new ArrayList<StepHandler>();
stepStart = Double.NaN;
stepSize = Double.NaN;
eventsHandlersManager = new CombinedEventsManager();
eventsStates = new ArrayList<EventState>();
statesInitialized = false;
setMaxEvaluations(-1);
resetEvaluations();
}
@ -101,22 +117,25 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
}
/** {@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 +204,109 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
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
@ -216,60 +338,4 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
}
/** Add an event handler for end time checking.
* <p>This method can be used to simplify handling of integration end time.
* It leverages the nominal stop condition with the exceptional stop
* conditions.</p>
* @param startTime integration start time
* @param endTime desired end time
* @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
*/
protected CombinedEventsManager addEndTimeChecker(final double startTime,
final double endTime,
final CombinedEventsManager manager) {
CombinedEventsManager newManager = new CombinedEventsManager();
for (final EventState state : manager.getEventsStates()) {
newManager.addEventHandler(state.getEventHandler(),
state.getMaxCheckInterval(),
state.getConvergence(),
state.getMaxIterationCount());
}
newManager.addEventHandler(new EndTimeChecker(endTime),
Double.POSITIVE_INFINITY,
FastMath.ulp(FastMath.max(FastMath.abs(startTime), FastMath.abs(endTime))),
100);
return newManager;
}
/** Specialized event handler to stop integration. */
private static class EndTimeChecker implements EventHandler {
/** Desired end time. */
private final double endTime;
/** Build an instance.
* @param endTime desired time
*/
public EndTimeChecker(final double endTime) {
this.endTime = endTime;
}
/** {@inheritDoc} */
public int eventOccurred(double t, double[] y, boolean increasing) {
return STOP;
}
/** {@inheritDoc} */
public double g(double t, double[] y) {
return t - endTime;
}
/** {@inheritDoc} */
public void resetState(double t, double[] y) {
}
}
}

6
src/main/java/org/apache/commons/math/ode/ODEIntegrator.java
View File

@ -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

244
src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java
View File

@ -1,244 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.events;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.List;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.sampling.StepInterpolator;
/** This class manages several {@link EventHandler event handlers} during integration.
*
* @see EventHandler
* @see EventState
* @version $Revision$ $Date$
* @since 1.2
*/
public class CombinedEventsManager {
/** Events states. */
private final List<EventState> states;
/** First active event. */
private EventState first;
/** Initialization indicator. */
private boolean initialized;
/** Simple constructor.
* Create an empty manager
*/
public CombinedEventsManager() {
states = new ArrayList<EventState>();
first = null;
initialized = false;
}
/** Add an events handler.
* @param handler event handler
* @param maxCheckInterval maximal time interval between events
* checks (this interval prevents missing sign changes in
* case the integration steps becomes very large)
* @param convergence convergence threshold in the event time search
* @param maxIterationCount upper limit of the iteration count in
* the event time search
* @see #getEventsHandlers()
* @see #clearEventsHandlers()
*/
public void addEventHandler(final EventHandler handler, final double maxCheckInterval,
final double convergence, final int maxIterationCount) {
states.add(new EventState(handler, maxCheckInterval,
convergence, maxIterationCount));
}
/** Get all the events handlers that have been added to the manager.
* @return an unmodifiable collection of the added event handlers
* @see #addEventHandler(EventHandler, double, double, int)
* @see #clearEventsHandlers()
* @see #getEventsStates()
*/
public Collection<EventHandler> getEventsHandlers() {
final List<EventHandler> list = new ArrayList<EventHandler>();
for (EventState state : states) {
list.add(state.getEventHandler());
}
return Collections.unmodifiableCollection(list);
}
/** Remove all the events handlers that have been added to the manager.
* @see #addEventHandler(EventHandler, double, double, int)
* @see #getEventsHandlers()
*/
public void clearEventsHandlers() {
states.clear();
}
/** Get all the events state wrapping the handlers that have been added to the manager.
* @return a collection of the events states
* @see #getEventsHandlers()
*/
public Collection<EventState> getEventsStates() {
return states;
}
/** Check if the manager does not manage any event handlers.
* @return true if manager is empty
*/
public boolean isEmpty() {
return states.isEmpty();
}
/** Evaluate the impact of the proposed step on all managed
* event handlers.
* @param interpolator step interpolator for the proposed step
* @return true if at least one event handler triggers an event
* before the end of the proposed step (this implies the step should
* be rejected)
* @exception MathUserException if the interpolator fails to
* compute the function somewhere within the step
* @exception IntegratorException if an event cannot be located
*/
public boolean evaluateStep(final StepInterpolator interpolator)
throws MathUserException, IntegratorException {
try {
first = null;
if (states.isEmpty()) {
// there is nothing to do, return now to avoid setting the
// interpolator time (and hence avoid unneeded calls to the
// user function due to interpolator finalization)
return false;
}
if (! initialized) {
// initialize the events states
for (EventState state : states) {
state.reinitializeBegin(interpolator);
}
initialized = true;
}
// check events occurrence
for (EventState state : states) {
if (state.evaluateStep(interpolator)) {
if (first == null) {
first = state;
} else {
if (interpolator.isForward()) {
if (state.getEventTime() < first.getEventTime()) {
first = state;
}
} else {
if (state.getEventTime() > first.getEventTime()) {
first = state;
}
}
}
}
}
return first != null;
} catch (EventException se) {
throw new IntegratorException(se);
} catch (ConvergenceException ce) {
throw new IntegratorException(ce);
}
}
/** Get the occurrence time of the first event triggered in the
* last evaluated step.
* @return occurrence time of the first event triggered in the last
* evaluated step, or </code>Double.NaN</code> if no event is
* triggered
*/
public double getEventTime() {
return (first == null) ? Double.NaN : first.getEventTime();
}
/** Inform the event handlers that the step has been accepted
* by the integrator.
* @param t value of the independent <i>time</i> variable at the
* end of the step
* @param y array containing the current value of the state vector
* at the end of the step
* @exception IntegratorException if the value of one of the
* events states cannot be evaluated
*/
public void stepAccepted(final double t, final double[] y)
throws IntegratorException {
try {
for (EventState state : states) {
state.stepAccepted(t, y);
}
} catch (EventException se) {
throw new IntegratorException(se);
}
}
/** Check if the integration should be stopped at the end of the
* current step.
* @return true if the integration should be stopped
*/
public boolean stop() {
for (EventState state : states) {
if (state.stop()) {
return true;
}
}
return false;
}
/** Let the event handlers reset the state if they want.
* @param t value of the independent <i>time</i> variable at the
* beginning of the next step
* @param y array were to put the desired state vector at the beginning
* of the next step
* @return true if the integrator should reset the derivatives too
* @exception IntegratorException if one of the events states
* that should reset the state fails to do it
*/
public boolean reset(final double t, final double[] y)
throws IntegratorException {
try {
boolean resetDerivatives = false;
for (EventState state : states) {
if (state.reset(t, y)) {
resetDerivatives = true;
}
}
return resetDerivatives;
} catch (EventException se) {
throw new IntegratorException(se);
}
}
}

109
src/main/java/org/apache/commons/math/ode/events/EventState.java
View File

@ -18,10 +18,10 @@
package org.apache.commons.math.ode.events;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.exception.MathUserException;
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;
@ -33,13 +33,12 @@ import org.apache.commons.math.util.FastMath;
* of one handler during one integration step, with references to the
* state at the end of the preceding step. This information is used to
* decide if the handler should trigger an event or not during the
* proposed step (and hence the step should be reduced to ensure the
* event occurs at a bound rather than inside the step).</p>
* proposed step.</p>
*
* @version $Revision$ $Date$
* @since 1.2
*/
public class EventState {
public class EventState implements Comparable<EventState> {
/** Event handler. */
private final EventHandler handler;
@ -187,8 +186,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 +200,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,26 +222,6 @@ public class EventState {
if (g0Positive ^ (gb >= 0)) {
// there is a sign change: an event is expected during this step
if (ga * gb > 0) {
// this is a corner case:
// - there was an event near ta,
// - there is another event between ta and tb
// - when ta was computed, convergence was reached on the "wrong side" of the interval
// this implies that the real sign of ga is the same as gb, so we need to slightly
// shift ta to make sure ga and gb get opposite signs and the solver won't complain
// about bracketing
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());
}
if (ga * gb > 0) {
// this should never happen
throw MathRuntimeException.createInternalError(null);
}
}
// variation direction, with respect to the integration direction
increasing = gb >= ga;
@ -253,28 +236,45 @@ public class EventState {
}
};
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
// - when ta was computed, convergence was reached on the "wrong side" of the interval
// this implies that the real sign of ga is the same as gb, so we need to slightly
// shift ta to make sure ga and gb get opposite signs and the solver won't complain
// about bracketing
final double epsilon = (forward ? 0.25 : -0.25) * convergence;
for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
ta += epsilon;
ga = f.value(ta);
}
if (ga * gb > 0) {
// this should never happen
throw new MathInternalError();
}
}
final double root = (ta <= tb) ?
solver.solve(maxIterationCount, f, ta, tb) :
solver.solve(maxIterationCount, f, tb, ta);
if ((FastMath.abs(root - ta) <= convergence) &&
(FastMath.abs(root - previousEventTime) <= convergence)) {
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 {
@ -323,7 +323,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;
@ -354,7 +354,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;
}
@ -369,4 +369,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);
}
}

138
src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
View File

@ -21,7 +21,6 @@ import org.apache.commons.math.exception.MathUserException;
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 extends AdamsIntegrator {
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 extends AdamsIntegrator {
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 extends AdamsIntegrator {
}
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();
// 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);
// provide the step data to the step handler
for (StepHandler handler : stepHandlers) {
interpolator.setInterpolatedTime(nextStep);
handler.handleStep(interpolator, lastStep);
// 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];
}
stepStart = nextStep;
final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
interpolator.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
if (!lastStep && manager.reset(stepStart, y)) {
// discrete events handling
interpolator.storeTime(stepEnd);
stepStart = acceptStep(interpolator, stepHandlers, y, yDot, t);
scaled = predictedScaled;
nordsieck = nordsieckTmp;
interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
// 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 (!isLastStep) {
}
// prepare next step
interpolator.storeTime(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);
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 extends AdamsIntegrator {
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;
}

128
src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
View File

@ -24,7 +24,6 @@ import org.apache.commons.math.linear.Array2DRowRealMatrix;
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 extends AdamsIntegrator {
}
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 extends AdamsIntegrator {
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 extends AdamsIntegrator {
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;
// 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];
}
updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
// discrete events handling
System.arraycopy(yTmp, 0, y, 0, n);
interpolator.storeTime(nextStep);
manager.stepAccepted(nextStep, y);
lastStep = manager.stop();
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;
// provide the step data to the step handler
for (StepHandler handler : stepHandlers) {
interpolator.setInterpolatedTime(nextStep);
handler.handleStep(interpolator, lastStep);
}
stepStart = nextStep;
if (!isLastStep) {
if (!lastStep && manager.reset(stepStart, y)) {
// prepare next step
interpolator.storeTime(stepStart);
// 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 (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 extends AdamsIntegrator {
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;

12
src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
View File

@ -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);

115
src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
View File

@ -20,7 +20,6 @@ package org.apache.commons.math.ode.nonstiff;
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 EmbeddedRungeKuttaIntegrator
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 EmbeddedRungeKuttaIntegrator
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 EmbeddedRungeKuttaIntegrator
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 EmbeddedRungeKuttaIntegrator
// 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) {
// prepare next step
interpolator.storeTime(stepStart);
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 (fsal) {
// save the last evaluation for the next step
System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
}
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 (! 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();

95
src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
View File

@ -170,7 +170,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
@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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
// 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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
}
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
break;
default :
if ((firstTime || lastStep) && (error <= 1.0)) {
if ((firstTime || isLastStep) && (error <= 1.0)) {
loop = false;
}
break;
@ -816,6 +810,11 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
}
}
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
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 GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
}
}
// 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;
// Discrete events handling
interpolator.storeTime(stepStart + stepSize);
stepStart = acceptStep(interpolator, stepHandlers, y1, yDot1, t);
// prepare next step
interpolator.storeTime(stepStart);
System.arraycopy(y1, 0, y, 0, y0.length);
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;
}
System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
firstStepAlreadyComputed = true;
int optimalIter;
if (k == 1) {
@ -987,15 +952,17 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
firstTime = false;
if (reject) {
lastStep = false;
isLastStep = false;
previousRejected = true;
} else {
previousRejected = false;
}
}
} while (!isLastStep);
return stepStart;
final double stopTime = stepStart;
resetInternalState();
return stopTime;
}

4
src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
View File

@ -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;

97
src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
View File

@ -22,7 +22,6 @@ import org.apache.commons.math.exception.MathUserException;
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 RungeKuttaIntegrator extends AbstractIntegrator {
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 RungeKuttaIntegrator extends AbstractIntegrator {
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);
// stepsize control for next step
final double nextT = stepStart + stepSize;
final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
if (nextIsLast) {
stepSize = t - stepStart;
}
}
// make sure step size is set to default before next step
stepSize = forward ? step : -step;
}
} while (!isLastStep);
final double stopTime = stepStart;
stepStart = Double.NaN;

135
src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
View File

@ -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.MathUserException;
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 AbstractStepInterpolator
* 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 AbstractStepInterpolator
*/
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 AbstractStepInterpolator
*/
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 AbstractStepInterpolator
*/
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 AbstractStepInterpolator
* 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 AbstractStepInterpolator
*/
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 AbstractStepInterpolator
}
/** {@inheritDoc} */
public double getPreviousTime() {
return previousTime;
/** 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;
}
/** {@inheritDoc} */
/**
* 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 softPreviousTime;
}
/**
* 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 AbstractStepInterpolator
// 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 AbstractStepInterpolator
// 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 AbstractStepInterpolator
} 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 AbstractStepInterpolator
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 AbstractStepInterpolator
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;

4
src/site/xdoc/changes.xml
View File

@ -169,6 +169,10 @@ The <action> type attribute can be add,update,fix,remove.
</action>
</release>
<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="luc" type="fix" issue="MATH-467">
Fixed an awkward statement that triggered a false positive warning.
</action>

14
src/test/java/org/apache/commons/math/ode/TestProblem4.java
View File

@ -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];

8
src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java
View File

@ -158,6 +158,14 @@ public abstract class TestProblemAbstract
return new EventHandler[0];
}
/**
* 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

8
src/test/java/org/apache/commons/math/ode/TestProblemHandler.java
View File

@ -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();
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
View File

@ -146,9 +146,9 @@ public class AdamsBashforthIntegratorTest {
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
if (nSteps < 4) {
assertTrue(integ.getEvaluations() > 160);
assertTrue(integ.getEvaluations() > 150);
} else {
assertTrue(integ.getEvaluations() < 80);
assertTrue(integ.getEvaluations() < 70);
}
}

15
src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
View File

@ -134,7 +134,8 @@ public class ClassicalRungeKuttaIntegratorTest
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
double previousError = Double.NaN;
double previousValueError = Double.NaN;
double previousTimeError = Double.NaN;
for (int i = 4; i < 10; ++i) {
TestProblemAbstract pb = problems[k].copy();
@ -157,10 +158,16 @@ public class ClassicalRungeKuttaIntegratorTest
double error = handler.getMaximalValueError();
if (i > 4) {
assertTrue(error < FastMath.abs(previousError));
assertTrue(error < FastMath.abs(previousValueError));
}
previousError = error;
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
previousValueError = error;
double timeError = handler.getMaximalTimeError();
if (i > 4) {
assertTrue(timeError <= FastMath.abs(previousTimeError));
}
previousTimeError = timeError;
integ.clearEventHandlers();
assertEquals(0, integ.getEventHandlers().size());
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
View File

@ -64,8 +64,8 @@ public class ClassicalRungeKuttaStepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 700000);
assertTrue(bos.size () < 701000);
assertTrue(bos.size () > 753000);
assertTrue(bos.size () < 754000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

7
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java
View File

@ -217,9 +217,10 @@ public class DormandPrince54IntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
EventHandler[] functions = pb.getEventsHandlers();
double convergence = 1.0e-8 * maxStep;
for (int l = 0; l < functions.length; ++l) {
integ.addEventHandler(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
Double.POSITIVE_INFINITY, convergence, 1000);
}
assertEquals(functions.length, integ.getEventHandlers().size());
integ.integrate(pb,
@ -227,8 +228,8 @@ public class DormandPrince54IntegratorTest
pb.getFinalTime(), new double[pb.getDimension()]);
assertTrue(handler.getMaximalValueError() < 5.0e-6);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
assertEquals(0, handler.getMaximalTimeError(), convergence);
assertEquals(12.0, handler.getLastTime(), convergence);
integ.clearEventHandlers();
assertEquals(0, integ.getEventHandlers().size());

4
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
View File

@ -77,8 +77,8 @@ public class DormandPrince54StepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 119500);
assertTrue(bos.size () < 120500);
assertTrue(bos.size () > 126000);
assertTrue(bos.size () < 127000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

9
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
View File

@ -205,7 +205,7 @@ public class DormandPrince853IntegratorTest
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
assertTrue(handler.getLastError() < 8.1e-8);
assertTrue(handler.getLastError() < 1.1e-7);
assertTrue(handler.getMaximalValueError() < 1.1e-7);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals("Dormand-Prince 8 (5, 3)", integ.getName());
@ -226,9 +226,10 @@ public class DormandPrince853IntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
EventHandler[] functions = pb.getEventsHandlers();
double convergence = 1.0e-8 * maxStep;
for (int l = 0; l < functions.length; ++l) {
integ.addEventHandler(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
Double.POSITIVE_INFINITY, convergence, 1000);
}
assertEquals(functions.length, integ.getEventHandlers().size());
integ.integrate(pb,
@ -236,8 +237,8 @@ public class DormandPrince853IntegratorTest
pb.getFinalTime(), new double[pb.getDimension()]);
assertEquals(0, handler.getMaximalValueError(), 1.1e-7);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
assertEquals(0, handler.getMaximalTimeError(), convergence);
assertEquals(12.0, handler.getLastTime(), convergence);
integ.clearEventHandlers();
assertEquals(0, integ.getEventHandlers().size());

4
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
View File

@ -77,8 +77,8 @@ public class DormandPrince853StepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 86000);
assertTrue(bos.size () < 87000);
assertTrue(bos.size () > 88000);
assertTrue(bos.size () < 89000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

21
src/test/java/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java
View File

@ -60,12 +60,12 @@ public class EulerIntegratorTest
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
double previousError = Double.NaN;
for (int i = 4; i < 10; ++i) {
double previousValueError = Double.NaN;
double previousTimeError = Double.NaN;
for (int i = 4; i < 8; ++i) {
TestProblemAbstract pb = problems[k].copy();
double step = (pb.getFinalTime() - pb.getInitialTime())
* FastMath.pow(2.0, -i);
double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i);
FirstOrderIntegrator integ = new EulerIntegrator(step);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
@ -81,12 +81,17 @@ public class EulerIntegratorTest
assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
}
double error = handler.getMaximalValueError();
double valueError = handler.getMaximalValueError();
if (i > 4) {
assertTrue(error < FastMath.abs(previousError));
assertTrue(valueError < FastMath.abs(previousValueError));
}
previousError = error;
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
previousValueError = valueError;
double timeError = handler.getMaximalTimeError();
if (i > 4) {
assertTrue(timeError <= FastMath.abs(previousTimeError));
}
previousTimeError = timeError;
}

19
src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
View File

@ -61,12 +61,12 @@ public class GillIntegratorTest
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
double previousError = Double.NaN;
double previousValueError = Double.NaN;
double previousTimeError = Double.NaN;
for (int i = 5; i < 10; ++i) {
TestProblemAbstract pb = problems[k].copy();
double step = (pb.getFinalTime() - pb.getInitialTime())
* FastMath.pow(2.0, -i);
double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i);
FirstOrderIntegrator integ = new GillIntegrator(step);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
@ -82,12 +82,17 @@ public class GillIntegratorTest
assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
}
double error = handler.getMaximalValueError();
double valueError = handler.getMaximalValueError();
if (i > 5) {
assertTrue(error < FastMath.abs(previousError));
assertTrue(valueError < FastMath.abs(previousValueError));
}
previousError = error;
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
previousValueError = valueError;
double timeError = handler.getMaximalTimeError();
if (i > 5) {
assertTrue(timeError <= FastMath.abs(previousTimeError));
}
previousTimeError = timeError;
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java
View File

@ -65,8 +65,8 @@ public class GillStepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 700000);
assertTrue(bos.size () < 701000);
assertTrue(bos.size () > 753000);
assertTrue(bos.size () < 754000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

12
src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java
View File

@ -113,8 +113,8 @@ public class GraggBulirschStoerIntegratorTest
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
assertTrue(handler.getLastError() < 9.0e-10);
assertTrue(handler.getMaximalValueError() < 9.0e-10);
assertTrue(handler.getLastError() < 7.5e-9);
assertTrue(handler.getMaximalValueError() < 8.1e-9);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals("Gragg-Bulirsch-Stoer", integ.getName());
}
@ -210,9 +210,9 @@ public class GraggBulirschStoerIntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
EventHandler[] functions = pb.getEventsHandlers();
double convergence = 1.0e-8 * maxStep;
for (int l = 0; l < functions.length; ++l) {
integ.addEventHandler(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
integ.addEventHandler(functions[l], Double.POSITIVE_INFINITY, convergence, 1000);
}
assertEquals(functions.length, integ.getEventHandlers().size());
integ.integrate(pb,
@ -220,8 +220,8 @@ public class GraggBulirschStoerIntegratorTest
pb.getFinalTime(), new double[pb.getDimension()]);
assertTrue(handler.getMaximalValueError() < 5.0e-8);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
assertEquals(0, handler.getMaximalTimeError(), convergence);
assertEquals(12.0, handler.getLastTime(), convergence);
integ.clearEventHandlers();
assertEquals(0, integ.getEventHandlers().size());

4
src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
View File

@ -79,8 +79,8 @@ public class GraggBulirschStoerStepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 33000);
assertTrue(bos.size () < 34000);
assertTrue(bos.size () > 34000);
assertTrue(bos.size () < 35000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

48
src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
View File

@ -32,8 +32,6 @@ import org.apache.commons.math.ode.TestProblem5;
import org.apache.commons.math.ode.TestProblemHandler;
import org.apache.commons.math.ode.events.EventException;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;
import org.apache.commons.math.util.FastMath;
public class HighamHall54IntegratorTest
@ -176,9 +174,10 @@ public class HighamHall54IntegratorTest
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
EventHandler[] functions = pb.getEventsHandlers();
double convergence = 1.0e-8 * maxStep;
for (int l = 0; l < functions.length; ++l) {
integ.addEventHandler(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
Double.POSITIVE_INFINITY, convergence, 1000);
}
assertEquals(functions.length, integ.getEventHandlers().size());
integ.integrate(pb,
@ -186,8 +185,8 @@ public class HighamHall54IntegratorTest
pb.getFinalTime(), new double[pb.getDimension()]);
assertTrue(handler.getMaximalValueError() < 1.0e-7);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
assertEquals(0, handler.getMaximalTimeError(), convergence);
assertEquals(12.0, handler.getLastTime(), convergence);
integ.clearEventHandlers();
assertEquals(0, integ.getEventHandlers().size());
@ -343,46 +342,13 @@ public class HighamHall54IntegratorTest
FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
vecAbsoluteTolerance,
vecRelativeTolerance);
integ.addStepHandler(new KeplerHandler(pb));
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
assertEquals(0.0, handler.getMaximalValueError(), 1.5e-4);
assertEquals("Higham-Hall 5(4)", integ.getName());
}
private static class KeplerHandler implements StepHandler {
public KeplerHandler(TestProblem3 pb) {
this.pb = pb;
nbSteps = 0;
maxError = 0;
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
nbSteps = 0;
maxError = 0;
}
public void handleStep(StepInterpolator interpolator,
boolean isLast) throws MathUserException {
++nbSteps;
double[] interpolatedY = interpolator.getInterpolatedState();
double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
double error = dx * dx + dy * dy;
if (error > maxError) {
maxError = error;
}
if (isLast) {
assertTrue(maxError < 4.2e-11);
assertTrue(nbSteps < 670);
}
}
private TestProblem3 pb;
private int nbSteps;
private double maxError;
}
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java
View File

@ -77,8 +77,8 @@ public class HighamHall54StepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 158000);
assertTrue(bos.size () < 159000);
assertTrue(bos.size () > 167000);
assertTrue(bos.size () < 168000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

19
src/test/java/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java
View File

@ -60,12 +60,12 @@ public class MidpointIntegratorTest
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
double previousError = Double.NaN;
double previousValueError = Double.NaN;
double previousTimeError = Double.NaN;
for (int i = 4; i < 10; ++i) {
TestProblemAbstract pb = problems[k].copy();
double step = (pb.getFinalTime() - pb.getInitialTime())
* FastMath.pow(2.0, -i);
double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i);
FirstOrderIntegrator integ = new MidpointIntegrator(step);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
@ -81,12 +81,17 @@ public class MidpointIntegratorTest
assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
}
double error = handler.getMaximalValueError();
double valueError = handler.getMaximalValueError();
if (i > 4) {
assertTrue(error < FastMath.abs(previousError));
assertTrue(valueError < FastMath.abs(previousValueError));
}
previousError = error;
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
previousValueError = valueError;
double timeError = handler.getMaximalTimeError();
if (i > 4) {
assertTrue(timeError <= FastMath.abs(previousTimeError));
}
previousTimeError = timeError;
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java
View File

@ -65,8 +65,8 @@ public class MidpointStepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 98000);
assertTrue(bos.size () < 99000);
assertTrue(bos.size () > 114000);
assertTrue(bos.size () < 115000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

17
src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
View File

@ -61,12 +61,12 @@ public class ThreeEighthesIntegratorTest
TestProblemAbstract[] problems = TestProblemFactory.getProblems();
for (int k = 0; k < problems.length; ++k) {
double previousError = Double.NaN;
double previousValueError = Double.NaN;
double previousTimeError = Double.NaN;
for (int i = 4; i < 10; ++i) {
TestProblemAbstract pb = problems[k].copy();
double step = (pb.getFinalTime() - pb.getInitialTime())
* FastMath.pow(2.0, -i);
double step = (pb.getFinalTime() - pb.getInitialTime()) * FastMath.pow(2.0, -i);
FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
@ -84,10 +84,15 @@ public class ThreeEighthesIntegratorTest
double error = handler.getMaximalValueError();
if (i > 4) {
assertTrue(error < FastMath.abs(previousError));
assertTrue(error < FastMath.abs(previousValueError));
}
previousError = error;
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
previousValueError = error;
double timeError = handler.getMaximalTimeError();
if (i > 4) {
assertTrue(timeError <= FastMath.abs(previousTimeError));
}
previousTimeError = timeError;
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
View File

@ -64,8 +64,8 @@ public class ThreeEighthesStepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 700000);
assertTrue(bos.size () < 701000);
assertTrue(bos.size () > 753000);
assertTrue(bos.size () < 754000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

4
src/test/java/org/apache/commons/math/ode/sampling/DummyStepInterpolatorTest.java
View File

@ -91,8 +91,8 @@ public class DummyStepInterpolatorTest {
ObjectOutputStream oos = new ObjectOutputStream(bos);
oos.writeObject(interpolator);
assertTrue(bos.size () > 150);
assertTrue(bos.size () < 250);
assertTrue(bos.size () > 200);
assertTrue(bos.size () < 300);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);

6
src/test/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
View File

@ -41,7 +41,7 @@ public class NordsieckStepInterpolatorTest {
throws MathUserException, IntegratorException {
TestProblem3 pb = new TestProblem3();
AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(4, 0.0, 1.0, 1.0e-10, 1.0e-10);
StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 7e-10);
StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 5e-9);
}
@Test
@ -62,8 +62,8 @@ public class NordsieckStepInterpolatorTest {
oos.writeObject(handler);
}
assertTrue(bos.size () > 20000);
assertTrue(bos.size () < 25000);
assertTrue(bos.size () > 25500);
assertTrue(bos.size () < 26500);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);