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

Fixed inconsistency preventing use of secondary states in ODE events. 

JIRA: MATH-965

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1465654 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2013-04-08 14:41:32 +00:00
parent 2cc9931e2b
commit fda2af9247
4 changed files with 129 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/changes/changes.xml
View File

@ -51,6 +51,10 @@ If the output is not quite correct, check for invisible trailing spaces!
</properties>
<body>
<release version="x.y" date="TBD" description="TBD">
<action dev="luc" type="fix" issue="MATH-965" >
Fixed inconsistent dimensions preventing use of secondary states
in ODE events.
</action>
</release>
<release version="3.2" date="2013-04-06" description="
This is a minor release: It combines bug fixes and new features.

13
src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java
View File

@ -188,6 +188,7 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
evaluations.resetCount();
for (final EventState state : eventsStates) {
state.setExpandable(expandable);
state.getEventHandler().init(t0, y0, t);
}
@ -356,7 +357,6 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
// get state at event time
interpolator.setInterpolatedTime(eventT);
final double[] eventYPrimary = interpolator.getInterpolatedState().clone();
final double[] eventYComplete = new double[y.length];
expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
eventYComplete);
@ -368,7 +368,7 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
// advance all event states to current time
for (final EventState state : eventsStates) {
state.stepAccepted(eventT, eventYPrimary);
state.stepAccepted(eventT, eventYComplete);
isLastStep = isLastStep || state.stop();
}
@ -412,7 +412,14 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
// last part of the step, after the last event
interpolator.setInterpolatedTime(currentT);
final double[] currentY = interpolator.getInterpolatedState();
final double[] currentY = new double[y.length];
expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
currentY);
int index = 0;
for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
currentY);
}
for (final EventState state : eventsStates) {
state.stepAccepted(currentT, currentY);
isLastStep = isLastStep || state.stop();

41
src/main/java/org/apache/commons/math3/ode/events/EventState.java
View File

@ -25,6 +25,8 @@ import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.ode.EquationsMapper;
import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.sampling.StepInterpolator;
import org.apache.commons.math3.util.FastMath;
@ -55,6 +57,9 @@ public class EventState {
/** Upper limit in the iteration count for event localization. */
private final int maxIterationCount;
/** Equation being integrated. */
private ExpandableStatefulODE expandable;
/** Time at the beginning of the step. */
private double t0;
@ -107,6 +112,7 @@ public class EventState {
this.solver = solver;
// some dummy values ...
expandable = null;
t0 = Double.NaN;
g0 = Double.NaN;
g0Positive = true;
@ -125,6 +131,13 @@ public class EventState {
return handler;
}
/** Set the equation.
* @param expandable equation being integrated
*/
public void setExpandable(final ExpandableStatefulODE expandable) {
this.expandable = expandable;
}
/** Get the maximal time interval between events handler checks.
* @return maximal time interval between events handler checks
*/
@ -156,7 +169,7 @@ public class EventState {
t0 = interpolator.getPreviousTime();
interpolator.setInterpolatedTime(t0);
g0 = handler.g(t0, interpolator.getInterpolatedState());
g0 = handler.g(t0, getCompleteState(interpolator));
if (g0 == 0) {
// excerpt from MATH-421 issue:
// If an ODE solver is setup with an EventHandler that return STOP
@ -175,12 +188,32 @@ public class EventState {
FastMath.abs(solver.getRelativeAccuracy() * t0));
final double tStart = t0 + 0.5 * epsilon;
interpolator.setInterpolatedTime(tStart);
g0 = handler.g(tStart, interpolator.getInterpolatedState());
g0 = handler.g(tStart, getCompleteState(interpolator));
}
g0Positive = g0 >= 0;
}
/** Get the complete state (primary and secondary).
* @param interpolator interpolator to use
* @return complete state
*/
private double[] getCompleteState(final StepInterpolator interpolator) {
final double[] complete = new double[expandable.getTotalDimension()];
expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
complete);
int index = 0;
for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
complete);
}
return complete;
}
/** 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
@ -207,7 +240,7 @@ public class EventState {
public double value(final double t) throws LocalMaxCountExceededException {
try {
interpolator.setInterpolatedTime(t);
return handler.g(t, interpolator.getInterpolatedState());
return handler.g(t, getCompleteState(interpolator));
} catch (MaxCountExceededException mcee) {
throw new LocalMaxCountExceededException(mcee);
}
@ -221,7 +254,7 @@ public class EventState {
// evaluate handler value at the end of the substep
final double tb = t0 + (i + 1) * h;
interpolator.setInterpolatedTime(tb);
final double gb = handler.g(tb, interpolator.getInterpolatedState());
final double gb = handler.g(tb, getCompleteState(interpolator));
// check events occurrence
if (g0Positive ^ (gb >= 0)) {

78
src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java
View File

@ -23,7 +23,9 @@ import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.SecondaryEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math3.ode.sampling.DummyStepInterpolator;
@ -56,6 +58,13 @@ public class EventStateTest {
EventState es = new EventState(closeEventsGenerator, 1.5 * gap,
tolerance, 100,
new BrentSolver(tolerance));
es.setExpandable(new ExpandableStatefulODE(new FirstOrderDifferentialEquations() {
public int getDimension() {
return 0;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
}
}));
AbstractStepInterpolator interpolator =
new DummyStepInterpolator(new double[0], new double[0], true);
@ -145,5 +154,74 @@ public class EventStateTest {
}
// Jira: MATH-965
@Test
public void testIssue965()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
ExpandableStatefulODE equation =
new ExpandableStatefulODE(new FirstOrderDifferentialEquations() {
public int getDimension() {
return 1;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = 2.0;
}
});
equation.setTime(0.0);
equation.setPrimaryState(new double[1]);
equation.addSecondaryEquations(new SecondaryEquations() {
public int getDimension() {
return 1;
}
public void computeDerivatives(double t, double[] primary,
double[] primaryDot, double[] secondary,
double[] secondaryDot) {
secondaryDot[0] = -3.0;
}
});
int index = equation.getSecondaryMappers()[0].getFirstIndex();
DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
integrator.addEventHandler(new SecondaryStateEvent(index, -3.0), 0.1, 1.0e-9, 1000);
integrator.setInitialStepSize(3.0);
integrator.integrate(equation, 30.0);
Assert.assertEquals( 1.0, equation.getTime(), 1.0e-10);
Assert.assertEquals( 2.0, equation.getPrimaryState()[0], 1.0e-10);
Assert.assertEquals(-3.0, equation.getSecondaryState(0)[0], 1.0e-10);
}
private static class SecondaryStateEvent implements EventHandler {
private int index;
private final double target;
public SecondaryStateEvent(final int index, final double target) {
this.index = index;
this.target = target;
}
public void init(double t0, double[] y0, double t) {
}
public double g(double t, double[] y) {
return y[index] - target;
}
public Action eventOccurred(double t, double[] y, boolean increasing) {
return Action.STOP;
}
public void resetState(double t, double[] y) {
}
}
}