001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.ode.events; 019 020 import org.apache.commons.math.ConvergenceException; 021 import org.apache.commons.math.FunctionEvaluationException; 022 import org.apache.commons.math.MathRuntimeException; 023 import org.apache.commons.math.analysis.UnivariateRealFunction; 024 import org.apache.commons.math.analysis.solvers.BrentSolver; 025 import org.apache.commons.math.ode.DerivativeException; 026 import org.apache.commons.math.ode.sampling.StepInterpolator; 027 028 /** This class handles the state for one {@link EventHandler 029 * event handler} during integration steps. 030 * 031 * <p>Each time the integrator proposes a step, the event handler 032 * switching function should be checked. This class handles the state 033 * of one handler during one integration step, with references to the 034 * state at the end of the preceding step. This information is used to 035 * decide if the handler should trigger an event or not during the 036 * proposed step (and hence the step should be reduced to ensure the 037 * event occurs at a bound rather than inside the step).</p> 038 * 039 * @version $Revision: 889006 $ $Date: 2009-12-09 17:46:36 -0500 (Wed, 09 Dec 2009) $ 040 * @since 1.2 041 */ 042 public class EventState { 043 044 /** Event handler. */ 045 private final EventHandler handler; 046 047 /** Maximal time interval between events handler checks. */ 048 private final double maxCheckInterval; 049 050 /** Convergence threshold for event localization. */ 051 private final double convergence; 052 053 /** Upper limit in the iteration count for event localization. */ 054 private final int maxIterationCount; 055 056 /** Time at the beginning of the step. */ 057 private double t0; 058 059 /** Value of the events handler at the beginning of the step. */ 060 private double g0; 061 062 /** Simulated sign of g0 (we cheat when crossing events). */ 063 private boolean g0Positive; 064 065 /** Indicator of event expected during the step. */ 066 private boolean pendingEvent; 067 068 /** Occurrence time of the pending event. */ 069 private double pendingEventTime; 070 071 /** Occurrence time of the previous event. */ 072 private double previousEventTime; 073 074 /** Integration direction. */ 075 private boolean forward; 076 077 /** Variation direction around pending event. 078 * (this is considered with respect to the integration direction) 079 */ 080 private boolean increasing; 081 082 /** Next action indicator. */ 083 private int nextAction; 084 085 /** Simple constructor. 086 * @param handler event handler 087 * @param maxCheckInterval maximal time interval between switching 088 * function checks (this interval prevents missing sign changes in 089 * case the integration steps becomes very large) 090 * @param convergence convergence threshold in the event time search 091 * @param maxIterationCount upper limit of the iteration count in 092 * the event time search 093 */ 094 public EventState(final EventHandler handler, final double maxCheckInterval, 095 final double convergence, final int maxIterationCount) { 096 this.handler = handler; 097 this.maxCheckInterval = maxCheckInterval; 098 this.convergence = Math.abs(convergence); 099 this.maxIterationCount = maxIterationCount; 100 101 // some dummy values ... 102 t0 = Double.NaN; 103 g0 = Double.NaN; 104 g0Positive = true; 105 pendingEvent = false; 106 pendingEventTime = Double.NaN; 107 previousEventTime = Double.NaN; 108 increasing = true; 109 nextAction = EventHandler.CONTINUE; 110 111 } 112 113 /** Get the underlying event handler. 114 * @return underlying event handler 115 */ 116 public EventHandler getEventHandler() { 117 return handler; 118 } 119 120 /** Get the maximal time interval between events handler checks. 121 * @return maximal time interval between events handler checks 122 */ 123 public double getMaxCheckInterval() { 124 return maxCheckInterval; 125 } 126 127 /** Get the convergence threshold for event localization. 128 * @return convergence threshold for event localization 129 */ 130 public double getConvergence() { 131 return convergence; 132 } 133 134 /** Get the upper limit in the iteration count for event localization. 135 * @return upper limit in the iteration count for event localization 136 */ 137 public int getMaxIterationCount() { 138 return maxIterationCount; 139 } 140 141 /** Reinitialize the beginning of the step. 142 * @param tStart value of the independent <i>time</i> variable at the 143 * beginning of the step 144 * @param yStart array containing the current value of the state vector 145 * at the beginning of the step 146 * @exception EventException if the event handler 147 * value cannot be evaluated at the beginning of the step 148 */ 149 public void reinitializeBegin(final double tStart, final double[] yStart) 150 throws EventException { 151 t0 = tStart; 152 g0 = handler.g(tStart, yStart); 153 g0Positive = g0 >= 0; 154 } 155 156 /** Evaluate the impact of the proposed step on the event handler. 157 * @param interpolator step interpolator for the proposed step 158 * @return true if the event handler triggers an event before 159 * the end of the proposed step (this implies the step should be 160 * rejected) 161 * @exception DerivativeException if the interpolator fails to 162 * compute the switching function somewhere within the step 163 * @exception EventException if the switching function 164 * cannot be evaluated 165 * @exception ConvergenceException if an event cannot be located 166 */ 167 public boolean evaluateStep(final StepInterpolator interpolator) 168 throws DerivativeException, EventException, ConvergenceException { 169 170 try { 171 172 forward = interpolator.isForward(); 173 final double t1 = interpolator.getCurrentTime(); 174 final int n = Math.max(1, (int) Math.ceil(Math.abs(t1 - t0) / maxCheckInterval)); 175 final double h = (t1 - t0) / n; 176 177 double ta = t0; 178 double ga = g0; 179 double tb = t0 + (interpolator.isForward() ? convergence : -convergence); 180 for (int i = 0; i < n; ++i) { 181 182 // evaluate handler value at the end of the substep 183 tb += h; 184 interpolator.setInterpolatedTime(tb); 185 final double gb = handler.g(tb, interpolator.getInterpolatedState()); 186 187 // check events occurrence 188 if (g0Positive ^ (gb >= 0)) { 189 // there is a sign change: an event is expected during this step 190 191 if (ga * gb > 0) { 192 // this is a corner case: 193 // - there was an event near ta, 194 // - there is another event between ta and tb 195 // - when ta was computed, convergence was reached on the "wrong side" of the interval 196 // this implies that the real sign of ga is the same as gb, so we need to slightly 197 // shift ta to make sure ga and gb get opposite signs and the solver won't complain 198 // about bracketing 199 final double epsilon = (forward ? 0.25 : -0.25) * convergence; 200 for (int k = 0; (k < 4) && (ga * gb > 0); ++k) { 201 ta += epsilon; 202 interpolator.setInterpolatedTime(ta); 203 ga = handler.g(ta, interpolator.getInterpolatedState()); 204 } 205 if (ga * gb > 0) { 206 // this should never happen 207 throw MathRuntimeException.createInternalError(null); 208 } 209 } 210 211 // variation direction, with respect to the integration direction 212 increasing = gb >= ga; 213 214 final UnivariateRealFunction f = new UnivariateRealFunction() { 215 public double value(final double t) throws FunctionEvaluationException { 216 try { 217 interpolator.setInterpolatedTime(t); 218 return handler.g(t, interpolator.getInterpolatedState()); 219 } catch (DerivativeException e) { 220 throw new FunctionEvaluationException(e, t); 221 } catch (EventException e) { 222 throw new FunctionEvaluationException(e, t); 223 } 224 } 225 }; 226 final BrentSolver solver = new BrentSolver(); 227 solver.setAbsoluteAccuracy(convergence); 228 solver.setMaximalIterationCount(maxIterationCount); 229 final double root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta); 230 if ((Math.abs(root - ta) <= convergence) && 231 (Math.abs(root - previousEventTime) <= convergence)) { 232 // we have either found nothing or found (again ?) a past event, we simply ignore it 233 ta = tb; 234 ga = gb; 235 } else if (Double.isNaN(previousEventTime) || 236 (Math.abs(previousEventTime - root) > convergence)) { 237 pendingEventTime = root; 238 if (pendingEvent && (Math.abs(t1 - pendingEventTime) <= convergence)) { 239 // we were already waiting for this event which was 240 // found during a previous call for a step that was 241 // rejected, this step must now be accepted since it 242 // properly ends exactly at the event occurrence 243 return false; 244 } 245 // either we were not waiting for the event or it has 246 // moved in such a way the step cannot be accepted 247 pendingEvent = true; 248 return true; 249 } 250 251 } else { 252 // no sign change: there is no event for now 253 ta = tb; 254 ga = gb; 255 } 256 257 } 258 259 // no event during the whole step 260 pendingEvent = false; 261 pendingEventTime = Double.NaN; 262 return false; 263 264 } catch (FunctionEvaluationException e) { 265 final Throwable cause = e.getCause(); 266 if ((cause != null) && (cause instanceof DerivativeException)) { 267 throw (DerivativeException) cause; 268 } else if ((cause != null) && (cause instanceof EventException)) { 269 throw (EventException) cause; 270 } 271 throw new EventException(e); 272 } 273 274 } 275 276 /** Get the occurrence time of the event triggered in the current 277 * step. 278 * @return occurrence time of the event triggered in the current 279 * step. 280 */ 281 public double getEventTime() { 282 return pendingEventTime; 283 } 284 285 /** Acknowledge the fact the step has been accepted by the integrator. 286 * @param t value of the independent <i>time</i> variable at the 287 * end of the step 288 * @param y array containing the current value of the state vector 289 * at the end of the step 290 * @exception EventException if the value of the event 291 * handler cannot be evaluated 292 */ 293 public void stepAccepted(final double t, final double[] y) 294 throws EventException { 295 296 t0 = t; 297 g0 = handler.g(t, y); 298 299 if (pendingEvent) { 300 // force the sign to its value "just after the event" 301 previousEventTime = t; 302 g0Positive = increasing; 303 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward)); 304 } else { 305 g0Positive = g0 >= 0; 306 nextAction = EventHandler.CONTINUE; 307 } 308 } 309 310 /** Check if the integration should be stopped at the end of the 311 * current step. 312 * @return true if the integration should be stopped 313 */ 314 public boolean stop() { 315 return nextAction == EventHandler.STOP; 316 } 317 318 /** Let the event handler reset the state if it wants. 319 * @param t value of the independent <i>time</i> variable at the 320 * beginning of the next step 321 * @param y array were to put the desired state vector at the beginning 322 * of the next step 323 * @return true if the integrator should reset the derivatives too 324 * @exception EventException if the state cannot be reseted by the event 325 * handler 326 */ 327 public boolean reset(final double t, final double[] y) 328 throws EventException { 329 330 if (! pendingEvent) { 331 return false; 332 } 333 334 if (nextAction == EventHandler.RESET_STATE) { 335 handler.resetState(t, y); 336 } 337 pendingEvent = false; 338 pendingEventTime = Double.NaN; 339 340 return (nextAction == EventHandler.RESET_STATE) || 341 (nextAction == EventHandler.RESET_DERIVATIVES); 342 343 } 344 345 }