View Javadoc

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 }