1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math.ode.events;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.FunctionEvaluationException;
22 import org.apache.commons.math.analysis.UnivariateRealFunction;
23 import org.apache.commons.math.analysis.solvers.BrentSolver;
24 import org.apache.commons.math.ode.DerivativeException;
25 import org.apache.commons.math.ode.sampling.StepInterpolator;
26
27 /** This class handles the state for one {@link EventHandler
28 * event handler} during integration steps.
29 *
30 * <p>Each time the integrator proposes a step, the event handler
31 * switching function should be checked. This class handles the state
32 * of one handler during one integration step, with references to the
33 * state at the end of the preceding step. This information is used to
34 * decide if the handler should trigger an event or not during the
35 * proposed step (and hence the step should be reduced to ensure the
36 * event occurs at a bound rather than inside the step).</p>
37 *
38 * @version $Revision: 786881 $ $Date: 2009-06-20 14:53:08 -0400 (Sat, 20 Jun 2009) $
39 * @since 1.2
40 */
41 public class EventState {
42
43 /** Event handler. */
44 private final EventHandler handler;
45
46 /** Maximal time interval between events handler checks. */
47 private final double maxCheckInterval;
48
49 /** Convergence threshold for event localization. */
50 private final double convergence;
51
52 /** Upper limit in the iteration count for event localization. */
53 private final int maxIterationCount;
54
55 /** Time at the beginning of the step. */
56 private double t0;
57
58 /** Value of the events handler at the beginning of the step. */
59 private double g0;
60
61 /** Simulated sign of g0 (we cheat when crossing events). */
62 private boolean g0Positive;
63
64 /** Indicator of event expected during the step. */
65 private boolean pendingEvent;
66
67 /** Occurrence time of the pending event. */
68 private double pendingEventTime;
69
70 /** Occurrence time of the previous event. */
71 private double previousEventTime;
72
73 /** Integration direction. */
74 private boolean forward;
75
76 /** Variation direction around pending event.
77 * (this is considered with respect to the integration direction)
78 */
79 private boolean increasing;
80
81 /** Next action indicator. */
82 private int nextAction;
83
84 /** Simple constructor.
85 * @param handler event handler
86 * @param maxCheckInterval maximal time interval between switching
87 * function checks (this interval prevents missing sign changes in
88 * case the integration steps becomes very large)
89 * @param convergence convergence threshold in the event time search
90 * @param maxIterationCount upper limit of the iteration count in
91 * the event time search
92 */
93 public EventState(final EventHandler handler, final double maxCheckInterval,
94 final double convergence, final int maxIterationCount) {
95 this.handler = handler;
96 this.maxCheckInterval = maxCheckInterval;
97 this.convergence = Math.abs(convergence);
98 this.maxIterationCount = maxIterationCount;
99
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 }