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.Field; 030import org.apache.commons.math3.RealFieldElement; 031import org.apache.commons.math3.analysis.solvers.BracketedRealFieldUnivariateSolver; 032import org.apache.commons.math3.analysis.solvers.FieldBracketingNthOrderBrentSolver; 033import org.apache.commons.math3.exception.DimensionMismatchException; 034import org.apache.commons.math3.exception.MaxCountExceededException; 035import org.apache.commons.math3.exception.NoBracketingException; 036import org.apache.commons.math3.exception.NumberIsTooSmallException; 037import org.apache.commons.math3.exception.util.LocalizedFormats; 038import org.apache.commons.math3.ode.events.FieldEventHandler; 039import org.apache.commons.math3.ode.events.FieldEventState; 040import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator; 041import org.apache.commons.math3.ode.sampling.FieldStepHandler; 042import org.apache.commons.math3.util.FastMath; 043import org.apache.commons.math3.util.IntegerSequence; 044 045/** 046 * Base class managing common boilerplate for all integrators. 047 * @param <T> the type of the field elements 048 * @since 3.6 049 */ 050public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FirstOrderFieldIntegrator<T> { 051 052 /** Default relative accuracy. */ 053 private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14; 054 055 /** Default function value accuracy. */ 056 private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15; 057 058 /** Step handler. */ 059 private Collection<FieldStepHandler<T>> stepHandlers; 060 061 /** Current step start. */ 062 private FieldODEStateAndDerivative<T> stepStart; 063 064 /** Current stepsize. */ 065 private T stepSize; 066 067 /** Indicator for last step. */ 068 private boolean isLastStep; 069 070 /** Indicator that a state or derivative reset was triggered by some event. */ 071 private boolean resetOccurred; 072 073 /** Field to which the time and state vector elements belong. */ 074 private final Field<T> field; 075 076 /** Events states. */ 077 private Collection<FieldEventState<T>> eventsStates; 078 079 /** Initialization indicator of events states. */ 080 private boolean statesInitialized; 081 082 /** Name of the method. */ 083 private final String name; 084 085 /** Counter for number of evaluations. */ 086 private IntegerSequence.Incrementor evaluations; 087 088 /** Differential equations to integrate. */ 089 private transient FieldExpandableODE<T> equations; 090 091 /** Build an instance. 092 * @param field field to which the time and state vector elements belong 093 * @param name name of the method 094 */ 095 protected AbstractFieldIntegrator(final Field<T> field, final String name) { 096 this.field = field; 097 this.name = name; 098 stepHandlers = new ArrayList<FieldStepHandler<T>>(); 099 stepStart = null; 100 stepSize = null; 101 eventsStates = new ArrayList<FieldEventState<T>>(); 102 statesInitialized = false; 103 evaluations = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE); 104 } 105 106 /** Get the field to which state vector elements belong. 107 * @return field to which state vector elements belong 108 */ 109 public Field<T> getField() { 110 return field; 111 } 112 113 /** {@inheritDoc} */ 114 public String getName() { 115 return name; 116 } 117 118 /** {@inheritDoc} */ 119 public void addStepHandler(final FieldStepHandler<T> handler) { 120 stepHandlers.add(handler); 121 } 122 123 /** {@inheritDoc} */ 124 public Collection<FieldStepHandler<T>> getStepHandlers() { 125 return Collections.unmodifiableCollection(stepHandlers); 126 } 127 128 /** {@inheritDoc} */ 129 public void clearStepHandlers() { 130 stepHandlers.clear(); 131 } 132 133 /** {@inheritDoc} */ 134 public void addEventHandler(final FieldEventHandler<T> handler, 135 final double maxCheckInterval, 136 final double convergence, 137 final int maxIterationCount) { 138 addEventHandler(handler, maxCheckInterval, convergence, 139 maxIterationCount, 140 new FieldBracketingNthOrderBrentSolver<T>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY), 141 field.getZero().add(convergence), 142 field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY), 143 5)); 144 } 145 146 /** {@inheritDoc} */ 147 public void addEventHandler(final FieldEventHandler<T> handler, 148 final double maxCheckInterval, 149 final double convergence, 150 final int maxIterationCount, 151 final BracketedRealFieldUnivariateSolver<T> solver) { 152 eventsStates.add(new FieldEventState<T>(handler, maxCheckInterval, field.getZero().add(convergence), 153 maxIterationCount, solver)); 154 } 155 156 /** {@inheritDoc} */ 157 public Collection<FieldEventHandler<T>> getEventHandlers() { 158 final List<FieldEventHandler<T>> list = new ArrayList<FieldEventHandler<T>>(eventsStates.size()); 159 for (FieldEventState<T> state : eventsStates) { 160 list.add(state.getEventHandler()); 161 } 162 return Collections.unmodifiableCollection(list); 163 } 164 165 /** {@inheritDoc} */ 166 public void clearEventHandlers() { 167 eventsStates.clear(); 168 } 169 170 /** {@inheritDoc} */ 171 public FieldODEStateAndDerivative<T> getCurrentStepStart() { 172 return stepStart; 173 } 174 175 /** {@inheritDoc} */ 176 public T getCurrentSignedStepsize() { 177 return stepSize; 178 } 179 180 /** {@inheritDoc} */ 181 public void setMaxEvaluations(int maxEvaluations) { 182 evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations); 183 } 184 185 /** {@inheritDoc} */ 186 public int getMaxEvaluations() { 187 return evaluations.getMaximalCount(); 188 } 189 190 /** {@inheritDoc} */ 191 public int getEvaluations() { 192 return evaluations.getCount(); 193 } 194 195 /** Prepare the start of an integration. 196 * @param eqn equations to integrate 197 * @param t0 start value of the independent <i>time</i> variable 198 * @param y0 array containing the start value of the state vector 199 * @param t target time for the integration 200 * @return initial state with derivatives added 201 */ 202 protected FieldODEStateAndDerivative<T> initIntegration(final FieldExpandableODE<T> eqn, 203 final T t0, final T[] y0, final T t) { 204 205 this.equations = eqn; 206 evaluations = evaluations.withStart(0); 207 208 // initialize ODE 209 eqn.init(t0, y0, t); 210 211 // set up derivatives of initial state 212 final T[] y0Dot = computeDerivatives(t0, y0); 213 final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<T>(t0, y0, y0Dot); 214 215 // initialize event handlers 216 for (final FieldEventState<T> state : eventsStates) { 217 state.getEventHandler().init(state0, t); 218 } 219 220 // initialize step handlers 221 for (FieldStepHandler<T> handler : stepHandlers) { 222 handler.init(state0, t); 223 } 224 225 setStateInitialized(false); 226 227 return state0; 228 229 } 230 231 /** Get the differential equations to integrate. 232 * @return differential equations to integrate 233 */ 234 protected FieldExpandableODE<T> getEquations() { 235 return equations; 236 } 237 238 /** Get the evaluations counter. 239 * @return evaluations counter 240 */ 241 protected IntegerSequence.Incrementor getEvaluationsCounter() { 242 return evaluations; 243 } 244 245 /** Compute the derivatives and check the number of evaluations. 246 * @param t current value of the independent <I>time</I> variable 247 * @param y array containing the current value of the state vector 248 * @return state completed with derivatives 249 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 250 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 251 * @exception NullPointerException if the ODE equations have not been set (i.e. if this method 252 * is called outside of a call to {@link #integrate(FieldExpandableODE, FieldODEState, 253 * RealFieldElement) integrate} 254 */ 255 public T[] computeDerivatives(final T t, final T[] y) 256 throws DimensionMismatchException, MaxCountExceededException, NullPointerException { 257 evaluations.increment(); 258 return equations.computeDerivatives(t, y); 259 } 260 261 /** Set the stateInitialized flag. 262 * <p>This method must be called by integrators with the value 263 * {@code false} before they start integration, so a proper lazy 264 * initialization is done automatically on the first step.</p> 265 * @param stateInitialized new value for the flag 266 */ 267 protected void setStateInitialized(final boolean stateInitialized) { 268 this.statesInitialized = stateInitialized; 269 } 270 271 /** Accept a step, triggering events and step handlers. 272 * @param interpolator step interpolator 273 * @param tEnd final integration time 274 * @return state at end of step 275 * @exception MaxCountExceededException if the interpolator throws one because 276 * the number of functions evaluations is exceeded 277 * @exception NoBracketingException if the location of an event cannot be bracketed 278 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 279 */ 280 protected FieldODEStateAndDerivative<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator, 281 final T tEnd) 282 throws MaxCountExceededException, DimensionMismatchException, NoBracketingException { 283 284 FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState(); 285 final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState(); 286 287 // initialize the events states if needed 288 if (! statesInitialized) { 289 for (FieldEventState<T> state : eventsStates) { 290 state.reinitializeBegin(interpolator); 291 } 292 statesInitialized = true; 293 } 294 295 // search for next events that may occur during the step 296 final int orderingSign = interpolator.isForward() ? +1 : -1; 297 SortedSet<FieldEventState<T>> occurringEvents = new TreeSet<FieldEventState<T>>(new Comparator<FieldEventState<T>>() { 298 299 /** {@inheritDoc} */ 300 public int compare(FieldEventState<T> es0, FieldEventState<T> es1) { 301 return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal()); 302 } 303 304 }); 305 306 for (final FieldEventState<T> state : eventsStates) { 307 if (state.evaluateStep(interpolator)) { 308 // the event occurs during the current step 309 occurringEvents.add(state); 310 } 311 } 312 313 AbstractFieldStepInterpolator<T> restricted = interpolator; 314 while (!occurringEvents.isEmpty()) { 315 316 // handle the chronologically first event 317 final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator(); 318 final FieldEventState<T> currentEvent = iterator.next(); 319 iterator.remove(); 320 321 // get state at event time 322 final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime()); 323 324 // restrict the interpolator to the first part of the step, up to the event 325 restricted = restricted.restrictStep(previousState, eventState); 326 327 // advance all event states to current time 328 for (final FieldEventState<T> state : eventsStates) { 329 state.stepAccepted(eventState); 330 isLastStep = isLastStep || state.stop(); 331 } 332 333 // handle the first part of the step, up to the event 334 for (final FieldStepHandler<T> handler : stepHandlers) { 335 handler.handleStep(restricted, isLastStep); 336 } 337 338 if (isLastStep) { 339 // the event asked to stop integration 340 return eventState; 341 } 342 343 FieldODEState<T> newState = null; 344 resetOccurred = false; 345 for (final FieldEventState<T> state : eventsStates) { 346 newState = state.reset(eventState); 347 if (newState != null) { 348 // some event handler has triggered changes that 349 // invalidate the derivatives, we need to recompute them 350 final T[] y = equations.getMapper().mapState(newState); 351 final T[] yDot = computeDerivatives(newState.getTime(), y); 352 resetOccurred = true; 353 return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot); 354 } 355 } 356 357 // prepare handling of the remaining part of the step 358 previousState = eventState; 359 restricted = restricted.restrictStep(eventState, currentState); 360 361 // check if the same event occurs again in the remaining part of the step 362 if (currentEvent.evaluateStep(restricted)) { 363 // the event occurs during the current step 364 occurringEvents.add(currentEvent); 365 } 366 367 } 368 369 // last part of the step, after the last event 370 for (final FieldEventState<T> state : eventsStates) { 371 state.stepAccepted(currentState); 372 isLastStep = isLastStep || state.stop(); 373 } 374 isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= FastMath.ulp(tEnd.getReal()); 375 376 // handle the remaining part of the step, after all events if any 377 for (FieldStepHandler<T> handler : stepHandlers) { 378 handler.handleStep(restricted, isLastStep); 379 } 380 381 return currentState; 382 383 } 384 385 /** Check the integration span. 386 * @param eqn set of differential equations 387 * @param t target time for the integration 388 * @exception NumberIsTooSmallException if integration span is too small 389 * @exception DimensionMismatchException if adaptive step size integrators 390 * tolerance arrays dimensions are not compatible with equations settings 391 */ 392 protected void sanityChecks(final FieldODEState<T> eqn, final T t) 393 throws NumberIsTooSmallException, DimensionMismatchException { 394 395 final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(eqn.getTime().getReal()), 396 FastMath.abs(t.getReal()))); 397 final double dt = eqn.getTime().subtract(t).abs().getReal(); 398 if (dt <= threshold) { 399 throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL, 400 dt, threshold, false); 401 } 402 403 } 404 405 /** Check if a reset occurred while last step was accepted. 406 * @return true if a reset occurred while last step was accepted 407 */ 408 protected boolean resetOccurred() { 409 return resetOccurred; 410 } 411 412 /** Set the current step size. 413 * @param stepSize step size to set 414 */ 415 protected void setStepSize(final T stepSize) { 416 this.stepSize = stepSize; 417 } 418 419 /** Get the current step size. 420 * @return current step size 421 */ 422 protected T getStepSize() { 423 return stepSize; 424 } 425 /** Set current step start. 426 * @param stepStart step start 427 */ 428 protected void setStepStart(final FieldODEStateAndDerivative<T> stepStart) { 429 this.stepStart = stepStart; 430 } 431 432 /** Getcurrent step start. 433 * @return current step start 434 */ 435 protected FieldODEStateAndDerivative<T> getStepStart() { 436 return stepStart; 437 } 438 439 /** Set the last state flag. 440 * @param isLastStep if true, this step is the last one 441 */ 442 protected void setIsLastStep(final boolean isLastStep) { 443 this.isLastStep = isLastStep; 444 } 445 446 /** Check if this step is the last one. 447 * @return true if this step is the last one 448 */ 449 protected boolean isLastStep() { 450 return isLastStep; 451 } 452 453}