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

allow switching functions to trigger FunctionEvaluationException 

let the low level errors in switching functions find their way to upper level

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@591674 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2007-11-03 20:23:39 +00:00
parent af62f7927f
commit 62e00938cb
18 changed files with 219 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/java/org/apache/commons/math/FunctionEvaluationException.java
View File

@ -17,7 +17,7 @@
package org.apache.commons.math;
/**
* Exeption thrown when an error occurs evaluating a function.
* Exception thrown when an error occurs evaluating a function.
* <p>
* Maintains an <code>argument</code> property holding the input value that
* caused the function evaluation to fail.

7
src/java/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.java
View File

@ -151,11 +151,14 @@ public abstract class AdaptiveStepsizeIntegrator
* function 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
*/
public void addSwitchingFunction(SwitchingFunction function,
double maxCheckInterval,
double convergence) {
switchesHandler.add(function, maxCheckInterval, convergence);
double convergence,
int maxIterationCount) {
switchesHandler.add(function, maxCheckInterval, convergence, maxIterationCount);
}
/** Perform some sanity checks on the integration parameters.

5
src/java/org/apache/commons/math/ode/FirstOrderIntegrator.java
View File

@ -58,10 +58,13 @@ public interface FirstOrderIntegrator {
* function 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
*/
public void addSwitchingFunction(SwitchingFunction function,
double maxCheckInterval,
double convergence);
double convergence,
int maxIterationCount);
/** Integrate the differential equations up to the given time.
* <p>This method solves an Initial Value Problem (IVP).</p>

7
src/java/org/apache/commons/math/ode/GraggBulirschStoerIntegrator.java
View File

@ -291,11 +291,14 @@ public class GraggBulirschStoerIntegrator
* function 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
*/
public void addSwitchingFunction(SwitchingFunction function,
double maxCheckInterval,
double convergence) {
super.addSwitchingFunction(function, maxCheckInterval, convergence);
double convergence,
int maxIterationCount) {
super.addSwitchingFunction(function, maxCheckInterval, convergence, maxIterationCount);
denseOutput = (handler.requiresDenseOutput()
|| (! switchesHandler.isEmpty()));

10
src/java/org/apache/commons/math/ode/IntegratorException.java
View File

@ -36,6 +36,14 @@ public class IntegratorException
super(specifier, parts);
}
private static final long serialVersionUID = -1390328069787882608L;
/**
* Create an exception with a given root cause.
* @param cause the exception or error that caused this exception to be thrown
*/
public IntegratorException(Throwable cause) {
super(cause);
}
private static final long serialVersionUID = -1215318282266670558L;
}

7
src/java/org/apache/commons/math/ode/RungeKuttaIntegrator.java
View File

@ -95,11 +95,14 @@ public abstract class RungeKuttaIntegrator
* function 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
*/
public void addSwitchingFunction(SwitchingFunction function,
double maxCheckInterval,
double convergence) {
switchesHandler.add(function, maxCheckInterval, convergence);
double convergence,
int maxIterationCount) {
switchesHandler.add(function, maxCheckInterval, convergence, maxIterationCount);
}
/** Perform some sanity checks on the integration parameters.

96
src/java/org/apache/commons/math/ode/SwitchState.java
View File

@ -42,7 +42,7 @@ import org.apache.commons.math.analysis.UnivariateRealSolver;
class SwitchState implements Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = 3256541562455482289L;
private static final long serialVersionUID = -7307007422156119622L;
/** Switching function. */
private SwitchingFunction function;
@ -53,6 +53,9 @@ class SwitchState implements Serializable {
/** Convergence threshold for event localisation. */
private double convergence;
/** Upper limit in the iteration count for event localisation. */
private int maxIterationCount;
/** Time at the beginning of the step. */
private double t0;
@ -85,12 +88,15 @@ class SwitchState implements Serializable {
* function 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
*/
public SwitchState(SwitchingFunction function,
double maxCheckInterval, double convergence) {
this.function = function;
this.maxCheckInterval = maxCheckInterval;
this.convergence = Math.abs(convergence);
public SwitchState(SwitchingFunction function, double maxCheckInterval,
double convergence, int maxIterationCount) {
this.function = function;
this.maxCheckInterval = maxCheckInterval;
this.convergence = Math.abs(convergence);
this.maxIterationCount = maxIterationCount;
// some dummy values ...
t0 = Double.NaN;
@ -109,8 +115,11 @@ class SwitchState implements Serializable {
* beginning of the step
* @param y0 array containing the current value of the state vector
* at the beginning of the step
* @exception FunctionEvaluationException if the switching function
* value cannot be evaluated at the beginning of the step
*/
public void reinitializeBegin(double t0, double[] y0) {
public void reinitializeBegin(double t0, double[] y0)
throws FunctionEvaluationException {
this.t0 = t0;
g0 = function.g(t0, y0);
g0Positive = (g0 >= 0);
@ -121,8 +130,14 @@ class SwitchState implements Serializable {
* @return true if the switching function triggers an event before
* the end of the proposed step (this implies the step should be
* rejected)
* @exception DerivativeException if the interpolator fails to
* compute the function somewhere within the step
* @exception FunctionEvaluationException if the switching function
* cannot be evaluated
* @exception ConvergenceException if an event cannot be located
*/
public boolean evaluateStep(final StepInterpolator interpolator) {
public boolean evaluateStep(final StepInterpolator interpolator)
throws DerivativeException, FunctionEvaluationException, ConvergenceException {
try {
@ -147,36 +162,32 @@ class SwitchState implements Serializable {
// variation direction, with respect to the integration direction
increasing = (gb >= ga);
try {
UnivariateRealSolver solver = new BrentSolver(new UnivariateRealFunction() {
public double value(double t) throws FunctionEvaluationException {
try {
interpolator.setInterpolatedTime(t);
return function.g(t, interpolator.getInterpolatedState());
} catch (DerivativeException e) {
throw new FunctionEvaluationException(t, e);
}
UnivariateRealSolver solver = new BrentSolver(new UnivariateRealFunction() {
public double value(double t) throws FunctionEvaluationException {
try {
interpolator.setInterpolatedTime(t);
return function.g(t, interpolator.getInterpolatedState());
} catch (DerivativeException e) {
throw new FunctionEvaluationException(t, e);
}
});
solver.setAbsoluteAccuracy(convergence);
solver.setMaximalIterationCount(1000);
double root = solver.solve(ta, tb);
if (Double.isNaN(previousEventTime) || (Math.abs(previousEventTime - root) > convergence)) {
pendingEventTime = root;
if (pendingEvent && (Math.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;
}
} catch (ConvergenceException ce) {
throw new RuntimeException("internal error");
});
solver.setAbsoluteAccuracy(convergence);
solver.setMaximalIterationCount(maxIterationCount);
double root = solver.solve(ta, tb);
if (Double.isNaN(previousEventTime) || (Math.abs(previousEventTime - root) > convergence)) {
pendingEventTime = root;
if (pendingEvent && (Math.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 {
@ -192,10 +203,12 @@ class SwitchState implements Serializable {
pendingEventTime = Double.NaN;
return false;
} catch (DerivativeException e) {
throw new RuntimeException("unexpected exception: " + e.getMessage());
} catch (FunctionEvaluationException e) {
throw new RuntimeException("unexpected exception: " + e.getMessage());
Throwable cause = e.getCause();
if ((cause != null) && (cause instanceof DerivativeException)) {
throw (DerivativeException) cause;
}
throw e;
}
}
@ -214,8 +227,11 @@ class SwitchState implements Serializable {
* end of the step
* @param y array containing the current value of the state vector
* at the end of the step
* @exception FunctionEvaluationException if the value of the switching
* function cannot be evaluated
*/
public void stepAccepted(double t, double[] y) {
public void stepAccepted(double t, double[] y)
throws FunctionEvaluationException {
t0 = t;
g0 = function.g(t, y);

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

@ -19,6 +19,8 @@ package org.apache.commons.math.ode;
import java.io.Serializable;
import org.apache.commons.math.FunctionEvaluationException;
/** This interface represents a switching function.
*
* <p>A switching function allows to handle discrete events in
@ -91,8 +93,10 @@ public interface SwitchingFunction extends Serializable {
* @param t current value of the independent <i>time</i> variable
* @param y array containing the current value of the state vector
* @return value of the g function
* @exception FunctionEvaluationException if the value of the function
* cannot be evaluated
*/
public double g(double t, double[] y);
public double g(double t, double[] y) throws FunctionEvaluationException;
/** Handle an event and choose what to do next.

36
src/java/org/apache/commons/math/ode/SwitchingFunctionsHandler.java
View File

@ -17,6 +17,8 @@
package org.apache.commons.math.ode;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.ode.DerivativeException;
import java.util.ArrayList;
@ -48,10 +50,13 @@ public class SwitchingFunctionsHandler {
* function 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
*/
public void add(SwitchingFunction function,
double maxCheckInterval, double convergence) {
functions.add(new SwitchState(function, maxCheckInterval, convergence));
public void add(SwitchingFunction function, double maxCheckInterval,
double convergence, int maxIterationCount) {
functions.add(new SwitchState(function, maxCheckInterval,
convergence, maxIterationCount));
}
/** Check if the handler does not have any condition.
@ -67,8 +72,12 @@ public class SwitchingFunctionsHandler {
* @return true if at least one switching function triggers an event
* before the end of the proposed step (this implies the step should
* be rejected)
* @exception DerivativeException if the interpolator fails to
* compute the function somewhere within the step
* @exception IntegratorException if an event cannot be located
*/
public boolean evaluateStep(StepInterpolator interpolator) {
public boolean evaluateStep(StepInterpolator interpolator)
throws DerivativeException, IntegratorException {
try {
@ -118,8 +127,10 @@ public class SwitchingFunctionsHandler {
return first != null;
} catch (DerivativeException e) {
throw new RuntimeException("unexpected exception: " + e.getMessage());
} catch (FunctionEvaluationException fee) {
throw new IntegratorException(fee);
} catch (ConvergenceException ce) {
throw new IntegratorException(ce);
}
}
@ -140,10 +151,17 @@ public class SwitchingFunctionsHandler {
* 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
* switching functions cannot be evaluated
*/
public void stepAccepted(double t, double[] y) {
for (Iterator iter = functions.iterator(); iter.hasNext();) {
((SwitchState) iter.next()).stepAccepted(t, y);
public void stepAccepted(double t, double[] y)
throws IntegratorException {
try {
for (Iterator iter = functions.iterator(); iter.hasNext();) {
((SwitchState) iter.next()).stepAccepted(t, y);
}
} catch (FunctionEvaluationException fee) {
throw new IntegratorException(fee);
}
}

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

@ -86,7 +86,7 @@ public class ClassicalRungeKuttaIntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);

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

@ -189,7 +189,7 @@ public class DormandPrince54IntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),

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

@ -142,7 +142,7 @@ public class DormandPrince853IntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
@ -227,7 +227,7 @@ public class DormandPrince853IntegratorTest
final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
FirstOrderIntegrator integ =
new DormandPrince853Integrator(0.1, 10, 1.0e-12, 0.0);
integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12, 1000);
double[] y = { Double.NaN };
integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
assertEquals(8.0, y[0], 1.0e-12);

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

@ -64,7 +64,7 @@ public class EulerIntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),

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

@ -66,7 +66,7 @@ public class GillIntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@ -138,7 +138,7 @@ public class GillIntegratorTest
throws DerivativeException, IntegratorException {
final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
FirstOrderIntegrator integ = new GillIntegrator(0.3);
integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12, 1000);
double[] y = { Double.NaN };
integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
assertEquals(8.0, y[0], 1.0e-12);

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

@ -182,7 +182,7 @@ public class GraggBulirschStoerIntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
@ -238,7 +238,7 @@ public class GraggBulirschStoerIntegratorTest
final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
FirstOrderIntegrator integ =
new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12, 1000);
double[] y = { Double.NaN };
integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
assertEquals(8.0, y[0], 1.0e-12);

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

@ -17,6 +17,8 @@
package org.apache.commons.math.ode;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.HighamHall54Integrator;
@ -149,7 +151,7 @@ public class HighamHall54IntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep);
Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
@ -161,6 +163,94 @@ public class HighamHall54IntegratorTest
}
public void testSwitchingFunctionsError()
throws DerivativeException, IntegratorException {
final TestProblem1 pb = new TestProblem1();
double minStep = 0;
double maxStep = pb.getFinalTime() - pb.getInitialTime();
double scalAbsoluteTolerance = 1.0e-8;
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
FirstOrderIntegrator integ =
new HighamHall54Integrator(minStep, maxStep,
scalAbsoluteTolerance, scalRelativeTolerance);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.setStepHandler(handler);
integ.addSwitchingFunction(new SwitchingFunction() {
public int eventOccurred(double t, double[] y) {
return SwitchingFunction.CONTINUE;
}
public double g(double t, double[] y) throws FunctionEvaluationException {
double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
double offset = t - middle;
if (offset > 0) {
throw new FunctionEvaluationException(t);
}
return offset;
}
public void resetState(double t, double[] y) {
}
private static final long serialVersionUID = 935652725339916361L;
}, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
try {
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
fail("an exception should have been thrown");
} catch (IntegratorException ie) {
// expected behavior
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testSwitchingFunctionsNoConvergence()
throws DerivativeException, IntegratorException {
final TestProblem1 pb = new TestProblem1();
double minStep = 0;
double maxStep = pb.getFinalTime() - pb.getInitialTime();
double scalAbsoluteTolerance = 1.0e-8;
double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
FirstOrderIntegrator integ =
new HighamHall54Integrator(minStep, maxStep,
scalAbsoluteTolerance, scalRelativeTolerance);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.setStepHandler(handler);
integ.addSwitchingFunction(new SwitchingFunction() {
public int eventOccurred(double t, double[] y) {
return SwitchingFunction.CONTINUE;
}
public double g(double t, double[] y) {
double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
double offset = t - middle;
return (offset > 0) ? (offset + 0.5) : (offset - 0.5);
}
public void resetState(double t, double[] y) {
}
private static final long serialVersionUID = 935652725339916361L;
}, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
try {
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
fail("an exception should have been thrown");
} catch (IntegratorException ie) {
assertTrue(ie.getCause() != null);
assertTrue(ie.getCause() instanceof ConvergenceException);
} catch (Exception e) {
fail("wrong exception type caught");
}
}
public void testSanityChecks() {
try {
final TestProblem3 pb = new TestProblem3(0.9);

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

@ -63,7 +63,7 @@ public class MidpointIntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
}
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),

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

@ -66,7 +66,7 @@ public class ThreeEighthesIntegratorTest
SwitchingFunction[] functions = pb.getSwitchingFunctions();
for (int l = 0; l < functions.length; ++l) {
integ.addSwitchingFunction(functions[l],
Double.POSITIVE_INFINITY, 1.0e-6 * step);
Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
}
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);