1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.ode.events; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ConvergenceException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.FunctionEvaluationException; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.UnivariateRealFunction; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.solvers.BrentSolver; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.MathInternalError; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.DerivativeException; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.StepInterpolator; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** This class handles the state for one {@link EventHandler 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * event handler} during integration steps. 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>Each time the integrator proposes a step, the event handler 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * switching function should be checked. This class handles the state 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * of one handler during one integration step, with references to the 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * state at the end of the preceding step. This information is used to 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * decide if the handler should trigger an event or not during the 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * proposed step (and hence the step should be reduced to ensure the 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * event occurs at a bound rather than inside the step).</p> 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class EventState { 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Event handler. */ 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final EventHandler handler; 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Maximal time interval between events handler checks. */ 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double maxCheckInterval; 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Convergence threshold for event localization. */ 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double convergence; 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Upper limit in the iteration count for event localization. */ 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final int maxIterationCount; 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Time at the beginning of the step. */ 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double t0; 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Value of the events handler at the beginning of the step. */ 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double g0; 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simulated sign of g0 (we cheat when crossing events). */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean g0Positive; 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Indicator of event expected during the step. */ 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean pendingEvent; 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Occurrence time of the pending event. */ 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double pendingEventTime; 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Occurrence time of the previous event. */ 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double previousEventTime; 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Integration direction. */ 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean forward; 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Variation direction around pending event. 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (this is considered with respect to the integration direction) 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean increasing; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Next action indicator. */ 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int nextAction; 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param handler event handler 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param maxCheckInterval maximal time interval between switching 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * function checks (this interval prevents missing sign changes in 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * case the integration steps becomes very large) 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param convergence convergence threshold in the event time search 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param maxIterationCount upper limit of the iteration count in 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the event time search 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public EventState(final EventHandler handler, final double maxCheckInterval, 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double convergence, final int maxIterationCount) { 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.handler = handler; 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.maxCheckInterval = maxCheckInterval; 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.convergence = FastMath.abs(convergence); 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.maxIterationCount = maxIterationCount; 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // some dummy values ... 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond t0 = Double.NaN; 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0 = Double.NaN; 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0Positive = true; 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEvent = false; 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEventTime = Double.NaN; 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond previousEventTime = Double.NaN; 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond increasing = true; 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond nextAction = EventHandler.CONTINUE; 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the underlying event handler. 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return underlying event handler 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public EventHandler getEventHandler() { 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return handler; 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the maximal time interval between events handler checks. 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return maximal time interval between events handler checks 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getMaxCheckInterval() { 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return maxCheckInterval; 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the convergence threshold for event localization. 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return convergence threshold for event localization 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getConvergence() { 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return convergence; 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the upper limit in the iteration count for event localization. 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return upper limit in the iteration count for event localization 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public int getMaxIterationCount() { 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return maxIterationCount; 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reinitialize the beginning of the step. 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param interpolator valid for the current step 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EventException if the event handler 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * value cannot be evaluated at the beginning of the step 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void reinitializeBegin(final StepInterpolator interpolator) 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EventException { 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // excerpt from MATH-421 issue: 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // If an ODE solver is setup with an EventHandler that return STOP 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // when the even is triggered, the integrator stops (which is exactly 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the expected behavior). If however the user want to restart the 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // solver from the final state reached at the event with the same 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // configuration (expecting the event to be triggered again at a 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // later time), then the integrator may fail to start. It can get stuck 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // at the previous event. 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // The use case for the bug MATH-421 is fairly general, so events occurring 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // less than epsilon after the solver start in the first step should be ignored, 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // where epsilon is the convergence threshold of the event. The sign of the g 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // function should be evaluated after this initial ignore zone, not exactly at 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // beginning (if there are no event at the very beginning g(t0) and g(t0+epsilon) 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // have the same sign, so this does not hurt ; if there is an event at the very 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // beginning, g(t0) and g(t0+epsilon) have opposite signs and we want to start 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // with the second one. Of course, the sign of epsilon depend on the integration 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // direction (forward or backward). This explains what is done below. 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double ignoreZone = interpolator.isForward() ? getConvergence() : -getConvergence(); 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond t0 = interpolator.getPreviousTime() + ignoreZone; 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setInterpolatedTime(t0); 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0 = handler.g(t0, interpolator.getInterpolatedState()); 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (g0 == 0) { 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // extremely rare case: there is a zero EXACTLY at end of ignore zone 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we will use the opposite of sign at step beginning to force ignoring this zero 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double tStart = interpolator.getPreviousTime(); 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setInterpolatedTime(tStart); 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0Positive = handler.g(tStart, interpolator.getInterpolatedState()) <= 0; 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0Positive = g0 >= 0; 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (DerivativeException mue) { 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new EventException(mue); 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Evaluate the impact of the proposed step on the event handler. 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param interpolator step interpolator for the proposed step 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return true if the event handler triggers an event before 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the end of the proposed step 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception DerivativeException if the interpolator fails to 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * compute the switching function somewhere within the step 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EventException if the switching function 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * cannot be evaluated 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception ConvergenceException if an event cannot be located 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public boolean evaluateStep(final StepInterpolator interpolator) 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DerivativeException, EventException, ConvergenceException { 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond forward = interpolator.isForward(); 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t1 = interpolator.getCurrentTime(); 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.abs(t1 - t0) < convergence) { 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we cannot do anything on such a small step, don't trigger any events 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return false; 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double start = forward ? (t0 + convergence) : t0 - convergence; 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dt = t1 - start; 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval)); 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double h = dt / n; 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double ta = t0; 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double ga = g0; 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // evaluate handler value at the end of the substep 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double tb = start + (i + 1) * h; 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setInterpolatedTime(tb); 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double gb = handler.g(tb, interpolator.getInterpolatedState()); 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // check events occurrence 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (g0Positive ^ (gb >= 0)) { 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // there is a sign change: an event is expected during this step 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // variation direction, with respect to the integration direction 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond increasing = gb >= ga; 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final UnivariateRealFunction f = new UnivariateRealFunction() { 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(final double t) { 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolator.setInterpolatedTime(t); 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return handler.g(t, interpolator.getInterpolatedState()); 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (DerivativeException e) { 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new EmbeddedDerivativeException(e); 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (EventException e) { 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new EmbeddedEventException(e); 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final BrentSolver solver = new BrentSolver(convergence); 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (ga * gb >= 0) { 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // this is a corner case: 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // - there was an event near ta, 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // - there is another event between ta and tb 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // - when ta was computed, convergence was reached on the "wrong side" of the interval 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // this implies that the real sign of ga is the same as gb, so we need to slightly 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // shift ta to make sure ga and gb get opposite signs and the solver won't complain 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // about bracketing 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double epsilon = (forward ? 0.25 : -0.25) * convergence; 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; (k < 4) && (ga * gb > 0); ++k) { 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ta += epsilon; 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ga = f.value(ta); 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (FunctionEvaluationException ex) { 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DerivativeException(ex); 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (ga * gb > 0) { 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // this should never happen 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathInternalError(); 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double root; 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond root = (ta <= tb) ? 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond solver.solve(maxIterationCount, f, ta, tb) : 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond solver.solve(maxIterationCount, f, tb, ta); 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (FunctionEvaluationException ex) { 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DerivativeException(ex); 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if ((!Double.isNaN(previousEventTime)) && 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond (FastMath.abs(root - ta) <= convergence) && 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond (FastMath.abs(root - previousEventTime) <= convergence)) { 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we have either found nothing or found (again ?) a past event, we simply ignore it 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ta = tb; 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ga = gb; 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else if (Double.isNaN(previousEventTime) || 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond (FastMath.abs(previousEventTime - root) > convergence)) { 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEventTime = root; 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEvent = true; 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return true; 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // no sign change: there is no event for now 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ta = tb; 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ga = gb; 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // no sign change: there is no event for now 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ta = tb; 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ga = gb; 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // no event during the whole step 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEvent = false; 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEventTime = Double.NaN; 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return false; 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (EmbeddedDerivativeException ede) { 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw ede.getDerivativeException(); 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (EmbeddedEventException eee) { 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw eee.getEventException(); 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the occurrence time of the event triggered in the current step. 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return occurrence time of the event triggered in the current 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * step or positive infinity if no events are triggered 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getEventTime() { 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return pendingEvent ? pendingEventTime : Double.POSITIVE_INFINITY; 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Acknowledge the fact the step has been accepted by the integrator. 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param t value of the independent <i>time</i> variable at the 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * end of the step 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y array containing the current value of the state vector 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * at the end of the step 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EventException if the value of the event 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * handler cannot be evaluated 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void stepAccepted(final double t, final double[] y) 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EventException { 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond t0 = t; 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0 = handler.g(t, y); 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) { 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // force the sign to its value "just after the event" 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond previousEventTime = t; 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0Positive = increasing; 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond nextAction = handler.eventOccurred(t, y, !(increasing ^ forward)); 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond g0Positive = g0 >= 0; 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond nextAction = EventHandler.CONTINUE; 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Check if the integration should be stopped at the end of the 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * current step. 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return true if the integration should be stopped 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public boolean stop() { 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return nextAction == EventHandler.STOP; 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Let the event handler reset the state if it wants. 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param t value of the independent <i>time</i> variable at the 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * beginning of the next step 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y array were to put the desired state vector at the beginning 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * of the next step 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return true if the integrator should reset the derivatives too 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception EventException if the state cannot be reseted by the event 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * handler 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public boolean reset(final double t, final double[] y) 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws EventException { 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) { 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return false; 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (nextAction == EventHandler.RESET_STATE) { 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond handler.resetState(t, y); 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEvent = false; 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond pendingEventTime = Double.NaN; 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return (nextAction == EventHandler.RESET_STATE) || 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond (nextAction == EventHandler.RESET_DERIVATIVES); 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 380dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 381dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 382dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Local exception for embedding DerivativeException. */ 383dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static class EmbeddedDerivativeException extends RuntimeException { 384dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 385dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Serializable UID. */ 386dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final long serialVersionUID = 3574188382434584610L; 387dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 388dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Embedded exception. */ 389dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final DerivativeException derivativeException; 390dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 391dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 392dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param derivativeException embedded exception 393dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 394dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public EmbeddedDerivativeException(final DerivativeException derivativeException) { 395dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.derivativeException = derivativeException; 396dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 397dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 398dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the embedded exception. 399dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return embedded exception 400dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 401dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public DerivativeException getDerivativeException() { 402dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return derivativeException; 403dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 404dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 405dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 406dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 407dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Local exception for embedding EventException. */ 408dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static class EmbeddedEventException extends RuntimeException { 409dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 410dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Serializable UID. */ 411dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final long serialVersionUID = -1337749250090455474L; 412dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 413dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Embedded exception. */ 414dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final EventException eventException; 415dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 416dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 417dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param eventException embedded exception 418dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 419dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public EmbeddedEventException(final EventException eventException) { 420dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.eventException = eventException; 421dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 422dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 423dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the embedded exception. 424dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return embedded exception 425dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 426dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public EventException getEventException() { 427dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return eventException; 428dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 429dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 430dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 431dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 432