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