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 }