in commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ode/events/EventState.java [221:334]
public boolean evaluateStep(final StepInterpolator interpolator)
throws MaxCountExceededException, NoBracketingException {
try {
forward = interpolator.isForward();
final double t1 = interpolator.getCurrentTime();
final double dt = t1 - t0;
if (JdkMath.abs(dt) < convergence) {
// we cannot do anything on such a small step, don't trigger any events
return false;
}
final int n = JdkMath.max(1, (int) JdkMath.ceil(JdkMath.abs(dt) / maxCheckInterval));
final double h = dt / n;
final UnivariateFunction f = new UnivariateFunction() {
/** {@inheritDoc} */
@Override
public double value(final double t) throws LocalMaxCountExceededException {
try {
interpolator.setInterpolatedTime(t);
return handler.g(t, getCompleteState(interpolator));
} catch (MaxCountExceededException mcee) {
throw new LocalMaxCountExceededException(mcee);
}
}
};
double ta = t0;
double ga = g0;
for (int i = 0; i < n; ++i) {
// evaluate handler value at the end of the substep
final double tb = (i == n - 1) ? t1 : t0 + (i + 1) * h;
interpolator.setInterpolatedTime(tb);
final double gb = handler.g(tb, getCompleteState(interpolator));
// check events occurrence
if (g0Positive ^ (gb >= 0)) {
// there is a sign change: an event is expected during this step
// variation direction, with respect to the integration direction
increasing = gb >= ga;
// find the event time making sure we select a solution just at or past the exact root
final double root;
if (solver instanceof BracketedUnivariateSolver<?>) {
@SuppressWarnings("unchecked")
BracketedUnivariateSolver<UnivariateFunction> bracketing =
(BracketedUnivariateSolver<UnivariateFunction>) solver;
root = forward ?
bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
} else {
final double baseRoot = forward ?
solver.solve(maxIterationCount, f, ta, tb) :
solver.solve(maxIterationCount, f, tb, ta);
final int remainingEval = maxIterationCount - solver.getEvaluations();
BracketedUnivariateSolver<UnivariateFunction> bracketing =
new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
root = forward ?
UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
}
if (!Double.isNaN(previousEventTime) &&
JdkMath.abs(root - ta) <= convergence &&
JdkMath.abs(root - previousEventTime) <= convergence) {
// we have either found nothing or found (again ?) a past event,
// retry the substep excluding this value, and taking care to have the
// required sign in case the g function is noisy around its zero and
// crosses the axis several times
do {
ta = forward ? ta + convergence : ta - convergence;
ga = f.value(ta);
} while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb)));
if (forward ^ (ta >= tb)) {
// we were able to skip this spurious root
--i;
} else {
// we can't avoid this root before the end of the step,
// we have to handle it despite it is close to the former one
// maybe we have two very close roots
pendingEventTime = root;
pendingEvent = true;
return true;
}
} else if (Double.isNaN(previousEventTime) ||
JdkMath.abs(previousEventTime - root) > convergence) {
pendingEventTime = root;
pendingEvent = true;
return true;
} else {
// no sign change: there is no event for now
ta = tb;
ga = gb;
}
} else {
// no sign change: there is no event for now
ta = tb;
ga = gb;
}
}
// no event during the whole step
pendingEvent = false;
pendingEventTime = Double.NaN;
return false;
} catch (LocalMaxCountExceededException lmcee) {
throw lmcee.getException();
}
}