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