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;
019
020 import java.util.ArrayList;
021 import java.util.Collection;
022 import java.util.Collections;
023 import java.util.Comparator;
024 import java.util.Iterator;
025 import java.util.List;
026 import java.util.SortedSet;
027 import java.util.TreeSet;
028
029 import org.apache.commons.math.ConvergenceException;
030 import org.apache.commons.math.MaxEvaluationsExceededException;
031 import org.apache.commons.math.exception.util.LocalizedFormats;
032 import org.apache.commons.math.ode.events.CombinedEventsManager;
033 import org.apache.commons.math.ode.events.EventException;
034 import org.apache.commons.math.ode.events.EventHandler;
035 import org.apache.commons.math.ode.events.EventState;
036 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
037 import org.apache.commons.math.ode.sampling.StepHandler;
038 import org.apache.commons.math.util.FastMath;
039 import org.apache.commons.math.util.MathUtils;
040
041 /**
042 * Base class managing common boilerplate for all integrators.
043 * @version $Revision: 1073267 $ $Date: 2011-02-22 10:06:20 +0100 (mar. 22 f??vr. 2011) $
044 * @since 2.0
045 */
046 public abstract class AbstractIntegrator implements FirstOrderIntegrator {
047
048 /** Step handler. */
049 protected Collection<StepHandler> stepHandlers;
050
051 /** Current step start time. */
052 protected double stepStart;
053
054 /** Current stepsize. */
055 protected double stepSize;
056
057 /** Indicator for last step. */
058 protected boolean isLastStep;
059
060 /** Indicator that a state or derivative reset was triggered by some event. */
061 protected boolean resetOccurred;
062
063 /** Events states. */
064 private Collection<EventState> eventsStates;
065
066 /** Initialization indicator of events states. */
067 private boolean statesInitialized;
068
069 /** Name of the method. */
070 private final String name;
071
072 /** Maximal number of evaluations allowed. */
073 private int maxEvaluations;
074
075 /** Number of evaluations already performed. */
076 private int evaluations;
077
078 /** Differential equations to integrate. */
079 private transient FirstOrderDifferentialEquations equations;
080
081 /** Build an instance.
082 * @param name name of the method
083 */
084 public AbstractIntegrator(final String name) {
085 this.name = name;
086 stepHandlers = new ArrayList<StepHandler>();
087 stepStart = Double.NaN;
088 stepSize = Double.NaN;
089 eventsStates = new ArrayList<EventState>();
090 statesInitialized = false;
091 setMaxEvaluations(-1);
092 resetEvaluations();
093 }
094
095 /** Build an instance with a null name.
096 */
097 protected AbstractIntegrator() {
098 this(null);
099 }
100
101 /** {@inheritDoc} */
102 public String getName() {
103 return name;
104 }
105
106 /** {@inheritDoc} */
107 public void addStepHandler(final StepHandler handler) {
108 stepHandlers.add(handler);
109 }
110
111 /** {@inheritDoc} */
112 public Collection<StepHandler> getStepHandlers() {
113 return Collections.unmodifiableCollection(stepHandlers);
114 }
115
116 /** {@inheritDoc} */
117 public void clearStepHandlers() {
118 stepHandlers.clear();
119 }
120
121 /** {@inheritDoc} */
122 public void addEventHandler(final EventHandler handler,
123 final double maxCheckInterval,
124 final double convergence,
125 final int maxIterationCount) {
126 eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount));
127 }
128
129 /** {@inheritDoc} */
130 public Collection<EventHandler> getEventHandlers() {
131 final List<EventHandler> list = new ArrayList<EventHandler>();
132 for (EventState state : eventsStates) {
133 list.add(state.getEventHandler());
134 }
135 return Collections.unmodifiableCollection(list);
136 }
137
138 /** {@inheritDoc} */
139 public void clearEventHandlers() {
140 eventsStates.clear();
141 }
142
143 /** Check if dense output is needed.
144 * @return true if there is at least one event handler or if
145 * one of the step handlers requires dense output
146 */
147 protected boolean requiresDenseOutput() {
148 if (!eventsStates.isEmpty()) {
149 return true;
150 }
151 for (StepHandler handler : stepHandlers) {
152 if (handler.requiresDenseOutput()) {
153 return true;
154 }
155 }
156 return false;
157 }
158
159 /** {@inheritDoc} */
160 public double getCurrentStepStart() {
161 return stepStart;
162 }
163
164 /** {@inheritDoc} */
165 public double getCurrentSignedStepsize() {
166 return stepSize;
167 }
168
169 /** {@inheritDoc} */
170 public void setMaxEvaluations(int maxEvaluations) {
171 this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
172 }
173
174 /** {@inheritDoc} */
175 public int getMaxEvaluations() {
176 return maxEvaluations;
177 }
178
179 /** {@inheritDoc} */
180 public int getEvaluations() {
181 return evaluations;
182 }
183
184 /** Reset the number of evaluations to zero.
185 */
186 protected void resetEvaluations() {
187 evaluations = 0;
188 }
189
190 /** Set the differential equations.
191 * @param equations differential equations to integrate
192 * @see #computeDerivatives(double, double[], double[])
193 */
194 protected void setEquations(final FirstOrderDifferentialEquations equations) {
195 this.equations = equations;
196 }
197
198 /** Compute the derivatives and check the number of evaluations.
199 * @param t current value of the independent <I>time</I> variable
200 * @param y array containing the current value of the state vector
201 * @param yDot placeholder array where to put the time derivative of the state vector
202 * @throws DerivativeException this user-defined exception should be used if an error is
203 * is triggered by user code
204 */
205 public void computeDerivatives(final double t, final double[] y, final double[] yDot)
206 throws DerivativeException {
207 if (++evaluations > maxEvaluations) {
208 throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
209 }
210 equations.computeDerivatives(t, y, yDot);
211 }
212
213 /** Set the stateInitialized flag.
214 * <p>This method must be called by integrators with the value
215 * {@code false} before they start integration, so a proper lazy
216 * initialization is done automatically on the first step.</p>
217 * @param stateInitialized new value for the flag
218 * @since 2.2
219 */
220 protected void setStateInitialized(final boolean stateInitialized) {
221 this.statesInitialized = stateInitialized;
222 }
223
224 /** Accept a step, triggering events and step handlers.
225 * @param interpolator step interpolator
226 * @param y state vector at step end time, must be reset if an event
227 * asks for resetting or if an events stops integration during the step
228 * @param yDot placeholder array where to put the time derivative of the state vector
229 * @param tEnd final integration time
230 * @return time at end of step
231 * @throws DerivativeException this exception is propagated to the caller if
232 * the underlying user function triggers one
233 * @exception IntegratorException if the value of one event state cannot be evaluated
234 * @since 2.2
235 */
236 protected double acceptStep(final AbstractStepInterpolator interpolator,
237 final double[] y, final double[] yDot, final double tEnd)
238 throws DerivativeException, IntegratorException {
239
240 try {
241 double previousT = interpolator.getGlobalPreviousTime();
242 final double currentT = interpolator.getGlobalCurrentTime();
243 resetOccurred = false;
244
245 // initialize the events states if needed
246 if (! statesInitialized) {
247 for (EventState state : eventsStates) {
248 state.reinitializeBegin(interpolator);
249 }
250 statesInitialized = true;
251 }
252
253 // search for next events that may occur during the step
254 final int orderingSign = interpolator.isForward() ? +1 : -1;
255 SortedSet<EventState> occuringEvents = new TreeSet<EventState>(new Comparator<EventState>() {
256
257 /** {@inheritDoc} */
258 public int compare(EventState es0, EventState es1) {
259 return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
260 }
261
262 });
263
264 for (final EventState state : eventsStates) {
265 if (state.evaluateStep(interpolator)) {
266 // the event occurs during the current step
267 occuringEvents.add(state);
268 }
269 }
270
271 while (!occuringEvents.isEmpty()) {
272
273 // handle the chronologically first event
274 final Iterator<EventState> iterator = occuringEvents.iterator();
275 final EventState currentEvent = iterator.next();
276 iterator.remove();
277
278 // restrict the interpolator to the first part of the step, up to the event
279 final double eventT = currentEvent.getEventTime();
280 interpolator.setSoftPreviousTime(previousT);
281 interpolator.setSoftCurrentTime(eventT);
282
283 // trigger the event
284 interpolator.setInterpolatedTime(eventT);
285 final double[] eventY = interpolator.getInterpolatedState();
286 currentEvent.stepAccepted(eventT, eventY);
287 isLastStep = currentEvent.stop();
288
289 // handle the first part of the step, up to the event
290 for (final StepHandler handler : stepHandlers) {
291 handler.handleStep(interpolator, isLastStep);
292 }
293
294 if (isLastStep) {
295 // the event asked to stop integration
296 System.arraycopy(eventY, 0, y, 0, y.length);
297 return eventT;
298 }
299
300 if (currentEvent.reset(eventT, eventY)) {
301 // some event handler has triggered changes that
302 // invalidate the derivatives, we need to recompute them
303 System.arraycopy(eventY, 0, y, 0, y.length);
304 computeDerivatives(eventT, y, yDot);
305 resetOccurred = true;
306 return eventT;
307 }
308
309 // prepare handling of the remaining part of the step
310 previousT = eventT;
311 interpolator.setSoftPreviousTime(eventT);
312 interpolator.setSoftCurrentTime(currentT);
313
314 // check if the same event occurs again in the remaining part of the step
315 if (currentEvent.evaluateStep(interpolator)) {
316 // the event occurs during the current step
317 occuringEvents.add(currentEvent);
318 }
319
320 }
321
322 interpolator.setInterpolatedTime(currentT);
323 final double[] currentY = interpolator.getInterpolatedState();
324 for (final EventState state : eventsStates) {
325 state.stepAccepted(currentT, currentY);
326 isLastStep = isLastStep || state.stop();
327 }
328 isLastStep = isLastStep || MathUtils.equals(currentT, tEnd, 1);
329
330 // handle the remaining part of the step, after all events if any
331 for (StepHandler handler : stepHandlers) {
332 handler.handleStep(interpolator, isLastStep);
333 }
334
335 return currentT;
336 } catch (EventException se) {
337 final Throwable cause = se.getCause();
338 if ((cause != null) && (cause instanceof DerivativeException)) {
339 throw (DerivativeException) cause;
340 }
341 throw new IntegratorException(se);
342 } catch (ConvergenceException ce) {
343 throw new IntegratorException(ce);
344 }
345
346 }
347
348 /** Perform some sanity checks on the integration parameters.
349 * @param ode differential equations set
350 * @param t0 start time
351 * @param y0 state vector at t0
352 * @param t target time for the integration
353 * @param y placeholder where to put the state vector
354 * @exception IntegratorException if some inconsistency is detected
355 */
356 protected void sanityChecks(final FirstOrderDifferentialEquations ode,
357 final double t0, final double[] y0,
358 final double t, final double[] y)
359 throws IntegratorException {
360
361 if (ode.getDimension() != y0.length) {
362 throw new IntegratorException(
363 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y0.length);
364 }
365
366 if (ode.getDimension() != y.length) {
367 throw new IntegratorException(
368 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y.length);
369 }
370
371 if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) {
372 throw new IntegratorException(
373 LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
374 FastMath.abs(t - t0));
375 }
376
377 }
378
379 /** Add an event handler for end time checking.
380 * <p>This method can be used to simplify handling of integration end time.
381 * It leverages the nominal stop condition with the exceptional stop
382 * conditions.</p>
383 * @param startTime integration start time
384 * @param endTime desired end time
385 * @param manager manager containing the user-defined handlers
386 * @return a new manager containing all the user-defined handlers plus a
387 * dedicated manager triggering a stop event at entTime
388 * @deprecated as of 2.2, this method is not used any more
389 */
390 @Deprecated
391 protected CombinedEventsManager addEndTimeChecker(final double startTime,
392 final double endTime,
393 final CombinedEventsManager manager) {
394 CombinedEventsManager newManager = new CombinedEventsManager();
395 for (final EventState state : manager.getEventsStates()) {
396 newManager.addEventHandler(state.getEventHandler(),
397 state.getMaxCheckInterval(),
398 state.getConvergence(),
399 state.getMaxIterationCount());
400 }
401 newManager.addEventHandler(new EndTimeChecker(endTime),
402 Double.POSITIVE_INFINITY,
403 FastMath.ulp(FastMath.max(FastMath.abs(startTime), FastMath.abs(endTime))),
404 100);
405 return newManager;
406 }
407
408 /** Specialized event handler to stop integration.
409 * @deprecated as of 2.2, this class is not used anymore
410 */
411 @Deprecated
412 private static class EndTimeChecker implements EventHandler {
413
414 /** Desired end time. */
415 private final double endTime;
416
417 /** Build an instance.
418 * @param endTime desired time
419 */
420 public EndTimeChecker(final double endTime) {
421 this.endTime = endTime;
422 }
423
424 /** {@inheritDoc} */
425 public int eventOccurred(double t, double[] y, boolean increasing) {
426 return STOP;
427 }
428
429 /** {@inheritDoc} */
430 public double g(double t, double[] y) {
431 return t - endTime;
432 }
433
434 /** {@inheritDoc} */
435 public void resetState(double t, double[] y) {
436 }
437
438 }
439
440 }