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 }