org.apache.commons.math3.ode.AbstractFieldIntegrator Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of virtdata-lib-realer Show documentation
Show all versions of virtdata-lib-realer Show documentation
With inspiration from other libraries
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.ode;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.Comparator;
import java.util.Iterator;
import java.util.List;
import java.util.SortedSet;
import java.util.TreeSet;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.analysis.solvers.BracketedRealFieldUnivariateSolver;
import org.apache.commons.math3.analysis.solvers.FieldBracketingNthOrderBrentSolver;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.ode.events.FieldEventHandler;
import org.apache.commons.math3.ode.events.FieldEventState;
import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator;
import org.apache.commons.math3.ode.sampling.FieldStepHandler;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.IntegerSequence;
/**
* Base class managing common boilerplate for all integrators.
* @param the type of the field elements
* @since 3.6
*/
public abstract class AbstractFieldIntegrator> implements FirstOrderFieldIntegrator {
/** Default relative accuracy. */
private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
/** Default function value accuracy. */
private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
/** Step handler. */
private Collection> stepHandlers;
/** Current step start. */
private FieldODEStateAndDerivative stepStart;
/** Current stepsize. */
private T stepSize;
/** Indicator for last step. */
private boolean isLastStep;
/** Indicator that a state or derivative reset was triggered by some event. */
private boolean resetOccurred;
/** Field to which the time and state vector elements belong. */
private final Field field;
/** Events states. */
private Collection> eventsStates;
/** Initialization indicator of events states. */
private boolean statesInitialized;
/** Name of the method. */
private final String name;
/** Counter for number of evaluations. */
private IntegerSequence.Incrementor evaluations;
/** Differential equations to integrate. */
private transient FieldExpandableODE equations;
/** Build an instance.
* @param field field to which the time and state vector elements belong
* @param name name of the method
*/
protected AbstractFieldIntegrator(final Field field, final String name) {
this.field = field;
this.name = name;
stepHandlers = new ArrayList>();
stepStart = null;
stepSize = null;
eventsStates = new ArrayList>();
statesInitialized = false;
evaluations = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
}
/** Get the field to which state vector elements belong.
* @return field to which state vector elements belong
*/
public Field getField() {
return field;
}
/** {@inheritDoc} */
public String getName() {
return name;
}
/** {@inheritDoc} */
public void addStepHandler(final FieldStepHandler handler) {
stepHandlers.add(handler);
}
/** {@inheritDoc} */
public Collection> getStepHandlers() {
return Collections.unmodifiableCollection(stepHandlers);
}
/** {@inheritDoc} */
public void clearStepHandlers() {
stepHandlers.clear();
}
/** {@inheritDoc} */
public void addEventHandler(final FieldEventHandler handler,
final double maxCheckInterval,
final double convergence,
final int maxIterationCount) {
addEventHandler(handler, maxCheckInterval, convergence,
maxIterationCount,
new FieldBracketingNthOrderBrentSolver(field.getZero().add(DEFAULT_RELATIVE_ACCURACY),
field.getZero().add(convergence),
field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY),
5));
}
/** {@inheritDoc} */
public void addEventHandler(final FieldEventHandler handler,
final double maxCheckInterval,
final double convergence,
final int maxIterationCount,
final BracketedRealFieldUnivariateSolver solver) {
eventsStates.add(new FieldEventState(handler, maxCheckInterval, field.getZero().add(convergence),
maxIterationCount, solver));
}
/** {@inheritDoc} */
public Collection> getEventHandlers() {
final List> list = new ArrayList>(eventsStates.size());
for (FieldEventState state : eventsStates) {
list.add(state.getEventHandler());
}
return Collections.unmodifiableCollection(list);
}
/** {@inheritDoc} */
public void clearEventHandlers() {
eventsStates.clear();
}
/** {@inheritDoc} */
public FieldODEStateAndDerivative getCurrentStepStart() {
return stepStart;
}
/** {@inheritDoc} */
public T getCurrentSignedStepsize() {
return stepSize;
}
/** {@inheritDoc} */
public void setMaxEvaluations(int maxEvaluations) {
evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
}
/** {@inheritDoc} */
public int getMaxEvaluations() {
return evaluations.getMaximalCount();
}
/** {@inheritDoc} */
public int getEvaluations() {
return evaluations.getCount();
}
/** Prepare the start of an integration.
* @param eqn equations to integrate
* @param t0 start value of the independent time variable
* @param y0 array containing the start value of the state vector
* @param t target time for the integration
* @return initial state with derivatives added
*/
protected FieldODEStateAndDerivative initIntegration(final FieldExpandableODE eqn,
final T t0, final T[] y0, final T t) {
this.equations = eqn;
evaluations = evaluations.withStart(0);
// initialize ODE
eqn.init(t0, y0, t);
// set up derivatives of initial state
final T[] y0Dot = computeDerivatives(t0, y0);
final FieldODEStateAndDerivative state0 = new FieldODEStateAndDerivative(t0, y0, y0Dot);
// initialize event handlers
for (final FieldEventState state : eventsStates) {
state.getEventHandler().init(state0, t);
}
// initialize step handlers
for (FieldStepHandler handler : stepHandlers) {
handler.init(state0, t);
}
setStateInitialized(false);
return state0;
}
/** Get the differential equations to integrate.
* @return differential equations to integrate
*/
protected FieldExpandableODE getEquations() {
return equations;
}
/** Get the evaluations counter.
* @return evaluations counter
*/
protected IntegerSequence.Incrementor getEvaluationsCounter() {
return evaluations;
}
/** Compute the derivatives and check the number of evaluations.
* @param t current value of the independent time variable
* @param y array containing the current value of the state vector
* @return state completed with derivatives
* @exception DimensionMismatchException if arrays dimensions do not match equations settings
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
* @exception NullPointerException if the ODE equations have not been set (i.e. if this method
* is called outside of a call to {@link #integrate(FieldExpandableODE, FieldODEState,
* RealFieldElement) integrate}
*/
public T[] computeDerivatives(final T t, final T[] y)
throws DimensionMismatchException, MaxCountExceededException, NullPointerException {
evaluations.increment();
return equations.computeDerivatives(t, y);
}
/** Set the stateInitialized flag.
* This method must be called by integrators with the value
* {@code false} before they start integration, so a proper lazy
* initialization is done automatically on the first step.
* @param stateInitialized new value for the flag
*/
protected void setStateInitialized(final boolean stateInitialized) {
this.statesInitialized = stateInitialized;
}
/** Accept a step, triggering events and step handlers.
* @param interpolator step interpolator
* @param tEnd final integration time
* @return state at end of step
* @exception MaxCountExceededException if the interpolator throws one because
* the number of functions evaluations is exceeded
* @exception NoBracketingException if the location of an event cannot be bracketed
* @exception DimensionMismatchException if arrays dimensions do not match equations settings
*/
protected FieldODEStateAndDerivative acceptStep(final AbstractFieldStepInterpolator interpolator,
final T tEnd)
throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
FieldODEStateAndDerivative previousState = interpolator.getGlobalPreviousState();
final FieldODEStateAndDerivative currentState = interpolator.getGlobalCurrentState();
// initialize the events states if needed
if (! statesInitialized) {
for (FieldEventState state : eventsStates) {
state.reinitializeBegin(interpolator);
}
statesInitialized = true;
}
// search for next events that may occur during the step
final int orderingSign = interpolator.isForward() ? +1 : -1;
SortedSet> occurringEvents = new TreeSet>(new Comparator>() {
/** {@inheritDoc} */
public int compare(FieldEventState es0, FieldEventState es1) {
return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal());
}
});
for (final FieldEventState state : eventsStates) {
if (state.evaluateStep(interpolator)) {
// the event occurs during the current step
occurringEvents.add(state);
}
}
AbstractFieldStepInterpolator restricted = interpolator;
while (!occurringEvents.isEmpty()) {
// handle the chronologically first event
final Iterator> iterator = occurringEvents.iterator();
final FieldEventState currentEvent = iterator.next();
iterator.remove();
// get state at event time
final FieldODEStateAndDerivative eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
// restrict the interpolator to the first part of the step, up to the event
restricted = restricted.restrictStep(previousState, eventState);
// advance all event states to current time
for (final FieldEventState state : eventsStates) {
state.stepAccepted(eventState);
isLastStep = isLastStep || state.stop();
}
// handle the first part of the step, up to the event
for (final FieldStepHandler handler : stepHandlers) {
handler.handleStep(restricted, isLastStep);
}
if (isLastStep) {
// the event asked to stop integration
return eventState;
}
FieldODEState newState = null;
resetOccurred = false;
for (final FieldEventState state : eventsStates) {
newState = state.reset(eventState);
if (newState != null) {
// some event handler has triggered changes that
// invalidate the derivatives, we need to recompute them
final T[] y = equations.getMapper().mapState(newState);
final T[] yDot = computeDerivatives(newState.getTime(), y);
resetOccurred = true;
return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot);
}
}
// prepare handling of the remaining part of the step
previousState = eventState;
restricted = restricted.restrictStep(eventState, currentState);
// check if the same event occurs again in the remaining part of the step
if (currentEvent.evaluateStep(restricted)) {
// the event occurs during the current step
occurringEvents.add(currentEvent);
}
}
// last part of the step, after the last event
for (final FieldEventState state : eventsStates) {
state.stepAccepted(currentState);
isLastStep = isLastStep || state.stop();
}
isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= FastMath.ulp(tEnd.getReal());
// handle the remaining part of the step, after all events if any
for (FieldStepHandler handler : stepHandlers) {
handler.handleStep(restricted, isLastStep);
}
return currentState;
}
/** Check the integration span.
* @param eqn set of differential equations
* @param t target time for the integration
* @exception NumberIsTooSmallException if integration span is too small
* @exception DimensionMismatchException if adaptive step size integrators
* tolerance arrays dimensions are not compatible with equations settings
*/
protected void sanityChecks(final FieldODEState eqn, final T t)
throws NumberIsTooSmallException, DimensionMismatchException {
final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(eqn.getTime().getReal()),
FastMath.abs(t.getReal())));
final double dt = eqn.getTime().subtract(t).abs().getReal();
if (dt <= threshold) {
throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
dt, threshold, false);
}
}
/** Check if a reset occurred while last step was accepted.
* @return true if a reset occurred while last step was accepted
*/
protected boolean resetOccurred() {
return resetOccurred;
}
/** Set the current step size.
* @param stepSize step size to set
*/
protected void setStepSize(final T stepSize) {
this.stepSize = stepSize;
}
/** Get the current step size.
* @return current step size
*/
protected T getStepSize() {
return stepSize;
}
/** Set current step start.
* @param stepStart step start
*/
protected void setStepStart(final FieldODEStateAndDerivative stepStart) {
this.stepStart = stepStart;
}
/** Getcurrent step start.
* @return current step start
*/
protected FieldODEStateAndDerivative getStepStart() {
return stepStart;
}
/** Set the last state flag.
* @param isLastStep if true, this step is the last one
*/
protected void setIsLastStep(final boolean isLastStep) {
this.isLastStep = isLastStep;
}
/** Check if this step is the last one.
* @return true if this step is the last one
*/
protected boolean isLastStep() {
return isLastStep;
}
}