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.analysis.UnivariateRealFunction; 023 import org.apache.commons.math.analysis.solvers.BrentSolver; 024 import org.apache.commons.math.ode.DerivativeException; 025 import org.apache.commons.math.ode.sampling.StepInterpolator; 026 027 /** This class handles the state for one {@link EventHandler 028 * event handler} during integration steps. 029 * 030 * <p>Each time the integrator proposes a step, the event handler 031 * switching function should be checked. This class handles the state 032 * of one handler during one integration step, with references to the 033 * state at the end of the preceding step. This information is used to 034 * decide if the handler should trigger an event or not during the 035 * proposed step (and hence the step should be reduced to ensure the 036 * event occurs at a bound rather than inside the step).</p> 037 * 038 * @version $Revision: 786881 $ $Date: 2009-06-20 14:53:08 -0400 (Sat, 20 Jun 2009) $ 039 * @since 1.2 040 */ 041 public class EventState { 042 043 /** Event handler. */ 044 private final EventHandler handler; 045 046 /** Maximal time interval between events handler checks. */ 047 private final double maxCheckInterval; 048 049 /** Convergence threshold for event localization. */ 050 private final double convergence; 051 052 /** Upper limit in the iteration count for event localization. */ 053 private final int maxIterationCount; 054 055 /** Time at the beginning of the step. */ 056 private double t0; 057 058 /** Value of the events handler at the beginning of the step. */ 059 private double g0; 060 061 /** Simulated sign of g0 (we cheat when crossing events). */ 062 private boolean g0Positive; 063 064 /** Indicator of event expected during the step. */ 065 private boolean pendingEvent; 066 067 /** Occurrence time of the pending event. */ 068 private double pendingEventTime; 069 070 /** Occurrence time of the previous event. */ 071 private double previousEventTime; 072 073 /** Integration direction. */ 074 private boolean forward; 075 076 /** Variation direction around pending event. 077 * (this is considered with respect to the integration direction) 078 */ 079 private boolean increasing; 080 081 /** Next action indicator. */ 082 private int nextAction; 083 084 /** Simple constructor. 085 * @param handler event handler 086 * @param maxCheckInterval maximal time interval between switching 087 * function checks (this interval prevents missing sign changes in 088 * case the integration steps becomes very large) 089 * @param convergence convergence threshold in the event time search 090 * @param maxIterationCount upper limit of the iteration count in 091 * the event time search 092 */ 093 public EventState(final EventHandler handler, final double maxCheckInterval, 094 final double convergence, final int maxIterationCount) { 095 this.handler = handler; 096 this.maxCheckInterval = maxCheckInterval; 097 this.convergence = Math.abs(convergence); 098 this.maxIterationCount = maxIterationCount; 099 100 // some dummy values ... 101 t0 = Double.NaN; 102 g0 = Double.NaN; 103 g0Positive = true; 104 pendingEvent = false; 105 pendingEventTime = Double.NaN; 106 previousEventTime = Double.NaN; 107 increasing = true; 108 nextAction = EventHandler.CONTINUE; 109 110 } 111 112 /** Get the underlying event handler. 113 * @return underlying event handler 114 */ 115 public EventHandler getEventHandler() { 116 return handler; 117 } 118 119 /** Get the maximal time interval between events handler checks. 120 * @return maximal time interval between events handler checks 121 */ 122 public double getMaxCheckInterval() { 123 return maxCheckInterval; 124 } 125 126 /** Get the convergence threshold for event localization. 127 * @return convergence threshold for event localization 128 */ 129 public double getConvergence() { 130 return convergence; 131 } 132 133 /** Get the upper limit in the iteration count for event localization. 134 * @return upper limit in the iteration count for event localization 135 */ 136 public int getMaxIterationCount() { 137 return maxIterationCount; 138 } 139 140 /** Reinitialize the beginning of the step. 141 * @param t0 value of the independent <i>time</i> variable at the 142 * beginning of the step 143 * @param y0 array containing the current value of the state vector 144 * at the beginning of the step 145 * @exception EventException if the event handler 146 * value cannot be evaluated at the beginning of the step 147 */ 148 public void reinitializeBegin(final double t0, final double[] y0) 149 throws EventException { 150 this.t0 = t0; 151 g0 = handler.g(t0, y0); 152 g0Positive = (g0 >= 0); 153 } 154 155 /** Evaluate the impact of the proposed step on the event handler. 156 * @param interpolator step interpolator for the proposed step 157 * @return true if the event handler triggers an event before 158 * the end of the proposed step (this implies the step should be 159 * rejected) 160 * @exception DerivativeException if the interpolator fails to 161 * compute the switching function somewhere within the step 162 * @exception EventException if the switching function 163 * cannot be evaluated 164 * @exception ConvergenceException if an event cannot be located 165 */ 166 public boolean evaluateStep(final StepInterpolator interpolator) 167 throws DerivativeException, EventException, ConvergenceException { 168 169 try { 170 171 forward = interpolator.isForward(); 172 final double t1 = interpolator.getCurrentTime(); 173 final int n = Math.max(1, (int) Math.ceil(Math.abs(t1 - t0) / maxCheckInterval)); 174 final double h = (t1 - t0) / n; 175 176 double ta = t0; 177 double ga = g0; 178 double tb = t0 + (interpolator.isForward() ? convergence : -convergence); 179 for (int i = 0; i < n; ++i) { 180 181 // evaluate handler value at the end of the substep 182 tb += h; 183 interpolator.setInterpolatedTime(tb); 184 final double gb = handler.g(tb, interpolator.getInterpolatedState()); 185 186 // check events occurrence 187 if (g0Positive ^ (gb >= 0)) { 188 // there is a sign change: an event is expected during this step 189 190 // variation direction, with respect to the integration direction 191 increasing = (gb >= ga); 192 193 final UnivariateRealFunction f = new UnivariateRealFunction() { 194 public double value(final double t) throws FunctionEvaluationException { 195 try { 196 interpolator.setInterpolatedTime(t); 197 return handler.g(t, interpolator.getInterpolatedState()); 198 } catch (DerivativeException e) { 199 throw new FunctionEvaluationException(e, t); 200 } catch (EventException e) { 201 throw new FunctionEvaluationException(e, t); 202 } 203 } 204 }; 205 final BrentSolver solver = new BrentSolver(); 206 solver.setAbsoluteAccuracy(convergence); 207 solver.setMaximalIterationCount(maxIterationCount); 208 double root; 209 try { 210 root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta); 211 } catch (IllegalArgumentException iae) { 212 // the interval did not really bracket a root 213 root = Double.NaN; 214 } 215 if (Double.isNaN(root) || 216 ((Math.abs(root - ta) <= convergence) && 217 (Math.abs(root - previousEventTime) <= convergence))) { 218 // we have either found nothing or found (again ?) a past event, we simply ignore it 219 ta = tb; 220 ga = gb; 221 } else if (Double.isNaN(previousEventTime) || 222 (Math.abs(previousEventTime - root) > convergence)) { 223 pendingEventTime = root; 224 if (pendingEvent && (Math.abs(t1 - pendingEventTime) <= convergence)) { 225 // we were already waiting for this event which was 226 // found during a previous call for a step that was 227 // rejected, this step must now be accepted since it 228 // properly ends exactly at the event occurrence 229 return false; 230 } 231 // either we were not waiting for the event or it has 232 // moved in such a way the step cannot be accepted 233 pendingEvent = true; 234 return true; 235 } 236 237 } else { 238 // no sign change: there is no event for now 239 ta = tb; 240 ga = gb; 241 } 242 243 } 244 245 // no event during the whole step 246 pendingEvent = false; 247 pendingEventTime = Double.NaN; 248 return false; 249 250 } catch (FunctionEvaluationException e) { 251 final Throwable cause = e.getCause(); 252 if ((cause != null) && (cause instanceof DerivativeException)) { 253 throw (DerivativeException) cause; 254 } else if ((cause != null) && (cause instanceof EventException)) { 255 throw (EventException) cause; 256 } 257 throw new EventException(e); 258 } 259 260 } 261 262 /** Get the occurrence time of the event triggered in the current 263 * step. 264 * @return occurrence time of the event triggered in the current 265 * step. 266 */ 267 public double getEventTime() { 268 return pendingEventTime; 269 } 270 271 /** Acknowledge the fact the step has been accepted by the integrator. 272 * @param t value of the independent <i>time</i> variable at the 273 * end of the step 274 * @param y array containing the current value of the state vector 275 * at the end of the step 276 * @exception EventException if the value of the event 277 * handler cannot be evaluated 278 */ 279 public void stepAccepted(final double t, final double[] y) 280 throws EventException { 281 282 t0 = t; 283 g0 = handler.g(t, y); 284 285 if (pendingEvent) { 286 // force the sign to its value "just after the event" 287 previousEventTime = t; 288 g0Positive = increasing; 289 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward)); 290 } else { 291 g0Positive = (g0 >= 0); 292 nextAction = EventHandler.CONTINUE; 293 } 294 } 295 296 /** Check if the integration should be stopped at the end of the 297 * current step. 298 * @return true if the integration should be stopped 299 */ 300 public boolean stop() { 301 return nextAction == EventHandler.STOP; 302 } 303 304 /** Let the event handler reset the state if it wants. 305 * @param t value of the independent <i>time</i> variable at the 306 * beginning of the next step 307 * @param y array were to put the desired state vector at the beginning 308 * of the next step 309 * @return true if the integrator should reset the derivatives too 310 * @exception EventException if the state cannot be reseted by the event 311 * handler 312 */ 313 public boolean reset(final double t, final double[] y) 314 throws EventException { 315 316 if (! pendingEvent) { 317 return false; 318 } 319 320 if (nextAction == EventHandler.RESET_STATE) { 321 handler.resetState(t, y); 322 } 323 pendingEvent = false; 324 pendingEventTime = Double.NaN; 325 326 return (nextAction == EventHandler.RESET_STATE) || 327 (nextAction == EventHandler.RESET_DERIVATIVES); 328 329 } 330 331 }