All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.math3.ode.events.EventState Maven / Gradle / Ivy

There is a newer version: 2.12.15
Show newest version
/*
 * 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.math3.ode.events;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.solvers.AllowedSolution;
import org.apache.commons.math3.analysis.solvers.BracketedUnivariateSolver;
import org.apache.commons.math3.analysis.solvers.PegasusSolver;
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;

/** This class handles the state for one {@link EventHandler
 * event handler} during integration steps.
 *
 * 

Each time the integrator proposes a step, the event handler * switching function should be checked. This class handles the state * 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.

* * @since 1.2 */ public class EventState { /** Event handler. */ private final EventHandler handler; /** Maximal time interval between events handler checks. */ private final double maxCheckInterval; /** Convergence threshold for event localization. */ private final double convergence; /** 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; /** Value of the events handler at the beginning of the step. */ private double g0; /** Simulated sign of g0 (we cheat when crossing events). */ private boolean g0Positive; /** Indicator of event expected during the step. */ private boolean pendingEvent; /** Occurrence time of the pending event. */ private double pendingEventTime; /** Occurrence time of the previous event. */ private double previousEventTime; /** Integration direction. */ private boolean forward; /** Variation direction around pending event. * (this is considered with respect to the integration direction) */ private boolean increasing; /** Next action indicator. */ private EventHandler.Action nextAction; /** Root-finding algorithm to use to detect state events. */ private final UnivariateSolver solver; /** Simple constructor. * @param handler event handler * @param maxCheckInterval maximal time interval between switching * 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 * @param solver Root-finding algorithm to use to detect state events */ public EventState(final EventHandler handler, final double maxCheckInterval, final double convergence, final int maxIterationCount, final UnivariateSolver solver) { this.handler = handler; this.maxCheckInterval = maxCheckInterval; this.convergence = FastMath.abs(convergence); this.maxIterationCount = maxIterationCount; this.solver = solver; // some dummy values ... expandable = null; t0 = Double.NaN; g0 = Double.NaN; g0Positive = true; pendingEvent = false; pendingEventTime = Double.NaN; previousEventTime = Double.NaN; increasing = true; nextAction = EventHandler.Action.CONTINUE; } /** Get the underlying event handler. * @return underlying event handler */ public EventHandler getEventHandler() { 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 */ public double getMaxCheckInterval() { return maxCheckInterval; } /** Get the convergence threshold for event localization. * @return convergence threshold for event localization */ public double getConvergence() { return convergence; } /** Get the upper limit in the iteration count for event localization. * @return upper limit in the iteration count for event localization */ public int getMaxIterationCount() { return maxIterationCount; } /** Reinitialize the beginning of the step. * @param interpolator valid for the current step * @exception MaxCountExceededException if the interpolator throws one because * the number of functions evaluations is exceeded */ public void reinitializeBegin(final StepInterpolator interpolator) throws MaxCountExceededException { t0 = interpolator.getPreviousTime(); interpolator.setInterpolatedTime(t0); 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 // when the even is triggered, the integrator stops (which is exactly // the expected behavior). If however the user wants to restart the // solver from the final state reached at the event with the same // configuration (expecting the event to be triggered again at a // later time), then the integrator may fail to start. It can get stuck // at the previous event. The use case for the bug MATH-421 is fairly // general, so events occurring exactly at start in the first step should // be ignored. // extremely rare case: there is a zero EXACTLY at interval start // we will use the sign slightly after step beginning to force ignoring this zero final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(), FastMath.abs(solver.getRelativeAccuracy() * t0)); final double tStart = t0 + 0.5 * epsilon; interpolator.setInterpolatedTime(tStart); 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 * the end of the proposed step * @exception MaxCountExceededException if the interpolator throws one because * the number of functions evaluations is exceeded * @exception NoBracketingException if the event cannot be bracketed */ public boolean evaluateStep(final StepInterpolator interpolator) throws MaxCountExceededException, NoBracketingException { try { forward = interpolator.isForward(); final double t1 = interpolator.getCurrentTime(); final double dt = t1 - t0; if (FastMath.abs(dt) < convergence) { // we cannot do anything on such a small step, don't trigger any events return false; } final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval)); final double h = dt / n; final UnivariateFunction f = new UnivariateFunction() { /** {@inheritDoc} */ 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 bracketing = (BracketedUnivariateSolver) 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 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)) && (FastMath.abs(root - ta) <= convergence) && (FastMath.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) || (FastMath.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(); } } /** Get the occurrence time of the event triggered in the current step. * @return occurrence time of the event triggered in the current * step or infinity if no events are triggered */ public double getEventTime() { return pendingEvent ? pendingEventTime : (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY); } /** Acknowledge the fact the step has been accepted by the integrator. * @param t value of the independent time variable at the * end of the step * @param y array containing the current value of the state vector * at the end of the step */ public void stepAccepted(final double t, final double[] y) { t0 = t; g0 = handler.g(t, y); if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) { // force the sign to its value "just after the event" previousEventTime = t; g0Positive = increasing; nextAction = handler.eventOccurred(t, y, !(increasing ^ forward)); } else { g0Positive = g0 >= 0; nextAction = EventHandler.Action.CONTINUE; } } /** 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() { return nextAction == EventHandler.Action.STOP; } /** Let the event handler reset the state if it wants. * @param t value of the independent time 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 */ public boolean reset(final double t, final double[] y) { if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) { return false; } if (nextAction == EventHandler.Action.RESET_STATE) { handler.resetState(t, y); } pendingEvent = false; pendingEventTime = Double.NaN; return (nextAction == EventHandler.Action.RESET_STATE) || (nextAction == EventHandler.Action.RESET_DERIVATIVES); } /** Local wrapper to propagate exceptions. */ private static class LocalMaxCountExceededException extends RuntimeException { /** Serializable UID. */ private static final long serialVersionUID = 20120901L; /** Wrapped exception. */ private final MaxCountExceededException wrapped; /** Simple constructor. * @param exception exception to wrap */ LocalMaxCountExceededException(final MaxCountExceededException exception) { wrapped = exception; } /** Get the wrapped exception. * @return wrapped exception */ public MaxCountExceededException getException() { return wrapped; } } }




© 2015 - 2024 Weber Informatics LLC | Privacy Policy