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.math3.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.math3.analysis.solvers.BracketingNthOrderBrentSolver;
030 import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
031 import org.apache.commons.math3.exception.DimensionMismatchException;
032 import org.apache.commons.math3.exception.MaxCountExceededException;
033 import org.apache.commons.math3.exception.NoBracketingException;
034 import org.apache.commons.math3.exception.NumberIsTooSmallException;
035 import org.apache.commons.math3.exception.util.LocalizedFormats;
036 import org.apache.commons.math3.ode.events.EventHandler;
037 import org.apache.commons.math3.ode.events.EventState;
038 import org.apache.commons.math3.ode.sampling.AbstractStepInterpolator;
039 import org.apache.commons.math3.ode.sampling.StepHandler;
040 import org.apache.commons.math3.util.FastMath;
041 import org.apache.commons.math3.util.Incrementor;
042 import org.apache.commons.math3.util.Precision;
043
044 /**
045 * Base class managing common boilerplate for all integrators.
046 * @version $Id: AbstractIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
047 * @since 2.0
048 */
049 public abstract class AbstractIntegrator implements FirstOrderIntegrator {
050
051 /** Step handler. */
052 protected Collection<StepHandler> stepHandlers;
053
054 /** Current step start time. */
055 protected double stepStart;
056
057 /** Current stepsize. */
058 protected double stepSize;
059
060 /** Indicator for last step. */
061 protected boolean isLastStep;
062
063 /** Indicator that a state or derivative reset was triggered by some event. */
064 protected boolean resetOccurred;
065
066 /** Events states. */
067 private Collection<EventState> eventsStates;
068
069 /** Initialization indicator of events states. */
070 private boolean statesInitialized;
071
072 /** Name of the method. */
073 private final String name;
074
075 /** Counter for number of evaluations. */
076 private Incrementor evaluations;
077
078 /** Differential equations to integrate. */
079 private transient ExpandableStatefulODE expandable;
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 evaluations = new Incrementor();
092 setMaxEvaluations(-1);
093 evaluations.resetCount();
094 }
095
096 /** Build an instance with a null name.
097 */
098 protected AbstractIntegrator() {
099 this(null);
100 }
101
102 /** {@inheritDoc} */
103 public String getName() {
104 return name;
105 }
106
107 /** {@inheritDoc} */
108 public void addStepHandler(final StepHandler handler) {
109 stepHandlers.add(handler);
110 }
111
112 /** {@inheritDoc} */
113 public Collection<StepHandler> getStepHandlers() {
114 return Collections.unmodifiableCollection(stepHandlers);
115 }
116
117 /** {@inheritDoc} */
118 public void clearStepHandlers() {
119 stepHandlers.clear();
120 }
121
122 /** {@inheritDoc} */
123 public void addEventHandler(final EventHandler handler,
124 final double maxCheckInterval,
125 final double convergence,
126 final int maxIterationCount) {
127 addEventHandler(handler, maxCheckInterval, convergence,
128 maxIterationCount,
129 new BracketingNthOrderBrentSolver(convergence, 5));
130 }
131
132 /** {@inheritDoc} */
133 public void addEventHandler(final EventHandler handler,
134 final double maxCheckInterval,
135 final double convergence,
136 final int maxIterationCount,
137 final UnivariateSolver solver) {
138 eventsStates.add(new EventState(handler, maxCheckInterval, convergence,
139 maxIterationCount, solver));
140 }
141
142 /** {@inheritDoc} */
143 public Collection<EventHandler> getEventHandlers() {
144 final List<EventHandler> list = new ArrayList<EventHandler>();
145 for (EventState state : eventsStates) {
146 list.add(state.getEventHandler());
147 }
148 return Collections.unmodifiableCollection(list);
149 }
150
151 /** {@inheritDoc} */
152 public void clearEventHandlers() {
153 eventsStates.clear();
154 }
155
156 /** {@inheritDoc} */
157 public double getCurrentStepStart() {
158 return stepStart;
159 }
160
161 /** {@inheritDoc} */
162 public double getCurrentSignedStepsize() {
163 return stepSize;
164 }
165
166 /** {@inheritDoc} */
167 public void setMaxEvaluations(int maxEvaluations) {
168 evaluations.setMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
169 }
170
171 /** {@inheritDoc} */
172 public int getMaxEvaluations() {
173 return evaluations.getMaximalCount();
174 }
175
176 /** {@inheritDoc} */
177 public int getEvaluations() {
178 return evaluations.getCount();
179 }
180
181 /** Prepare the start of an integration.
182 * @param t0 start value of the independent <i>time</i> variable
183 * @param y0 array containing the start value of the state vector
184 * @param t target time for the integration
185 */
186 protected void initIntegration(final double t0, final double[] y0, final double t) {
187
188 evaluations.resetCount();
189
190 for (final EventState state : eventsStates) {
191 state.getEventHandler().init(t0, y0, t);
192 }
193
194 for (StepHandler handler : stepHandlers) {
195 handler.init(t0, y0, t);
196 }
197
198 setStateInitialized(false);
199
200 }
201
202 /** Set the equations.
203 * @param equations equations to set
204 */
205 protected void setEquations(final ExpandableStatefulODE equations) {
206 this.expandable = equations;
207 }
208
209 /** {@inheritDoc} */
210 public double integrate(final FirstOrderDifferentialEquations equations,
211 final double t0, final double[] y0, final double t, final double[] y)
212 throws DimensionMismatchException, NumberIsTooSmallException,
213 MaxCountExceededException, NoBracketingException {
214
215 if (y0.length != equations.getDimension()) {
216 throw new DimensionMismatchException(y0.length, equations.getDimension());
217 }
218 if (y.length != equations.getDimension()) {
219 throw new DimensionMismatchException(y.length, equations.getDimension());
220 }
221
222 // prepare expandable stateful equations
223 final ExpandableStatefulODE expandableODE = new ExpandableStatefulODE(equations);
224 expandableODE.setTime(t0);
225 expandableODE.setPrimaryState(y0);
226
227 // perform integration
228 integrate(expandableODE, t);
229
230 // extract results back from the stateful equations
231 System.arraycopy(expandableODE.getPrimaryState(), 0, y, 0, y.length);
232 return expandableODE.getTime();
233
234 }
235
236 /** Integrate a set of differential equations up to the given time.
237 * <p>This method solves an Initial Value Problem (IVP).</p>
238 * <p>The set of differential equations is composed of a main set, which
239 * can be extended by some sets of secondary equations. The set of
240 * equations must be already set up with initial time and partial states.
241 * At integration completion, the final time and partial states will be
242 * available in the same object.</p>
243 * <p>Since this method stores some internal state variables made
244 * available in its public interface during integration ({@link
245 * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
246 * @param equations complete set of differential equations to integrate
247 * @param t target time for the integration
248 * (can be set to a value smaller than <code>t0</code> for backward integration)
249 * @exception NumberIsTooSmallException if integration step is too small
250 * @throws DimensionMismatchException if the dimension of the complete state does not
251 * match the complete equations sets dimension
252 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
253 * @exception NoBracketingException if the location of an event cannot be bracketed
254 */
255 public abstract void integrate(ExpandableStatefulODE equations, double t)
256 throws NumberIsTooSmallException, DimensionMismatchException,
257 MaxCountExceededException, NoBracketingException;
258
259 /** Compute the derivatives and check the number of evaluations.
260 * @param t current value of the independent <I>time</I> variable
261 * @param y array containing the current value of the state vector
262 * @param yDot placeholder array where to put the time derivative of the state vector
263 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
264 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
265 */
266 public void computeDerivatives(final double t, final double[] y, final double[] yDot)
267 throws MaxCountExceededException, DimensionMismatchException {
268 evaluations.incrementCount();
269 expandable.computeDerivatives(t, y, yDot);
270 }
271
272 /** Set the stateInitialized flag.
273 * <p>This method must be called by integrators with the value
274 * {@code false} before they start integration, so a proper lazy
275 * initialization is done automatically on the first step.</p>
276 * @param stateInitialized new value for the flag
277 * @since 2.2
278 */
279 protected void setStateInitialized(final boolean stateInitialized) {
280 this.statesInitialized = stateInitialized;
281 }
282
283 /** Accept a step, triggering events and step handlers.
284 * @param interpolator step interpolator
285 * @param y state vector at step end time, must be reset if an event
286 * asks for resetting or if an events stops integration during the step
287 * @param yDot placeholder array where to put the time derivative of the state vector
288 * @param tEnd final integration time
289 * @return time at end of step
290 * @exception MaxCountExceededException if the interpolator throws one because
291 * the number of functions evaluations is exceeded
292 * @exception NoBracketingException if the location of an event cannot be bracketed
293 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
294 * @since 2.2
295 */
296 protected double acceptStep(final AbstractStepInterpolator interpolator,
297 final double[] y, final double[] yDot, final double tEnd)
298 throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
299
300 double previousT = interpolator.getGlobalPreviousTime();
301 final double currentT = interpolator.getGlobalCurrentTime();
302
303 // initialize the events states if needed
304 if (! statesInitialized) {
305 for (EventState state : eventsStates) {
306 state.reinitializeBegin(interpolator);
307 }
308 statesInitialized = true;
309 }
310
311 // search for next events that may occur during the step
312 final int orderingSign = interpolator.isForward() ? +1 : -1;
313 SortedSet<EventState> occuringEvents = new TreeSet<EventState>(new Comparator<EventState>() {
314
315 /** {@inheritDoc} */
316 public int compare(EventState es0, EventState es1) {
317 return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
318 }
319
320 });
321
322 for (final EventState state : eventsStates) {
323 if (state.evaluateStep(interpolator)) {
324 // the event occurs during the current step
325 occuringEvents.add(state);
326 }
327 }
328
329 while (!occuringEvents.isEmpty()) {
330
331 // handle the chronologically first event
332 final Iterator<EventState> iterator = occuringEvents.iterator();
333 final EventState currentEvent = iterator.next();
334 iterator.remove();
335
336 // restrict the interpolator to the first part of the step, up to the event
337 final double eventT = currentEvent.getEventTime();
338 interpolator.setSoftPreviousTime(previousT);
339 interpolator.setSoftCurrentTime(eventT);
340
341 // trigger the event
342 interpolator.setInterpolatedTime(eventT);
343 final double[] eventY = interpolator.getInterpolatedState().clone();
344 currentEvent.stepAccepted(eventT, eventY);
345 isLastStep = currentEvent.stop();
346
347 // handle the first part of the step, up to the event
348 for (final StepHandler handler : stepHandlers) {
349 handler.handleStep(interpolator, isLastStep);
350 }
351
352 if (isLastStep) {
353 // the event asked to stop integration
354 System.arraycopy(eventY, 0, y, 0, y.length);
355 for (final EventState remaining : occuringEvents) {
356 remaining.stepAccepted(eventT, eventY);
357 }
358 return eventT;
359 }
360
361 if (currentEvent.reset(eventT, eventY)) {
362 // some event handler has triggered changes that
363 // invalidate the derivatives, we need to recompute them
364 System.arraycopy(eventY, 0, y, 0, y.length);
365 computeDerivatives(eventT, y, yDot);
366 resetOccurred = true;
367 for (final EventState remaining : occuringEvents) {
368 remaining.stepAccepted(eventT, eventY);
369 }
370 return eventT;
371 }
372
373 // prepare handling of the remaining part of the step
374 previousT = eventT;
375 interpolator.setSoftPreviousTime(eventT);
376 interpolator.setSoftCurrentTime(currentT);
377
378 // check if the same event occurs again in the remaining part of the step
379 if (currentEvent.evaluateStep(interpolator)) {
380 // the event occurs during the current step
381 occuringEvents.add(currentEvent);
382 }
383
384 }
385
386 interpolator.setInterpolatedTime(currentT);
387 final double[] currentY = interpolator.getInterpolatedState();
388 for (final EventState state : eventsStates) {
389 state.stepAccepted(currentT, currentY);
390 isLastStep = isLastStep || state.stop();
391 }
392 isLastStep = isLastStep || Precision.equals(currentT, tEnd, 1);
393
394 // handle the remaining part of the step, after all events if any
395 for (StepHandler handler : stepHandlers) {
396 handler.handleStep(interpolator, isLastStep);
397 }
398
399 return currentT;
400
401 }
402
403 /** Check the integration span.
404 * @param equations set of differential equations
405 * @param t target time for the integration
406 * @exception NumberIsTooSmallException if integration span is too small
407 * @exception DimensionMismatchException if adaptive step size integrators
408 * tolerance arrays dimensions are not compatible with equations settings
409 */
410 protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
411 throws NumberIsTooSmallException, DimensionMismatchException {
412
413 final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
414 FastMath.abs(t)));
415 final double dt = FastMath.abs(equations.getTime() - t);
416 if (dt <= threshold) {
417 throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
418 dt, threshold, false);
419 }
420
421 }
422
423 }