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 018package org.apache.commons.math3.ode.events; 019 020import org.apache.commons.math3.RealFieldElement; 021import org.apache.commons.math3.analysis.RealFieldUnivariateFunction; 022import org.apache.commons.math3.analysis.solvers.AllowedSolution; 023import org.apache.commons.math3.analysis.solvers.BracketedRealFieldUnivariateSolver; 024import org.apache.commons.math3.exception.MaxCountExceededException; 025import org.apache.commons.math3.exception.NoBracketingException; 026import org.apache.commons.math3.ode.FieldODEState; 027import org.apache.commons.math3.ode.FieldODEStateAndDerivative; 028import org.apache.commons.math3.ode.sampling.FieldStepInterpolator; 029import org.apache.commons.math3.util.FastMath; 030 031/** This class handles the state for one {@link EventHandler 032 * event handler} during integration steps. 033 * 034 * <p>Each time the integrator proposes a step, the event handler 035 * switching function should be checked. This class handles the state 036 * of one handler during one integration step, with references to the 037 * state at the end of the preceding step. This information is used to 038 * decide if the handler should trigger an event or not during the 039 * proposed step.</p> 040 * 041 * @param <T> the type of the field elements 042 * @since 3.6 043 */ 044public class FieldEventState<T extends RealFieldElement<T>> { 045 046 /** Event handler. */ 047 private final FieldEventHandler<T> handler; 048 049 /** Maximal time interval between events handler checks. */ 050 private final double maxCheckInterval; 051 052 /** Convergence threshold for event localization. */ 053 private final T convergence; 054 055 /** Upper limit in the iteration count for event localization. */ 056 private final int maxIterationCount; 057 058 /** Time at the beginning of the step. */ 059 private T t0; 060 061 /** Value of the events handler at the beginning of the step. */ 062 private T g0; 063 064 /** Simulated sign of g0 (we cheat when crossing events). */ 065 private boolean g0Positive; 066 067 /** Indicator of event expected during the step. */ 068 private boolean pendingEvent; 069 070 /** Occurrence time of the pending event. */ 071 private T pendingEventTime; 072 073 /** Occurrence time of the previous event. */ 074 private T previousEventTime; 075 076 /** Integration direction. */ 077 private boolean forward; 078 079 /** Variation direction around pending event. 080 * (this is considered with respect to the integration direction) 081 */ 082 private boolean increasing; 083 084 /** Next action indicator. */ 085 private Action nextAction; 086 087 /** Root-finding algorithm to use to detect state events. */ 088 private final BracketedRealFieldUnivariateSolver<T> solver; 089 090 /** Simple constructor. 091 * @param handler event handler 092 * @param maxCheckInterval maximal time interval between switching 093 * function checks (this interval prevents missing sign changes in 094 * case the integration steps becomes very large) 095 * @param convergence convergence threshold in the event time search 096 * @param maxIterationCount upper limit of the iteration count in 097 * the event time search 098 * @param solver Root-finding algorithm to use to detect state events 099 */ 100 public FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval, 101 final T convergence, final int maxIterationCount, 102 final BracketedRealFieldUnivariateSolver<T> solver) { 103 this.handler = handler; 104 this.maxCheckInterval = maxCheckInterval; 105 this.convergence = convergence.abs(); 106 this.maxIterationCount = maxIterationCount; 107 this.solver = solver; 108 109 // some dummy values ... 110 t0 = null; 111 g0 = null; 112 g0Positive = true; 113 pendingEvent = false; 114 pendingEventTime = null; 115 previousEventTime = null; 116 increasing = true; 117 nextAction = Action.CONTINUE; 118 119 } 120 121 /** Get the underlying event handler. 122 * @return underlying event handler 123 */ 124 public FieldEventHandler<T> getEventHandler() { 125 return handler; 126 } 127 128 /** Get the maximal time interval between events handler checks. 129 * @return maximal time interval between events handler checks 130 */ 131 public double getMaxCheckInterval() { 132 return maxCheckInterval; 133 } 134 135 /** Get the convergence threshold for event localization. 136 * @return convergence threshold for event localization 137 */ 138 public T getConvergence() { 139 return convergence; 140 } 141 142 /** Get the upper limit in the iteration count for event localization. 143 * @return upper limit in the iteration count for event localization 144 */ 145 public int getMaxIterationCount() { 146 return maxIterationCount; 147 } 148 149 /** Reinitialize the beginning of the step. 150 * @param interpolator valid for the current step 151 * @exception MaxCountExceededException if the interpolator throws one because 152 * the number of functions evaluations is exceeded 153 */ 154 public void reinitializeBegin(final FieldStepInterpolator<T> interpolator) 155 throws MaxCountExceededException { 156 157 final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState(); 158 t0 = s0.getTime(); 159 g0 = handler.g(s0); 160 if (g0.getReal() == 0) { 161 // excerpt from MATH-421 issue: 162 // If an ODE solver is setup with an EventHandler that return STOP 163 // when the even is triggered, the integrator stops (which is exactly 164 // the expected behavior). If however the user wants to restart the 165 // solver from the final state reached at the event with the same 166 // configuration (expecting the event to be triggered again at a 167 // later time), then the integrator may fail to start. It can get stuck 168 // at the previous event. The use case for the bug MATH-421 is fairly 169 // general, so events occurring exactly at start in the first step should 170 // be ignored. 171 172 // extremely rare case: there is a zero EXACTLY at interval start 173 // we will use the sign slightly after step beginning to force ignoring this zero 174 final double epsilon = FastMath.max(solver.getAbsoluteAccuracy().getReal(), 175 FastMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal())); 176 final T tStart = t0.add(0.5 * epsilon); 177 g0 = handler.g(interpolator.getInterpolatedState(tStart)); 178 } 179 g0Positive = g0.getReal() >= 0; 180 181 } 182 183 /** Evaluate the impact of the proposed step on the event handler. 184 * @param interpolator step interpolator for the proposed step 185 * @return true if the event handler triggers an event before 186 * the end of the proposed step 187 * @exception MaxCountExceededException if the interpolator throws one because 188 * the number of functions evaluations is exceeded 189 * @exception NoBracketingException if the event cannot be bracketed 190 */ 191 public boolean evaluateStep(final FieldStepInterpolator<T> interpolator) 192 throws MaxCountExceededException, NoBracketingException { 193 194 forward = interpolator.isForward(); 195 final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState(); 196 final T t1 = s1.getTime(); 197 final T dt = t1.subtract(t0); 198 if (dt.abs().subtract(convergence).getReal() < 0) { 199 // we cannot do anything on such a small step, don't trigger any events 200 return false; 201 } 202 final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt.getReal()) / maxCheckInterval)); 203 final T h = dt.divide(n); 204 205 final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() { 206 /** {@inheritDoc} */ 207 public T value(final T t) { 208 return handler.g(interpolator.getInterpolatedState(t)); 209 } 210 }; 211 212 T ta = t0; 213 T ga = g0; 214 for (int i = 0; i < n; ++i) { 215 216 // evaluate handler value at the end of the substep 217 final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1)); 218 final T gb = handler.g(interpolator.getInterpolatedState(tb)); 219 220 // check events occurrence 221 if (g0Positive ^ (gb.getReal() >= 0)) { 222 // there is a sign change: an event is expected during this step 223 224 // variation direction, with respect to the integration direction 225 increasing = gb.subtract(ga).getReal() >= 0; 226 227 // find the event time making sure we select a solution just at or past the exact root 228 final T root = forward ? 229 solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) : 230 solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE); 231 232 if (previousEventTime != null && 233 root.subtract(ta).abs().subtract(convergence).getReal() <= 0 && 234 root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) { 235 // we have either found nothing or found (again ?) a past event, 236 // retry the substep excluding this value, and taking care to have the 237 // required sign in case the g function is noisy around its zero and 238 // crosses the axis several times 239 do { 240 ta = forward ? ta.add(convergence) : ta.subtract(convergence); 241 ga = f.value(ta); 242 } while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0))); 243 244 if (forward ^ (ta.subtract(tb).getReal() >= 0)) { 245 // we were able to skip this spurious root 246 --i; 247 } else { 248 // we can't avoid this root before the end of the step, 249 // we have to handle it despite it is close to the former one 250 // maybe we have two very close roots 251 pendingEventTime = root; 252 pendingEvent = true; 253 return true; 254 } 255 } else if (previousEventTime == null || 256 previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) { 257 pendingEventTime = root; 258 pendingEvent = true; 259 return true; 260 } else { 261 // no sign change: there is no event for now 262 ta = tb; 263 ga = gb; 264 } 265 266 } else { 267 // no sign change: there is no event for now 268 ta = tb; 269 ga = gb; 270 } 271 272 } 273 274 // no event during the whole step 275 pendingEvent = false; 276 pendingEventTime = null; 277 return false; 278 279 } 280 281 /** Get the occurrence time of the event triggered in the current step. 282 * @return occurrence time of the event triggered in the current 283 * step or infinity if no events are triggered 284 */ 285 public T getEventTime() { 286 return pendingEvent ? 287 pendingEventTime : 288 t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY); 289 } 290 291 /** Acknowledge the fact the step has been accepted by the integrator. 292 * @param state state at the end of the step 293 */ 294 public void stepAccepted(final FieldODEStateAndDerivative<T> state) { 295 296 t0 = state.getTime(); 297 g0 = handler.g(state); 298 299 if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) { 300 // force the sign to its value "just after the event" 301 previousEventTime = state.getTime(); 302 g0Positive = increasing; 303 nextAction = handler.eventOccurred(state, !(increasing ^ forward)); 304 } else { 305 g0Positive = g0.getReal() >= 0; 306 nextAction = Action.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 == Action.STOP; 316 } 317 318 /** Let the event handler reset the state if it wants. 319 * @param state state at the beginning of the next step 320 * @return reset state (may by the same as initial state if only 321 * derivatives should be reset), or null if nothing is reset 322 */ 323 public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) { 324 325 if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) { 326 return null; 327 } 328 329 final FieldODEState<T> newState; 330 if (nextAction == Action.RESET_STATE) { 331 newState = handler.resetState(state); 332 } else if (nextAction == Action.RESET_DERIVATIVES) { 333 newState = state; 334 } else { 335 newState = null; 336 } 337 pendingEvent = false; 338 pendingEventTime = null; 339 340 return newState; 341 342 } 343 344}