Implementation of the top level abstract classes in field ode.
This layer implements boilerplate code, mainly step handling and events handling. It is independent of the type of integrator used. Below this layer will be the real implementations (Runge-Kutta, embedded runge-Kutta, Gragg-Bulirsch-Stoer, Adams, ...).
This commit is contained in:
Luc Maisonobe
parent
7bf9d3dbc5
commit
6da8a0eba0
|
|
@ -0,0 +1,390 @@
|
|||
/*
|
||||
* 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.math4.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.math4.Field;
|
||||
import org.apache.commons.math4.RealFieldElement;
|
||||
import org.apache.commons.math4.analysis.solvers.FieldBracketingNthOrderBrentSolver;
|
||||
import org.apache.commons.math4.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.exception.NoBracketingException;
|
||||
import org.apache.commons.math4.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math4.exception.util.LocalizedFormats;
|
||||
import org.apache.commons.math4.ode.events.FieldEventHandler;
|
||||
import org.apache.commons.math4.ode.events.FieldEventState;
|
||||
import org.apache.commons.math4.ode.sampling.AbstractFieldStepInterpolator;
|
||||
import org.apache.commons.math4.ode.sampling.FieldStepHandler;
|
||||
import org.apache.commons.math4.util.FastMath;
|
||||
import org.apache.commons.math4.util.IntegerSequence;
|
||||
|
||||
/**
|
||||
* Base class managing common boilerplate for all integrators.
|
||||
* @param <T> the type of the field elements
|
||||
* @since 3.6
|
||||
*/
|
||||
public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FieldFirstOrderIntegrator<T> {
|
||||
|
||||
/** 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. */
|
||||
protected Collection<FieldStepHandler<T>> stepHandlers;
|
||||
|
||||
/** Current step start time. */
|
||||
protected T stepStart;
|
||||
|
||||
/** Current stepsize. */
|
||||
protected T stepSize;
|
||||
|
||||
/** Indicator for last step. */
|
||||
protected boolean isLastStep;
|
||||
|
||||
/** Indicator that a state or derivative reset was triggered by some event. */
|
||||
protected boolean resetOccurred;
|
||||
|
||||
/** Field to which the time and state vector elements belong. */
|
||||
private final Field<T> field;
|
||||
|
||||
/** Events states. */
|
||||
private Collection<FieldEventState<T>> 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<T> 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<T> field, final String name) {
|
||||
this.field = field;
|
||||
this.name = name;
|
||||
stepHandlers = new ArrayList<FieldStepHandler<T>>();
|
||||
stepStart = null;
|
||||
stepSize = null;
|
||||
eventsStates = new ArrayList<FieldEventState<T>>();
|
||||
statesInitialized = false;
|
||||
evaluations = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public String getName() {
|
||||
return name;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void addStepHandler(final FieldStepHandler<T> handler) {
|
||||
stepHandlers.add(handler);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public Collection<FieldStepHandler<T>> getStepHandlers() {
|
||||
return Collections.unmodifiableCollection(stepHandlers);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void clearStepHandlers() {
|
||||
stepHandlers.clear();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void addEventHandler(final FieldEventHandler<T> handler,
|
||||
final double maxCheckInterval,
|
||||
final double convergence,
|
||||
final int maxIterationCount) {
|
||||
addEventHandler(handler, maxCheckInterval, convergence,
|
||||
maxIterationCount,
|
||||
new FieldBracketingNthOrderBrentSolver<T>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY),
|
||||
field.getZero().add(convergence),
|
||||
field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY),
|
||||
5));
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void addEventHandler(final FieldEventHandler<T> handler,
|
||||
final double maxCheckInterval,
|
||||
final double convergence,
|
||||
final int maxIterationCount,
|
||||
final FieldBracketingNthOrderBrentSolver<T> solver) {
|
||||
eventsStates.add(new FieldEventState<T>(handler, maxCheckInterval, field.getZero().add(convergence),
|
||||
maxIterationCount, solver));
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public Collection<FieldEventHandler<T>> getEventHandlers() {
|
||||
final List<FieldEventHandler<T>> list = new ArrayList<FieldEventHandler<T>>(eventsStates.size());
|
||||
for (FieldEventState<T> state : eventsStates) {
|
||||
list.add(state.getEventHandler());
|
||||
}
|
||||
return Collections.unmodifiableCollection(list);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void clearEventHandlers() {
|
||||
eventsStates.clear();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public T 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 <i>time</i> variable
|
||||
* @param y0 array containing the start value of the state vector
|
||||
* @param t target time for the integration
|
||||
* @return derivative of the state at t0 (y0Dot)
|
||||
*/
|
||||
protected T[] initIntegration(final FieldExpandableODE<T> eqn,
|
||||
final T t0, final T[] y0, final T t) {
|
||||
|
||||
this.equations = eqn;
|
||||
evaluations = evaluations.withStart(0);
|
||||
final T[] y0Dot = computeDerivatives(t0, y0);
|
||||
final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<T>(t0, y0, y0Dot);
|
||||
|
||||
for (final FieldEventState<T> state : eventsStates) {
|
||||
state.getEventHandler().init(state0, t);
|
||||
}
|
||||
|
||||
for (FieldStepHandler<T> handler : stepHandlers) {
|
||||
handler.init(state0, t);
|
||||
}
|
||||
|
||||
setStateInitialized(false);
|
||||
|
||||
return y0Dot;
|
||||
|
||||
}
|
||||
|
||||
/** Get the differential equations to integrate.
|
||||
* @return differential equations to integrate
|
||||
*/
|
||||
protected FieldExpandableODE<T> 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 <I>time</I> 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(ExpandableStatefulODE, double)} or {@link
|
||||
* #integrate(FirstOrderDifferentialEquations, double, double[], double, double[])})
|
||||
*/
|
||||
public T[] computeDerivatives(final T t, final T[] y)
|
||||
throws DimensionMismatchException, MaxCountExceededException, NullPointerException {
|
||||
evaluations.increment();
|
||||
return equations.computeDerivatives(t, y);
|
||||
}
|
||||
|
||||
/** Set the stateInitialized flag.
|
||||
* <p>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.</p>
|
||||
* @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<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator,
|
||||
final T tEnd)
|
||||
throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
|
||||
|
||||
FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState();
|
||||
final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState();
|
||||
|
||||
// initialize the events states if needed
|
||||
if (! statesInitialized) {
|
||||
for (FieldEventState<T> 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<FieldEventState<T>> occurringEvents = new TreeSet<FieldEventState<T>>(new Comparator<FieldEventState<T>>() {
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int compare(FieldEventState<T> es0, FieldEventState<T> es1) {
|
||||
return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal());
|
||||
}
|
||||
|
||||
});
|
||||
|
||||
for (final FieldEventState<T> state : eventsStates) {
|
||||
if (state.evaluateStep(interpolator)) {
|
||||
// the event occurs during the current step
|
||||
occurringEvents.add(state);
|
||||
}
|
||||
}
|
||||
|
||||
while (!occurringEvents.isEmpty()) {
|
||||
|
||||
// handle the chronologically first event
|
||||
final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator();
|
||||
final FieldEventState<T> currentEvent = iterator.next();
|
||||
iterator.remove();
|
||||
|
||||
// get state at event time
|
||||
final FieldODEStateAndDerivative<T> eventState = interpolator.getInterpolatedState(currentEvent.getEventTime());
|
||||
|
||||
// restrict the interpolator to the first part of the step, up to the event
|
||||
interpolator.setSoftPreviousState(previousState);
|
||||
interpolator.setSoftCurrentState(eventState);
|
||||
|
||||
// advance all event states to current time
|
||||
for (final FieldEventState<T> state : eventsStates) {
|
||||
state.stepAccepted(eventState);
|
||||
isLastStep = isLastStep || state.stop();
|
||||
}
|
||||
|
||||
// handle the first part of the step, up to the event
|
||||
for (final FieldStepHandler<T> handler : stepHandlers) {
|
||||
handler.handleStep(interpolator, isLastStep);
|
||||
}
|
||||
|
||||
if (isLastStep) {
|
||||
// the event asked to stop integration
|
||||
return eventState;
|
||||
}
|
||||
|
||||
boolean needReset = false;
|
||||
for (final FieldEventState<T> state : eventsStates) {
|
||||
needReset = needReset || state.reset(eventState);
|
||||
}
|
||||
if (needReset) {
|
||||
// some event handler has triggered changes that
|
||||
// invalidate the derivatives, we need to recompute them
|
||||
final T[] y = equations.mapState(eventState);
|
||||
final T[] yDot = computeDerivatives(eventState.getTime(), y);
|
||||
resetOccurred = true;
|
||||
return equations.mapStateAndDerivative(eventState.getTime(), y, yDot);
|
||||
}
|
||||
|
||||
// prepare handling of the remaining part of the step
|
||||
previousState = eventState;
|
||||
interpolator.setSoftPreviousState(eventState);
|
||||
interpolator.setSoftCurrentState(currentState);
|
||||
|
||||
// check if the same event occurs again in the remaining part of the step
|
||||
if (currentEvent.evaluateStep(interpolator)) {
|
||||
// the event occurs during the current step
|
||||
occurringEvents.add(currentEvent);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
// last part of the step, after the last event
|
||||
for (final FieldEventState<T> 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<T> handler : stepHandlers) {
|
||||
handler.handleStep(interpolator, 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<T> 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);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -17,8 +17,6 @@
|
|||
|
||||
package org.apache.commons.math4.ode;
|
||||
|
||||
import java.io.Serializable;
|
||||
|
||||
import org.apache.commons.math4.RealFieldElement;
|
||||
import org.apache.commons.math4.util.MathArrays;
|
||||
|
||||
|
|
@ -32,10 +30,7 @@ import org.apache.commons.math4.util.MathArrays;
|
|||
* @param <T> the type of the field elements
|
||||
* @since 3.6
|
||||
*/
|
||||
class FieldEquationsMapper<T extends RealFieldElement<T>> implements Serializable {
|
||||
|
||||
/** Serializable UID. */
|
||||
private static final long serialVersionUID = 20151111L;
|
||||
public class FieldEquationsMapper<T extends RealFieldElement<T>> {
|
||||
|
||||
/** Index of the first equation element in complete state arrays. */
|
||||
private final int firstIndex;
|
||||
|
|
|
|||
|
|
@ -121,7 +121,7 @@ public class FieldExpandableODE<T extends RealFieldElement<T>> {
|
|||
|
||||
/** Map a state to a complete flat array.
|
||||
* @param state state to map
|
||||
* @return flat array containing the mapped state
|
||||
* @return flat array containing the mapped state, including primary and secondary components
|
||||
*/
|
||||
public T[] mapState(final FieldODEState<T> state) {
|
||||
final T[] y = MathArrays.buildArray(state.getTime().getField(), getTotalDimension());
|
||||
|
|
@ -134,7 +134,7 @@ public class FieldExpandableODE<T extends RealFieldElement<T>> {
|
|||
|
||||
/** Map a state derivative to a complete flat array.
|
||||
* @param state state to map
|
||||
* @return flat array containing the mapped state derivative
|
||||
* @return flat array containing the mapped state derivative, including primary and secondary components
|
||||
*/
|
||||
public T[] mapDerivative(final FieldODEStateAndDerivative<T> state) {
|
||||
final T[] yDot = MathArrays.buildArray(state.getTime().getField(), getTotalDimension());
|
||||
|
|
@ -147,7 +147,7 @@ public class FieldExpandableODE<T extends RealFieldElement<T>> {
|
|||
|
||||
/** Map a flat array to a state.
|
||||
* @param t time
|
||||
* @param y array to map
|
||||
* @param y array to map, including primary and secondary components
|
||||
* @return mapped state
|
||||
*/
|
||||
public FieldODEState<T> mapState(final T t, final T[] y) {
|
||||
|
|
@ -166,8 +166,8 @@ public class FieldExpandableODE<T extends RealFieldElement<T>> {
|
|||
|
||||
/** Map flat arrays to a state and derivative.
|
||||
* @param t time
|
||||
* @param y state array to map
|
||||
* @param yDot state derivative array to map
|
||||
* @param y state array to map, including primary and secondary components
|
||||
* @param yDot state derivative array to map, including primary and secondary components
|
||||
* @return mapped state
|
||||
*/
|
||||
public FieldODEStateAndDerivative<T> mapStateAndDerivative(final T t, final T[] y, final T[] yDot) {
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@ package org.apache.commons.math4.ode;
|
|||
import java.util.Collection;
|
||||
|
||||
import org.apache.commons.math4.RealFieldElement;
|
||||
import org.apache.commons.math4.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math4.analysis.solvers.FieldBracketingNthOrderBrentSolver;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.exception.NoBracketingException;
|
||||
import org.apache.commons.math4.exception.NumberIsTooSmallException;
|
||||
|
|
@ -64,11 +64,29 @@ public interface FieldFirstOrderIntegrator<T extends RealFieldElement<T>> {
|
|||
Collection<FieldStepHandler<T>> getStepHandlers();
|
||||
|
||||
/** Remove all the step handlers that have been added to the integrator.
|
||||
* @see #addStepHandler(StepHandler)
|
||||
* @see #addStepHandler(FieldStepHandler)
|
||||
* @see #getStepHandlers()
|
||||
*/
|
||||
void clearStepHandlers();
|
||||
|
||||
/** Add an event handler to the integrator.
|
||||
* <p>
|
||||
* The default solver is a 5<sup>th</sup> order {@link FieldBracketingNthOrderBrentSolver}.
|
||||
* </p>
|
||||
* @param handler event handler
|
||||
* @param maxCheckInterval maximal time interval between switching
|
||||
* function checks (this interval prevents missing sign changes in
|
||||
* case the integration steps becomes very large)
|
||||
* @param convergence convergence threshold in the event time search
|
||||
* @param maxIterationCount upper limit of the iteration count in
|
||||
* the event time search events.
|
||||
* @see #addEventHandler(FieldEventHandler, double, double, int, FieldBracketingNthOrderBrentSolver)
|
||||
* @see #getEventHandlers()
|
||||
* @see #clearEventHandlers()
|
||||
*/
|
||||
void addEventHandler(FieldEventHandler<T> handler, double maxCheckInterval,
|
||||
double convergence, int maxIterationCount);
|
||||
|
||||
/** Add an event handler to the integrator.
|
||||
* @param handler event handler
|
||||
* @param maxCheckInterval maximal time interval between switching
|
||||
|
|
@ -76,23 +94,25 @@ public interface FieldFirstOrderIntegrator<T extends RealFieldElement<T>> {
|
|||
* case the integration steps becomes very large)
|
||||
* @param convergence convergence threshold in the event time search
|
||||
* @param maxIterationCount upper limit of the iteration count in
|
||||
* the event time search
|
||||
* events.
|
||||
* the event time search events.
|
||||
* @param solver solver to use to locate the event
|
||||
* @see #addEventHandler(FieldEventHandler, double, double, int)
|
||||
* @see #getEventHandlers()
|
||||
* @see #clearEventHandlers()
|
||||
*/
|
||||
void addEventHandler(FieldEventHandler<T> handler, double maxCheckInterval,
|
||||
double convergence, int maxIterationCount);
|
||||
double convergence, int maxIterationCount,
|
||||
FieldBracketingNthOrderBrentSolver<T> solver);
|
||||
|
||||
/** Get all the event handlers that have been added to the integrator.
|
||||
* @return an unmodifiable collection of the added events handlers
|
||||
* @see #addEventHandler(EventHandler, double, double, int)
|
||||
* @see #addEventHandler(FieldEventHandler, double, double, int)
|
||||
* @see #clearEventHandlers()
|
||||
*/
|
||||
Collection<FieldEventHandler<T> > getEventHandlers();
|
||||
|
||||
/** Remove all the event handlers that have been added to the integrator.
|
||||
* @see #addEventHandler(EventHandler, double, double, int)
|
||||
* @see #addEventHandler(FieldEventHandler, double, double, int)
|
||||
* @see #getEventHandlers()
|
||||
*/
|
||||
void clearEventHandlers();
|
||||
|
|
@ -155,14 +175,12 @@ public interface FieldFirstOrderIntegrator<T extends RealFieldElement<T>> {
|
|||
* @return final state, its time will be the same as {@code finalTime} if
|
||||
* integration reached its target, but may be different if some {@link
|
||||
* org.apache.commons.math4.ode.events.FieldEventHandler} stops it at some point.
|
||||
* @exception DimensionMismatchException if arrays dimension do not match equations settings
|
||||
* @exception NumberIsTooSmallException if integration step is too small
|
||||
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
|
||||
* @exception NoBracketingException if the location of an event cannot be bracketed
|
||||
*/
|
||||
FieldODEStateAndDerivative<T> integrate (FieldExpandableODE<T> equations,
|
||||
FieldODEState<T>initialState, T finalTime)
|
||||
throws DimensionMismatchException, NumberIsTooSmallException,
|
||||
MaxCountExceededException, NoBracketingException;
|
||||
FieldODEStateAndDerivative<T> integrate(FieldExpandableODE<T> equations,
|
||||
FieldODEState<T> initialState, T finalTime)
|
||||
throws NumberIsTooSmallException, MaxCountExceededException, NoBracketingException;
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,372 @@
|
|||
/*
|
||||
* 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.math4.ode.events;
|
||||
|
||||
import org.apache.commons.math4.RealFieldElement;
|
||||
import org.apache.commons.math4.analysis.FieldUnivariateFunction;
|
||||
import org.apache.commons.math4.analysis.solvers.AllowedSolution;
|
||||
import org.apache.commons.math4.analysis.solvers.FieldBracketingNthOrderBrentSolver;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.exception.NoBracketingException;
|
||||
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
|
||||
import org.apache.commons.math4.ode.sampling.FieldStepInterpolator;
|
||||
import org.apache.commons.math4.util.FastMath;
|
||||
|
||||
/** This class handles the state for one {@link EventHandler
|
||||
* event handler} during integration steps.
|
||||
*
|
||||
* <p>Each time the integrator proposes a step, the event handler
|
||||
* switching function should be checked. This class handles the state
|
||||
* of one handler during one integration step, with references to the
|
||||
* state at the end of the preceding step. This information is used to
|
||||
* decide if the handler should trigger an event or not during the
|
||||
* proposed step.</p>
|
||||
*
|
||||
* @param <T> the type of the field elements
|
||||
* @since 3.6
|
||||
*/
|
||||
public class FieldEventState<T extends RealFieldElement<T>> {
|
||||
|
||||
/** Event handler. */
|
||||
private final FieldEventHandler<T> handler;
|
||||
|
||||
/** Maximal time interval between events handler checks. */
|
||||
private final double maxCheckInterval;
|
||||
|
||||
/** Convergence threshold for event localization. */
|
||||
private final T convergence;
|
||||
|
||||
/** Upper limit in the iteration count for event localization. */
|
||||
private final int maxIterationCount;
|
||||
|
||||
/** Time at the beginning of the step. */
|
||||
private T t0;
|
||||
|
||||
/** Value of the events handler at the beginning of the step. */
|
||||
private T g0;
|
||||
|
||||
/** Simulated sign of g0 (we cheat when crossing events). */
|
||||
private boolean g0Positive;
|
||||
|
||||
/** Indicator of event expected during the step. */
|
||||
private boolean pendingEvent;
|
||||
|
||||
/** Occurrence time of the pending event. */
|
||||
private T pendingEventTime;
|
||||
|
||||
/** Occurrence time of the previous event. */
|
||||
private T previousEventTime;
|
||||
|
||||
/** Integration direction. */
|
||||
private boolean forward;
|
||||
|
||||
/** Variation direction around pending event.
|
||||
* (this is considered with respect to the integration direction)
|
||||
*/
|
||||
private boolean increasing;
|
||||
|
||||
/** Next action indicator. */
|
||||
private Action nextAction;
|
||||
|
||||
/** Root-finding algorithm to use to detect state events. */
|
||||
private final FieldBracketingNthOrderBrentSolver<T> solver;
|
||||
|
||||
/** Simple constructor.
|
||||
* @param handler event handler
|
||||
* @param maxCheckInterval maximal time interval between switching
|
||||
* function checks (this interval prevents missing sign changes in
|
||||
* case the integration steps becomes very large)
|
||||
* @param convergence convergence threshold in the event time search
|
||||
* @param maxIterationCount upper limit of the iteration count in
|
||||
* the event time search
|
||||
* @param solver Root-finding algorithm to use to detect state events
|
||||
*/
|
||||
public FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval,
|
||||
final T convergence, final int maxIterationCount,
|
||||
final FieldBracketingNthOrderBrentSolver<T> solver) {
|
||||
this.handler = handler;
|
||||
this.maxCheckInterval = maxCheckInterval;
|
||||
this.convergence = convergence.abs();
|
||||
this.maxIterationCount = maxIterationCount;
|
||||
this.solver = solver;
|
||||
|
||||
// some dummy values ...
|
||||
t0 = null;
|
||||
g0 = null;
|
||||
g0Positive = true;
|
||||
pendingEvent = false;
|
||||
pendingEventTime = null;
|
||||
previousEventTime = null;
|
||||
increasing = true;
|
||||
nextAction = Action.CONTINUE;
|
||||
|
||||
}
|
||||
|
||||
/** Get the underlying event handler.
|
||||
* @return underlying event handler
|
||||
*/
|
||||
public FieldEventHandler<T> getEventHandler() {
|
||||
return handler;
|
||||
}
|
||||
|
||||
/** Get the maximal time interval between events handler checks.
|
||||
* @return maximal time interval between events handler checks
|
||||
*/
|
||||
public double getMaxCheckInterval() {
|
||||
return maxCheckInterval;
|
||||
}
|
||||
|
||||
/** Get the convergence threshold for event localization.
|
||||
* @return convergence threshold for event localization
|
||||
*/
|
||||
public T getConvergence() {
|
||||
return convergence;
|
||||
}
|
||||
|
||||
/** Get the upper limit in the iteration count for event localization.
|
||||
* @return upper limit in the iteration count for event localization
|
||||
*/
|
||||
public int getMaxIterationCount() {
|
||||
return maxIterationCount;
|
||||
}
|
||||
|
||||
/** Reinitialize the beginning of the step.
|
||||
* @param interpolator valid for the current step
|
||||
* @exception MaxCountExceededException if the interpolator throws one because
|
||||
* the number of functions evaluations is exceeded
|
||||
*/
|
||||
public void reinitializeBegin(final FieldStepInterpolator<T> interpolator)
|
||||
throws MaxCountExceededException {
|
||||
|
||||
final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
|
||||
t0 = s0.getTime();
|
||||
g0 = handler.g(s0);
|
||||
if (g0.getReal() == 0) {
|
||||
// excerpt from MATH-421 issue:
|
||||
// If an ODE solver is setup with an EventHandler that return STOP
|
||||
// when the even is triggered, the integrator stops (which is exactly
|
||||
// the expected behavior). If however the user wants to restart the
|
||||
// solver from the final state reached at the event with the same
|
||||
// configuration (expecting the event to be triggered again at a
|
||||
// later time), then the integrator may fail to start. It can get stuck
|
||||
// at the previous event. The use case for the bug MATH-421 is fairly
|
||||
// general, so events occurring exactly at start in the first step should
|
||||
// be ignored.
|
||||
|
||||
// extremely rare case: there is a zero EXACTLY at interval start
|
||||
// we will use the sign slightly after step beginning to force ignoring this zero
|
||||
final double epsilon = FastMath.max(solver.getAbsoluteAccuracy().getReal(),
|
||||
FastMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
|
||||
final T tStart = t0.add(0.5 * epsilon);
|
||||
g0 = handler.g(interpolator.getInterpolatedState(tStart));
|
||||
}
|
||||
g0Positive = g0.getReal() >= 0;
|
||||
|
||||
}
|
||||
|
||||
/** Evaluate the impact of the proposed step on the event handler.
|
||||
* @param interpolator step interpolator for the proposed step
|
||||
* @return true if the event handler triggers an event before
|
||||
* the end of the proposed step
|
||||
* @exception MaxCountExceededException if the interpolator throws one because
|
||||
* the number of functions evaluations is exceeded
|
||||
* @exception NoBracketingException if the event cannot be bracketed
|
||||
*/
|
||||
public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
|
||||
throws MaxCountExceededException, NoBracketingException {
|
||||
|
||||
try {
|
||||
forward = interpolator.isForward();
|
||||
final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
|
||||
final T t1 = s1.getTime();
|
||||
final T dt = t1.subtract(t0);
|
||||
if (dt.abs().subtract(convergence).getReal() < 0) {
|
||||
// we cannot do anything on such a small step, don't trigger any events
|
||||
return false;
|
||||
}
|
||||
final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt.getReal()) / maxCheckInterval));
|
||||
final T h = dt.divide(n);
|
||||
|
||||
final FieldUnivariateFunction<T> f = new FieldUnivariateFunction<T>() {
|
||||
/** {@inheritDoc} */
|
||||
public T value(final T t) throws LocalMaxCountExceededException {
|
||||
try {
|
||||
return handler.g(interpolator.getInterpolatedState(t));
|
||||
} catch (MaxCountExceededException mcee) {
|
||||
throw new LocalMaxCountExceededException(mcee);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
T ta = t0;
|
||||
T ga = g0;
|
||||
for (int i = 0; i < n; ++i) {
|
||||
|
||||
// evaluate handler value at the end of the substep
|
||||
final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
|
||||
final T gb = handler.g(interpolator.getInterpolatedState(tb));
|
||||
|
||||
// check events occurrence
|
||||
if (g0Positive ^ (gb.getReal() >= 0)) {
|
||||
// there is a sign change: an event is expected during this step
|
||||
|
||||
// variation direction, with respect to the integration direction
|
||||
increasing = gb.subtract(ga).getReal() >= 0;
|
||||
|
||||
// find the event time making sure we select a solution just at or past the exact root
|
||||
final T root = forward ?
|
||||
solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
|
||||
solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
|
||||
|
||||
if (previousEventTime != null &&
|
||||
root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
|
||||
root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
|
||||
// we have either found nothing or found (again ?) a past event,
|
||||
// retry the substep excluding this value, and taking care to have the
|
||||
// required sign in case the g function is noisy around its zero and
|
||||
// crosses the axis several times
|
||||
do {
|
||||
ta = forward ? ta.add(convergence) : ta.subtract(convergence);
|
||||
ga = f.value(ta);
|
||||
} while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
|
||||
|
||||
if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
|
||||
// we were able to skip this spurious root
|
||||
--i;
|
||||
} else {
|
||||
// we can't avoid this root before the end of the step,
|
||||
// we have to handle it despite it is close to the former one
|
||||
// maybe we have two very close roots
|
||||
pendingEventTime = root;
|
||||
pendingEvent = true;
|
||||
return true;
|
||||
}
|
||||
} else if (previousEventTime == null ||
|
||||
previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
|
||||
pendingEventTime = root;
|
||||
pendingEvent = true;
|
||||
return true;
|
||||
} else {
|
||||
// no sign change: there is no event for now
|
||||
ta = tb;
|
||||
ga = gb;
|
||||
}
|
||||
|
||||
} else {
|
||||
// no sign change: there is no event for now
|
||||
ta = tb;
|
||||
ga = gb;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
// no event during the whole step
|
||||
pendingEvent = false;
|
||||
pendingEventTime = null;
|
||||
return false;
|
||||
|
||||
} catch (LocalMaxCountExceededException lmcee) {
|
||||
throw lmcee.getException();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Get the occurrence time of the event triggered in the current step.
|
||||
* @return occurrence time of the event triggered in the current
|
||||
* step or infinity if no events are triggered
|
||||
*/
|
||||
public T getEventTime() {
|
||||
return pendingEvent ?
|
||||
pendingEventTime :
|
||||
t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
|
||||
}
|
||||
|
||||
/** Acknowledge the fact the step has been accepted by the integrator.
|
||||
* @param state state at the end of the step
|
||||
*/
|
||||
public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
|
||||
|
||||
t0 = state.getTime();
|
||||
g0 = handler.g(state);
|
||||
|
||||
if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
|
||||
// force the sign to its value "just after the event"
|
||||
previousEventTime = state.getTime();
|
||||
g0Positive = increasing;
|
||||
nextAction = handler.eventOccurred(state, !(increasing ^ forward));
|
||||
} else {
|
||||
g0Positive = g0.getReal() >= 0;
|
||||
nextAction = Action.CONTINUE;
|
||||
}
|
||||
}
|
||||
|
||||
/** Check if the integration should be stopped at the end of the
|
||||
* current step.
|
||||
* @return true if the integration should be stopped
|
||||
*/
|
||||
public boolean stop() {
|
||||
return nextAction == Action.STOP;
|
||||
}
|
||||
|
||||
/** Let the event handler reset the state if it wants.
|
||||
* @param state state at the beginning of the next step
|
||||
* @return true if the integrator should reset the derivatives too
|
||||
*/
|
||||
public boolean reset(final FieldODEStateAndDerivative<T> state) {
|
||||
|
||||
if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
if (nextAction == Action.RESET_STATE) {
|
||||
handler.resetState(state);
|
||||
}
|
||||
pendingEvent = false;
|
||||
pendingEventTime = null;
|
||||
|
||||
return (nextAction == Action.RESET_STATE) ||
|
||||
(nextAction == Action.RESET_DERIVATIVES);
|
||||
|
||||
}
|
||||
|
||||
/** Local wrapper to propagate exceptions. */
|
||||
private static class LocalMaxCountExceededException extends RuntimeException {
|
||||
|
||||
/** Serializable UID. */
|
||||
private static final long serialVersionUID = 20151113L;
|
||||
|
||||
/** Wrapped exception. */
|
||||
private final MaxCountExceededException wrapped;
|
||||
|
||||
/** Simple constructor.
|
||||
* @param exception exception to wrap
|
||||
*/
|
||||
LocalMaxCountExceededException(final MaxCountExceededException exception) {
|
||||
wrapped = exception;
|
||||
}
|
||||
|
||||
/** Get the wrapped exception.
|
||||
* @return wrapped exception
|
||||
*/
|
||||
public MaxCountExceededException getException() {
|
||||
return wrapped;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,355 @@
|
|||
/*
|
||||
* 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.math4.ode.sampling;
|
||||
|
||||
import org.apache.commons.math4.RealFieldElement;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.ode.FieldExpandableODE;
|
||||
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
|
||||
|
||||
/** This abstract class represents an interpolator over the last step
|
||||
* during an ODE integration.
|
||||
*
|
||||
* <p>The various ODE integrators provide objects extending this class
|
||||
* to the step handlers. The handlers can use these objects to
|
||||
* retrieve the state vector at intermediate times between the
|
||||
* previous and the current grid points (dense output).</p>
|
||||
*
|
||||
* @see org.apache.commons.math4.ode.FieldFirstOrderIntegrator
|
||||
* @see StepHandler
|
||||
*
|
||||
* @param <T> the type of the field elements
|
||||
* @since 3.6
|
||||
*/
|
||||
|
||||
public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T>>
|
||||
implements FieldStepInterpolator<T> {
|
||||
|
||||
/** Current time step. */
|
||||
protected T h;
|
||||
|
||||
/** Current state. */
|
||||
protected T[] currentState;
|
||||
|
||||
/** Global previous state. */
|
||||
private FieldODEStateAndDerivative<T> globalPreviousState;
|
||||
|
||||
/** Global current state. */
|
||||
private FieldODEStateAndDerivative<T> globalCurrentState;
|
||||
|
||||
/** Soft previous state. */
|
||||
private FieldODEStateAndDerivative<T> softPreviousState;
|
||||
|
||||
/** Soft current state. */
|
||||
private FieldODEStateAndDerivative<T> softCurrentState;
|
||||
|
||||
/** indicate if the step has been finalized or not. */
|
||||
private boolean finalized;
|
||||
|
||||
/** integration direction. */
|
||||
private boolean forward;
|
||||
|
||||
/** ODE equations with primary and secondary components. */
|
||||
private FieldExpandableODE<T> equations;
|
||||
|
||||
/** Simple constructor.
|
||||
* This constructor builds an instance that is not usable yet, the
|
||||
* {@link #reinitialize} method should be called before using the
|
||||
* instance in order to initialize the internal arrays. This
|
||||
* constructor is used only in order to delay the initialization in
|
||||
* some cases. As an example, the {@link
|
||||
* org.apache.commons.math4.ode.nonstiff.EmbeddedRungeKuttaIntegrator}
|
||||
* class uses the prototyping design pattern to create the step
|
||||
* interpolators by cloning an uninitialized model and latter
|
||||
* initializing the copy.
|
||||
*/
|
||||
protected AbstractFieldStepInterpolator() {
|
||||
globalPreviousState = null;
|
||||
globalCurrentState = null;
|
||||
softPreviousState = null;
|
||||
softCurrentState = null;
|
||||
h = null;
|
||||
currentState = null;
|
||||
finalized = false;
|
||||
forward = true;
|
||||
equations = null;
|
||||
}
|
||||
|
||||
/** Simple constructor.
|
||||
* @param y reference to the integrator array holding the state at
|
||||
* the end of the step
|
||||
* @param forward integration direction indicator
|
||||
* @param equations primary and secondary equations set
|
||||
*/
|
||||
protected AbstractFieldStepInterpolator(final T[] y, final boolean forward,
|
||||
final FieldExpandableODE<T> equations) {
|
||||
|
||||
globalPreviousState = null;
|
||||
globalCurrentState = null;
|
||||
softPreviousState = null;
|
||||
softCurrentState = null;
|
||||
h = null;
|
||||
currentState = y;
|
||||
finalized = false;
|
||||
this.forward = forward;
|
||||
this.equations = equations;
|
||||
|
||||
}
|
||||
|
||||
/** Copy constructor.
|
||||
|
||||
* <p>The copied interpolator should have been finalized before the
|
||||
* copy, otherwise the copy will not be able to perform correctly
|
||||
* any derivative computation and will throw a {@link
|
||||
* NullPointerException} later. Since we don't want this constructor
|
||||
* to throw the exceptions finalization may involve and since we
|
||||
* don't want this method to modify the state of the copied
|
||||
* interpolator, finalization is <strong>not</strong> done
|
||||
* automatically, it remains under user control.</p>
|
||||
|
||||
* <p>The copy is a deep copy: its arrays are separated from the
|
||||
* original arrays of the instance.</p>
|
||||
|
||||
* @param interpolator interpolator to copy from.
|
||||
|
||||
*/
|
||||
protected AbstractFieldStepInterpolator(final AbstractFieldStepInterpolator<T> interpolator) {
|
||||
|
||||
globalPreviousState = interpolator.globalPreviousState;
|
||||
globalCurrentState = interpolator.globalCurrentState;
|
||||
softPreviousState = interpolator.softPreviousState;
|
||||
softCurrentState = interpolator.softCurrentState;
|
||||
h = interpolator.h;
|
||||
|
||||
if (interpolator.currentState == null) {
|
||||
currentState = null;
|
||||
equations = null;
|
||||
} else {
|
||||
currentState = interpolator.currentState.clone();
|
||||
}
|
||||
|
||||
finalized = interpolator.finalized;
|
||||
forward = interpolator.forward;
|
||||
equations = interpolator.equations;
|
||||
|
||||
}
|
||||
|
||||
/** Reinitialize the instance
|
||||
* @param y reference to the integrator array holding the state at the end of the step
|
||||
* @param isForward integration direction indicator
|
||||
* @param eqn primary and secondary equations set
|
||||
*/
|
||||
protected void reinitialize(final T[] y, final boolean isForward, final FieldExpandableODE<T> eqn) {
|
||||
|
||||
globalPreviousState = null;
|
||||
globalCurrentState = null;
|
||||
softPreviousState = null;
|
||||
softCurrentState = null;
|
||||
h = null;
|
||||
currentState = y.clone();
|
||||
finalized = false;
|
||||
this.forward = isForward;
|
||||
this.equations = eqn;
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public FieldStepInterpolator<T> copy() throws MaxCountExceededException {
|
||||
|
||||
// finalize the step before performing copy
|
||||
finalizeStep();
|
||||
|
||||
// create the new independent instance
|
||||
return doCopy();
|
||||
|
||||
}
|
||||
|
||||
/** Really copy the finalized instance.
|
||||
* <p>This method is called by {@link #copy()} after the
|
||||
* step has been finalized. It must perform a deep copy
|
||||
* to have an new instance completely independent for the
|
||||
* original instance.
|
||||
* @return a copy of the finalized instance
|
||||
*/
|
||||
protected abstract FieldStepInterpolator<T> doCopy();
|
||||
|
||||
/** Shift one step forward.
|
||||
* Copy the current time into the previous time, hence preparing the
|
||||
* interpolator for future calls to {@link #storeTime storeTime}
|
||||
*/
|
||||
public void shift() {
|
||||
globalPreviousState = globalCurrentState;
|
||||
softPreviousState = globalPreviousState;
|
||||
softCurrentState = globalCurrentState;
|
||||
}
|
||||
|
||||
/** Store the current step state.
|
||||
* @param state current state
|
||||
*/
|
||||
public void storeState(final FieldODEStateAndDerivative<T> state) {
|
||||
|
||||
globalCurrentState = state;
|
||||
softCurrentState = globalCurrentState;
|
||||
h = globalCurrentState.getTime().subtract(globalPreviousState.getTime());
|
||||
|
||||
// the step is not finalized anymore
|
||||
finalized = false;
|
||||
|
||||
}
|
||||
|
||||
/** Restrict step range to a limited part of the global step.
|
||||
* <p>
|
||||
* This method can be used to restrict a step and make it appear
|
||||
* as if the original step was smaller. Calling this method
|
||||
* <em>only</em> changes the value returned by {@link #getPreviousState()},
|
||||
* it does not change any other property
|
||||
* </p>
|
||||
* @param softPreviousState start of the restricted step
|
||||
*/
|
||||
public void setSoftPreviousState(final FieldODEStateAndDerivative<T> softPreviousState) {
|
||||
this.softPreviousState = softPreviousState;
|
||||
}
|
||||
|
||||
/** Restrict step range to a limited part of the global step.
|
||||
* <p>
|
||||
* This method can be used to restrict a step and make it appear
|
||||
* as if the original step was smaller. Calling this method
|
||||
* <em>only</em> changes the value returned by {@link #getCurrentState()},
|
||||
* it does not change any other property
|
||||
* </p>
|
||||
* @param softCurrentState end of the restricted step
|
||||
*/
|
||||
public void setSoftCurrentState(final FieldODEStateAndDerivative<T> softCurrentState) {
|
||||
this.softCurrentState = softCurrentState;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the previous global grid point state.
|
||||
* @return previous global grid point state
|
||||
*/
|
||||
public FieldODEStateAndDerivative<T> getGlobalPreviousState() {
|
||||
return globalPreviousState;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the current global grid point state.
|
||||
* @return current global grid point state
|
||||
*/
|
||||
public FieldODEStateAndDerivative<T> getGlobalCurrentState() {
|
||||
return globalCurrentState;
|
||||
}
|
||||
|
||||
/** {@inheritDoc}
|
||||
* @see #setSoftPreviousState(FieldODEStateAndDerivative)
|
||||
*/
|
||||
public FieldODEStateAndDerivative<T> getPreviousState() {
|
||||
return softPreviousState;
|
||||
}
|
||||
|
||||
/** {@inheritDoc}
|
||||
* @see #setSoftCurrentState(FieldODEStateAndDerivative)
|
||||
*/
|
||||
public FieldODEStateAndDerivative<T> getCurrentState() {
|
||||
return softCurrentState;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public FieldODEStateAndDerivative<T> getInterpolatedState(final T time) {
|
||||
final T oneMinusThetaH = globalCurrentState.getTime().subtract(time);
|
||||
final T theta = (h.getReal() == 0) ? h.getField().getZero() : h.subtract(oneMinusThetaH).divide(h);
|
||||
return computeInterpolatedStateAndDerivatives(equations, theta, oneMinusThetaH);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public boolean isForward() {
|
||||
return forward;
|
||||
}
|
||||
|
||||
/** Compute the state and derivatives at the interpolated time.
|
||||
* This is the main processing method that should be implemented by
|
||||
* the derived classes to perform the interpolation.
|
||||
* @param eqn ODE equations with primary and secondary components
|
||||
* @param theta normalized interpolation abscissa within the step
|
||||
* (theta is zero at the previous time step and one at the current time step)
|
||||
* @param oneMinusThetaH time gap between the interpolated time and
|
||||
* the current time
|
||||
* @return interpolated state and derivatives
|
||||
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
|
||||
*/
|
||||
protected abstract FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(FieldExpandableODE<T> eqn,
|
||||
T theta,
|
||||
T oneMinusThetaH)
|
||||
throws MaxCountExceededException;
|
||||
|
||||
/**
|
||||
* Finalize the step.
|
||||
|
||||
* <p>Some embedded Runge-Kutta integrators need fewer functions
|
||||
* evaluations than their counterpart step interpolators. These
|
||||
* interpolators should perform the last evaluations they need by
|
||||
* themselves only if they need them. This method triggers these
|
||||
* extra evaluations. It can be called directly by the user step
|
||||
* handler and it is called automatically if {@link
|
||||
* #setInterpolatedTime} is called.</p>
|
||||
|
||||
* <p>Once this method has been called, <strong>no</strong> other
|
||||
* evaluation will be performed on this step. If there is a need to
|
||||
* have some side effects between the step handler and the
|
||||
* differential equations (for example update some data in the
|
||||
* equations once the step has been done), it is advised to call
|
||||
* this method explicitly from the step handler before these side
|
||||
* effects are set up. If the step handler induces no side effect,
|
||||
* then this method can safely be ignored, it will be called
|
||||
* transparently as needed.</p>
|
||||
|
||||
* <p><strong>Warning</strong>: since the step interpolator provided
|
||||
* to the step handler as a parameter of the {@link
|
||||
* StepHandler#handleStep handleStep} is valid only for the duration
|
||||
* of the {@link StepHandler#handleStep handleStep} call, one cannot
|
||||
* simply store a reference and reuse it later. One should first
|
||||
* finalize the instance, then copy this finalized instance into a
|
||||
* new object that can be kept.</p>
|
||||
|
||||
* <p>This method calls the protected <code>doFinalize</code> method
|
||||
* if it has never been called during this step and set a flag
|
||||
* indicating that it has been called once. It is the <code>
|
||||
* doFinalize</code> method which should perform the evaluations.
|
||||
* This wrapping prevents from calling <code>doFinalize</code> several
|
||||
* times and hence evaluating the differential equations too often.
|
||||
* Therefore, subclasses are not allowed not reimplement it, they
|
||||
* should rather reimplement <code>doFinalize</code>.</p>
|
||||
|
||||
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
|
||||
|
||||
*/
|
||||
public final void finalizeStep() throws MaxCountExceededException {
|
||||
if (! finalized) {
|
||||
doFinalize();
|
||||
finalized = true;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Really finalize the step.
|
||||
* The default implementation of this method does nothing.
|
||||
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
|
||||
*/
|
||||
protected void doFinalize() throws MaxCountExceededException {
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -17,8 +17,6 @@
|
|||
|
||||
package org.apache.commons.math4.ode.sampling;
|
||||
|
||||
import java.io.Externalizable;
|
||||
|
||||
import org.apache.commons.math4.RealFieldElement;
|
||||
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
|
||||
|
||||
|
|
@ -38,7 +36,7 @@ import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
|
|||
* @since 3.6
|
||||
*/
|
||||
|
||||
public interface FieldStepInterpolator<T extends RealFieldElement<T>> extends Externalizable {
|
||||
public interface FieldStepInterpolator<T extends RealFieldElement<T>> {
|
||||
|
||||
/**
|
||||
* Get the state at previous grid point time.
|
||||
|
|
|
|||
|
|
@ -92,6 +92,7 @@ import org.apache.commons.math4.util.Precision;
|
|||
*/
|
||||
|
||||
public class FieldStepNormalizer<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
|
||||
|
||||
/** Fixed time step. */
|
||||
private double h;
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue