diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/ode/sampling')
15 files changed, 2492 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java new file mode 100644 index 0000000..e674752 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java @@ -0,0 +1,171 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.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.math3.ode.FirstOrderFieldIntegrator + * @see StepHandler + * + * @param <T> the type of the field elements + * @since 3.6 + */ + +public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T>> + implements FieldStepInterpolator<T> { + + /** Global previous state. */ + private final FieldODEStateAndDerivative<T> globalPreviousState; + + /** Global current state. */ + private final FieldODEStateAndDerivative<T> globalCurrentState; + + /** Soft previous state. */ + private final FieldODEStateAndDerivative<T> softPreviousState; + + /** Soft current state. */ + private final FieldODEStateAndDerivative<T> softCurrentState; + + /** integration direction. */ + private final boolean forward; + + /** Mapper for ODE equations primary and secondary components. */ + private FieldEquationsMapper<T> mapper; + + /** Simple constructor. + * @param isForward integration direction indicator + * @param globalPreviousState start of the global step + * @param globalCurrentState end of the global step + * @param softPreviousState start of the restricted step + * @param softCurrentState end of the restricted step + * @param equationsMapper mapper for ODE equations primary and secondary components + */ + protected AbstractFieldStepInterpolator(final boolean isForward, + final FieldODEStateAndDerivative<T> globalPreviousState, + final FieldODEStateAndDerivative<T> globalCurrentState, + final FieldODEStateAndDerivative<T> softPreviousState, + final FieldODEStateAndDerivative<T> softCurrentState, + final FieldEquationsMapper<T> equationsMapper) { + this.forward = isForward; + this.globalPreviousState = globalPreviousState; + this.globalCurrentState = globalCurrentState; + this.softPreviousState = softPreviousState; + this.softCurrentState = softCurrentState; + this.mapper = equationsMapper; + } + + /** Create a new restricted version of the instance. + * <p> + * The instance is not changed at all. + * </p> + * @param previousState start of the restricted step + * @param currentState end of the restricted step + * @return restricted version of the instance + * @see #getPreviousState() + * @see #getCurrentState() + */ + public AbstractFieldStepInterpolator<T> restrictStep(final FieldODEStateAndDerivative<T> previousState, + final FieldODEStateAndDerivative<T> currentState) { + return create(forward, globalPreviousState, globalCurrentState, previousState, currentState, mapper); + } + + /** Create a new instance. + * @param newForward integration direction indicator + * @param newGlobalPreviousState start of the global step + * @param newGlobalCurrentState end of the global step + * @param newSoftPreviousState start of the restricted step + * @param newSoftCurrentState end of the restricted step + * @param newMapper equations mapper for the all equations + * @return a new instance + */ + protected abstract AbstractFieldStepInterpolator<T> create(boolean newForward, + FieldODEStateAndDerivative<T> newGlobalPreviousState, + FieldODEStateAndDerivative<T> newGlobalCurrentState, + FieldODEStateAndDerivative<T> newSoftPreviousState, + FieldODEStateAndDerivative<T> newSoftCurrentState, + FieldEquationsMapper<T> newMapper); + + /** + * 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} */ + public FieldODEStateAndDerivative<T> getPreviousState() { + return softPreviousState; + } + + /** {@inheritDoc} */ + public FieldODEStateAndDerivative<T> getCurrentState() { + return softCurrentState; + } + + /** {@inheritDoc} */ + public FieldODEStateAndDerivative<T> getInterpolatedState(final T time) { + final T thetaH = time.subtract(globalPreviousState.getTime()); + final T oneMinusThetaH = globalCurrentState.getTime().subtract(time); + final T theta = thetaH.divide(globalCurrentState.getTime().subtract(globalPreviousState.getTime())); + return computeInterpolatedStateAndDerivatives(mapper, time, theta, thetaH, 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 equationsMapper mapper for ODE equations primary and secondary components + * @param time interpolation time + * @param theta normalized interpolation abscissa within the step + * (theta is zero at the previous time step and one at the current time step) + * @param thetaH time gap between the previous time and the interpolated time + * @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(FieldEquationsMapper<T> equationsMapper, + T time, T theta, + T thetaH, T oneMinusThetaH) + throws MaxCountExceededException; + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/AbstractStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractStepInterpolator.java new file mode 100644 index 0000000..1fbc04a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractStepInterpolator.java @@ -0,0 +1,605 @@ +/* + * 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.sampling; + +import java.io.IOException; +import java.io.ObjectInput; +import java.io.ObjectOutput; + +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.ode.EquationsMapper; + +/** 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.math3.ode.FirstOrderIntegrator + * @see org.apache.commons.math3.ode.SecondOrderIntegrator + * @see StepHandler + * + * @since 1.2 + * + */ + +public abstract class AbstractStepInterpolator + implements StepInterpolator { + + /** current time step */ + protected double h; + + /** current state */ + protected double[] currentState; + + /** interpolated time */ + protected double interpolatedTime; + + /** interpolated state */ + protected double[] interpolatedState; + + /** interpolated derivatives */ + protected double[] interpolatedDerivatives; + + /** interpolated primary state */ + protected double[] interpolatedPrimaryState; + + /** interpolated primary derivatives */ + protected double[] interpolatedPrimaryDerivatives; + + /** interpolated secondary state */ + protected double[][] interpolatedSecondaryState; + + /** interpolated secondary derivatives */ + protected double[][] interpolatedSecondaryDerivatives; + + /** global previous time */ + private double globalPreviousTime; + + /** global current time */ + private double globalCurrentTime; + + /** soft previous time */ + private double softPreviousTime; + + /** soft current time */ + private double softCurrentTime; + + /** indicate if the step has been finalized or not. */ + private boolean finalized; + + /** integration direction. */ + private boolean forward; + + /** indicator for dirty state. */ + private boolean dirtyState; + + /** Equations mapper for the primary equations set. */ + private EquationsMapper primaryMapper; + + /** Equations mappers for the secondary equations sets. */ + private EquationsMapper[] secondaryMappers; + + /** 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.math3.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 AbstractStepInterpolator() { + globalPreviousTime = Double.NaN; + globalCurrentTime = Double.NaN; + softPreviousTime = Double.NaN; + softCurrentTime = Double.NaN; + h = Double.NaN; + interpolatedTime = Double.NaN; + currentState = null; + finalized = false; + this.forward = true; + this.dirtyState = true; + primaryMapper = null; + secondaryMappers = null; + allocateInterpolatedArrays(-1); + } + + /** Simple constructor. + * @param y reference to the integrator array holding the state at + * the end of the step + * @param forward integration direction indicator + * @param primaryMapper equations mapper for the primary equations set + * @param secondaryMappers equations mappers for the secondary equations sets + */ + protected AbstractStepInterpolator(final double[] y, final boolean forward, + final EquationsMapper primaryMapper, + final EquationsMapper[] secondaryMappers) { + + globalPreviousTime = Double.NaN; + globalCurrentTime = Double.NaN; + softPreviousTime = Double.NaN; + softCurrentTime = Double.NaN; + h = Double.NaN; + interpolatedTime = Double.NaN; + currentState = y; + finalized = false; + this.forward = forward; + this.dirtyState = true; + this.primaryMapper = primaryMapper; + this.secondaryMappers = (secondaryMappers == null) ? null : secondaryMappers.clone(); + allocateInterpolatedArrays(y.length); + + } + + /** 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 AbstractStepInterpolator(final AbstractStepInterpolator interpolator) { + + globalPreviousTime = interpolator.globalPreviousTime; + globalCurrentTime = interpolator.globalCurrentTime; + softPreviousTime = interpolator.softPreviousTime; + softCurrentTime = interpolator.softCurrentTime; + h = interpolator.h; + interpolatedTime = interpolator.interpolatedTime; + + if (interpolator.currentState == null) { + currentState = null; + primaryMapper = null; + secondaryMappers = null; + allocateInterpolatedArrays(-1); + } else { + currentState = interpolator.currentState.clone(); + interpolatedState = interpolator.interpolatedState.clone(); + interpolatedDerivatives = interpolator.interpolatedDerivatives.clone(); + interpolatedPrimaryState = interpolator.interpolatedPrimaryState.clone(); + interpolatedPrimaryDerivatives = interpolator.interpolatedPrimaryDerivatives.clone(); + interpolatedSecondaryState = new double[interpolator.interpolatedSecondaryState.length][]; + interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][]; + for (int i = 0; i < interpolatedSecondaryState.length; ++i) { + interpolatedSecondaryState[i] = interpolator.interpolatedSecondaryState[i].clone(); + interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone(); + } + } + + finalized = interpolator.finalized; + forward = interpolator.forward; + dirtyState = interpolator.dirtyState; + primaryMapper = interpolator.primaryMapper; + secondaryMappers = (interpolator.secondaryMappers == null) ? + null : interpolator.secondaryMappers.clone(); + + } + + /** Allocate the various interpolated states arrays. + * @param dimension total dimension (negative if arrays should be set to null) + */ + private void allocateInterpolatedArrays(final int dimension) { + if (dimension < 0) { + interpolatedState = null; + interpolatedDerivatives = null; + interpolatedPrimaryState = null; + interpolatedPrimaryDerivatives = null; + interpolatedSecondaryState = null; + interpolatedSecondaryDerivatives = null; + } else { + interpolatedState = new double[dimension]; + interpolatedDerivatives = new double[dimension]; + interpolatedPrimaryState = new double[primaryMapper.getDimension()]; + interpolatedPrimaryDerivatives = new double[primaryMapper.getDimension()]; + if (secondaryMappers == null) { + interpolatedSecondaryState = null; + interpolatedSecondaryDerivatives = null; + } else { + interpolatedSecondaryState = new double[secondaryMappers.length][]; + interpolatedSecondaryDerivatives = new double[secondaryMappers.length][]; + for (int i = 0; i < secondaryMappers.length; ++i) { + interpolatedSecondaryState[i] = new double[secondaryMappers[i].getDimension()]; + interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()]; + } + } + } + } + + /** 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 primary equations mapper for the primary equations set + * @param secondary equations mappers for the secondary equations sets + */ + protected void reinitialize(final double[] y, final boolean isForward, + final EquationsMapper primary, + final EquationsMapper[] secondary) { + + globalPreviousTime = Double.NaN; + globalCurrentTime = Double.NaN; + softPreviousTime = Double.NaN; + softCurrentTime = Double.NaN; + h = Double.NaN; + interpolatedTime = Double.NaN; + currentState = y; + finalized = false; + this.forward = isForward; + this.dirtyState = true; + this.primaryMapper = primary; + this.secondaryMappers = secondary.clone(); + allocateInterpolatedArrays(y.length); + + } + + /** {@inheritDoc} */ + public StepInterpolator 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 StepInterpolator 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() { + globalPreviousTime = globalCurrentTime; + softPreviousTime = globalPreviousTime; + softCurrentTime = globalCurrentTime; + } + + /** Store the current step time. + * @param t current time + */ + public void storeTime(final double t) { + + globalCurrentTime = t; + softCurrentTime = globalCurrentTime; + h = globalCurrentTime - globalPreviousTime; + setInterpolatedTime(t); + + // 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 #getPreviousTime()}, + * it does not change any other property + * </p> + * @param softPreviousTime start of the restricted step + * @since 2.2 + */ + public void setSoftPreviousTime(final double softPreviousTime) { + this.softPreviousTime = softPreviousTime; + } + + /** 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 #getCurrentTime()}, + * it does not change any other property + * </p> + * @param softCurrentTime end of the restricted step + * @since 2.2 + */ + public void setSoftCurrentTime(final double softCurrentTime) { + this.softCurrentTime = softCurrentTime; + } + + /** + * Get the previous global grid point time. + * @return previous global grid point time + */ + public double getGlobalPreviousTime() { + return globalPreviousTime; + } + + /** + * Get the current global grid point time. + * @return current global grid point time + */ + public double getGlobalCurrentTime() { + return globalCurrentTime; + } + + /** + * Get the previous soft grid point time. + * @return previous soft grid point time + * @see #setSoftPreviousTime(double) + */ + public double getPreviousTime() { + return softPreviousTime; + } + + /** + * Get the current soft grid point time. + * @return current soft grid point time + * @see #setSoftCurrentTime(double) + */ + public double getCurrentTime() { + return softCurrentTime; + } + + /** {@inheritDoc} */ + public double getInterpolatedTime() { + return interpolatedTime; + } + + /** {@inheritDoc} */ + public void setInterpolatedTime(final double time) { + interpolatedTime = time; + dirtyState = true; + } + + /** {@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 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 + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + protected abstract void computeInterpolatedStateAndDerivatives(double theta, + double oneMinusThetaH) + throws MaxCountExceededException; + + /** Lazy evaluation of complete interpolated state. + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + private void evaluateCompleteInterpolatedState() + throws MaxCountExceededException { + // lazy evaluation of the state + if (dirtyState) { + final double oneMinusThetaH = globalCurrentTime - interpolatedTime; + final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h; + computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH); + dirtyState = false; + } + } + + /** {@inheritDoc} */ + public double[] getInterpolatedState() throws MaxCountExceededException { + evaluateCompleteInterpolatedState(); + primaryMapper.extractEquationData(interpolatedState, + interpolatedPrimaryState); + return interpolatedPrimaryState; + } + + /** {@inheritDoc} */ + public double[] getInterpolatedDerivatives() throws MaxCountExceededException { + evaluateCompleteInterpolatedState(); + primaryMapper.extractEquationData(interpolatedDerivatives, + interpolatedPrimaryDerivatives); + return interpolatedPrimaryDerivatives; + } + + /** {@inheritDoc} */ + public double[] getInterpolatedSecondaryState(final int index) throws MaxCountExceededException { + evaluateCompleteInterpolatedState(); + secondaryMappers[index].extractEquationData(interpolatedState, + interpolatedSecondaryState[index]); + return interpolatedSecondaryState[index]; + } + + /** {@inheritDoc} */ + public double[] getInterpolatedSecondaryDerivatives(final int index) throws MaxCountExceededException { + evaluateCompleteInterpolatedState(); + secondaryMappers[index].extractEquationData(interpolatedDerivatives, + interpolatedSecondaryDerivatives[index]); + return interpolatedSecondaryDerivatives[index]; + } + + /** + * 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 { + } + + /** {@inheritDoc} */ + public abstract void writeExternal(ObjectOutput out) + throws IOException; + + /** {@inheritDoc} */ + public abstract void readExternal(ObjectInput in) + throws IOException, ClassNotFoundException; + + /** Save the base state of the instance. + * This method performs step finalization if it has not been done + * before. + * @param out stream where to save the state + * @exception IOException in case of write error + */ + protected void writeBaseExternal(final ObjectOutput out) + throws IOException { + + if (currentState == null) { + out.writeInt(-1); + } else { + out.writeInt(currentState.length); + } + out.writeDouble(globalPreviousTime); + out.writeDouble(globalCurrentTime); + out.writeDouble(softPreviousTime); + out.writeDouble(softCurrentTime); + out.writeDouble(h); + out.writeBoolean(forward); + out.writeObject(primaryMapper); + out.write(secondaryMappers.length); + for (final EquationsMapper mapper : secondaryMappers) { + out.writeObject(mapper); + } + + if (currentState != null) { + for (int i = 0; i < currentState.length; ++i) { + out.writeDouble(currentState[i]); + } + } + + out.writeDouble(interpolatedTime); + + // we do not store the interpolated state, + // it will be recomputed as needed after reading + + try { + // finalize the step (and don't bother saving the now true flag) + finalizeStep(); + } catch (MaxCountExceededException mcee) { + final IOException ioe = new IOException(mcee.getLocalizedMessage()); + ioe.initCause(mcee); + throw ioe; + } + + } + + /** Read the base state of the instance. + * This method does <strong>neither</strong> set the interpolated + * time nor state. It is up to the derived class to reset it + * properly calling the {@link #setInterpolatedTime} method later, + * once all rest of the object state has been set up properly. + * @param in stream where to read the state from + * @return interpolated time to be set later by the caller + * @exception IOException in case of read error + * @exception ClassNotFoundException if an equation mapper class + * cannot be found + */ + protected double readBaseExternal(final ObjectInput in) + throws IOException, ClassNotFoundException { + + final int dimension = in.readInt(); + globalPreviousTime = in.readDouble(); + globalCurrentTime = in.readDouble(); + softPreviousTime = in.readDouble(); + softCurrentTime = in.readDouble(); + h = in.readDouble(); + forward = in.readBoolean(); + primaryMapper = (EquationsMapper) in.readObject(); + secondaryMappers = new EquationsMapper[in.read()]; + for (int i = 0; i < secondaryMappers.length; ++i) { + secondaryMappers[i] = (EquationsMapper) in.readObject(); + } + dirtyState = true; + + if (dimension < 0) { + currentState = null; + } else { + currentState = new double[dimension]; + for (int i = 0; i < currentState.length; ++i) { + currentState[i] = in.readDouble(); + } + } + + // we do NOT handle the interpolated time and state here + interpolatedTime = Double.NaN; + allocateInterpolatedArrays(dimension); + + finalized = true; + + return in.readDouble(); + + } + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/DummyStepHandler.java b/src/main/java/org/apache/commons/math3/ode/sampling/DummyStepHandler.java new file mode 100644 index 0000000..c1d4fe3 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/DummyStepHandler.java @@ -0,0 +1,89 @@ +/* + * 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.sampling; + +/** + * This class is a step handler that does nothing. + + * <p>This class is provided as a convenience for users who are only + * interested in the final state of an integration and not in the + * intermediate steps. Its handleStep method does nothing.</p> + * + * <p>Since this class has no internal state, it is implemented using + * the Singleton design pattern. This means that only one instance is + * ever created, which can be retrieved using the getInstance + * method. This explains why there is no public constructor.</p> + * + * @see StepHandler + * @since 1.2 + */ + +public class DummyStepHandler implements StepHandler { + + /** Private constructor. + * The constructor is private to prevent users from creating + * instances (Singleton design-pattern). + */ + private DummyStepHandler() { + } + + /** Get the only instance. + * @return the only instance + */ + public static DummyStepHandler getInstance() { + return LazyHolder.INSTANCE; + } + + /** {@inheritDoc} */ + public void init(double t0, double[] y0, double t) { + } + + /** + * Handle the last accepted step. + * This method does nothing in this class. + * @param interpolator interpolator for the last accepted step. For + * efficiency purposes, the various integrators reuse the same + * object on each call, so if the instance wants to keep it across + * all calls (for example to provide at the end of the integration a + * continuous model valid throughout the integration range), it + * should build a local copy using the clone method and store this + * copy. + * @param isLast true if the step is the last one + */ + public void handleStep(final StepInterpolator interpolator, final boolean isLast) { + } + + // CHECKSTYLE: stop HideUtilityClassConstructor + /** Holder for the instance. + * <p>We use here the Initialization On Demand Holder Idiom.</p> + */ + private static class LazyHolder { + /** Cached field instance. */ + private static final DummyStepHandler INSTANCE = new DummyStepHandler(); + } + // CHECKSTYLE: resume HideUtilityClassConstructor + + /** Handle deserialization of the singleton. + * @return the singleton instance + */ + private Object readResolve() { + // return the singleton instance + return LazyHolder.INSTANCE; + } + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/FieldFixedStepHandler.java b/src/main/java/org/apache/commons/math3/ode/sampling/FieldFixedStepHandler.java new file mode 100644 index 0000000..9d4fd70 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/FieldFixedStepHandler.java @@ -0,0 +1,69 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; + +/** + * This interface represents a handler that should be called after + * each successful fixed step. + + * <p>This interface should be implemented by anyone who is interested + * in getting the solution of an ordinary differential equation at + * fixed time steps. Objects implementing this interface should be + * wrapped within an instance of {@link FieldStepNormalizer} that itself + * is used as the general {@link FieldStepHandler} by the integrator. The + * {@link FieldStepNormalizer} object is called according to the integrator + * internal algorithms and it calls objects implementing this + * interface as necessary at fixed time steps.</p> + * + * @see FieldStepHandler + * @see FieldStepNormalizer + * @see FieldStepInterpolator + * @param <T> the type of the field elements + * @since 3.6 + */ + +public interface FieldFixedStepHandler<T extends RealFieldElement<T>> { + + /** Initialize step handler at the start of an ODE integration. + * <p> + * This method is called once at the start of the integration. It + * may be used by the step handler to initialize some internal data + * if needed. + * </p> + * @param initialState initial time, state vector and derivative + * @param finalTime target time for the integration + */ + void init(FieldODEStateAndDerivative<T> initialState, T finalTime); + + /** + * Handle the last accepted step + * @param state current value of the independent <i>time</i> variable, + * state vector and derivative + * For efficiency purposes, the {@link FieldStepNormalizer} class reuses + * the same array on each call, so if + * the instance wants to keep it across all calls (for example to + * provide at the end of the integration a complete array of all + * steps), it should build a local copy store this copy. + * @param isLast true if the step is the last one + */ + void handleStep(FieldODEStateAndDerivative<T> state, boolean isLast); + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepHandler.java b/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepHandler.java new file mode 100644 index 0000000..c56911b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepHandler.java @@ -0,0 +1,75 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; + +/** + * This interface represents a handler that should be called after + * each successful step. + * + * <p>The ODE integrators compute the evolution of the state vector at + * some grid points that depend on their own internal algorithm. Once + * they have found a new grid point (possibly after having computed + * several evaluation of the derivative at intermediate points), they + * provide it to objects implementing this interface. These objects + * typically either ignore the intermediate steps and wait for the + * last one, store the points in an ephemeris, or forward them to + * specialized processing or output methods.</p> + * + * @see org.apache.commons.math3.ode.FirstOrderFieldIntegrator + * @see FieldStepInterpolator + * @param <T> the type of the field elements + * @since 3.6 + */ + +public interface FieldStepHandler<T extends RealFieldElement<T>> { + + /** Initialize step handler at the start of an ODE integration. + * <p> + * This method is called once at the start of the integration. It + * may be used by the step handler to initialize some internal data + * if needed. + * </p> + * @param initialState initial time, state vector and derivative + * @param finalTime target time for the integration + */ + void init(FieldODEStateAndDerivative<T> initialState, T finalTime); + + /** + * Handle the last accepted step + * @param interpolator interpolator for the last accepted step. For + * efficiency purposes, the various integrators reuse the same + * object on each call, so if the instance wants to keep it across + * all calls (for example to provide at the end of the integration a + * continuous model valid throughout the integration range, as the + * {@link org.apache.commons.math3.ode.ContinuousOutputModel + * ContinuousOutputModel} class does), it should build a local copy + * using the clone method of the interpolator and store this copy. + * Keeping only a reference to the interpolator and reusing it will + * result in unpredictable behavior (potentially crashing the application). + * @param isLast true if the step is the last one + * @exception MaxCountExceededException if the interpolator throws one because + * the number of functions evaluations is exceeded + */ + void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) + throws MaxCountExceededException; + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepInterpolator.java new file mode 100644 index 0000000..a005fb1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepInterpolator.java @@ -0,0 +1,76 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; + +/** This interface represents an interpolator over the last step + * during an ODE integration. + * + * <p>The various ODE integrators provide objects implementing this + * interface to the step handlers. These objects are often custom + * objects tightly bound to the integrator internal algorithms. The + * handlers can use these objects to retrieve the state vector at + * intermediate times between the previous and the current grid points + * (this feature is often called dense output).</p> + * + * @param <T> the type of the field elements + * @see org.apache.commons.math3.ode.FirstOrderFieldIntegrator + * @see FieldStepHandler + * @since 3.6 + */ + +public interface FieldStepInterpolator<T extends RealFieldElement<T>> { + + /** + * Get the state at previous grid point time. + * @return state at previous grid point time + */ + FieldODEStateAndDerivative<T> getPreviousState(); + + /** + * Get the state at current grid point time. + * @return state at current grid point time + */ + FieldODEStateAndDerivative<T> getCurrentState(); + + /** + * Get the state at interpolated time. + * <p>Setting the time outside of the current step is allowed, but + * should be used with care since the accuracy of the interpolator will + * probably be very poor far from this step. This allowance has been + * added to simplify implementation of search algorithms near the + * step endpoints.</p> + * @param time time of the interpolated point + * @return state at interpolated time + */ + FieldODEStateAndDerivative<T> getInterpolatedState(T time); + + /** Check if the natural integration direction is forward. + * <p>This method provides the integration direction as specified by + * the integrator itself, it avoid some nasty problems in + * degenerated cases like null steps due to cancellation at step + * initialization, step control or discrete events + * triggering.</p> + * @return true if the integration variable (time) increases during + * integration + */ + boolean isForward(); + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepNormalizer.java b/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepNormalizer.java new file mode 100644 index 0000000..892fcf7 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/FieldStepNormalizer.java @@ -0,0 +1,273 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; + +/** + * This class wraps an object implementing {@link FieldFixedStepHandler} + * into a {@link FieldStepHandler}. + + * <p>This wrapper allows to use fixed step handlers with general + * integrators which cannot guaranty their integration steps will + * remain constant and therefore only accept general step + * handlers.</p> + * + * <p>The stepsize used is selected at construction time. The {@link + * FieldFixedStepHandler#handleStep handleStep} method of the underlying + * {@link FieldFixedStepHandler} object is called at normalized times. The + * normalized times can be influenced by the {@link StepNormalizerMode} and + * {@link StepNormalizerBounds}.</p> + * + * <p>There is no constraint on the integrator, it can use any time step + * it needs (time steps longer or shorter than the fixed time step and + * non-integer ratios are all allowed).</p> + * + * <p> + * <table border="1" align="center"> + * <tr BGCOLOR="#CCCCFF"><td colspan=6><font size="+2">Examples (step size = 0.5)</font></td></tr> + * <tr BGCOLOR="#EEEEFF"><font size="+1"><td>Start time</td><td>End time</td> + * <td>Direction</td><td>{@link StepNormalizerMode Mode}</td> + * <td>{@link StepNormalizerBounds Bounds}</td><td>Output</td></font></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.8, 1.3, 1.8, 2.3, 2.8</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.3, 0.8, 1.3, 1.8, 2.3, 2.8</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.8, 1.3, 1.8, 2.3, 2.8, 3.1</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.3, 0.8, 1.3, 1.8, 2.3, 2.8, 3.1</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.1</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.1</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.6, 2.1, 1.6, 1.1, 0.6</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.1, 2.6, 2.1, 1.6, 1.1, 0.6</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.6, 2.1, 1.6, 1.1, 0.6, 0.3</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.1, 2.6, 2.1, 1.6, 1.1, 0.6, 0.3</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.1, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.3</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.1, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.3</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * </table> + * </p> + * + * @param <T> the type of the field elements + * @see FieldStepHandler + * @see FieldFixedStepHandler + * @see StepNormalizerMode + * @see StepNormalizerBounds + * @since 3.6 + */ + +public class FieldStepNormalizer<T extends RealFieldElement<T>> implements FieldStepHandler<T> { + + /** Fixed time step. */ + private double h; + + /** Underlying step handler. */ + private final FieldFixedStepHandler<T> handler; + + /** First step state. */ + private FieldODEStateAndDerivative<T> first; + + /** Last step step. */ + private FieldODEStateAndDerivative<T> last; + + /** Integration direction indicator. */ + private boolean forward; + + /** The step normalizer bounds settings to use. */ + private final StepNormalizerBounds bounds; + + /** The step normalizer mode to use. */ + private final StepNormalizerMode mode; + + /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT} + * mode, and {@link StepNormalizerBounds#FIRST FIRST} bounds setting, for + * backwards compatibility. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + */ + public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler) { + this(h, handler, StepNormalizerMode.INCREMENT, + StepNormalizerBounds.FIRST); + } + + /** Simple constructor. Uses {@link StepNormalizerBounds#FIRST FIRST} + * bounds setting. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + * @param mode step normalizer mode to use + * @since 3.0 + */ + public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler, + final StepNormalizerMode mode) { + this(h, handler, mode, StepNormalizerBounds.FIRST); + } + + /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT} + * mode. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + * @param bounds step normalizer bounds setting to use + * @since 3.0 + */ + public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler, + final StepNormalizerBounds bounds) { + this(h, handler, StepNormalizerMode.INCREMENT, bounds); + } + + /** Simple constructor. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + * @param mode step normalizer mode to use + * @param bounds step normalizer bounds setting to use + * @since 3.0 + */ + public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler, + final StepNormalizerMode mode, final StepNormalizerBounds bounds) { + this.h = FastMath.abs(h); + this.handler = handler; + this.mode = mode; + this.bounds = bounds; + first = null; + last = null; + forward = true; + } + + /** {@inheritDoc} */ + public void init(final FieldODEStateAndDerivative<T> initialState, final T finalTime) { + + first = null; + last = null; + forward = true; + + // initialize the underlying handler + handler.init(initialState, finalTime); + + } + + /** + * Handle the last accepted step + * @param interpolator interpolator for the last accepted step. For + * efficiency purposes, the various integrators reuse the same + * object on each call, so if the instance wants to keep it across + * all calls (for example to provide at the end of the integration a + * continuous model valid throughout the integration range), it + * should build a local copy using the clone method and store this + * copy. + * @param isLast true if the step is the last one + * @exception MaxCountExceededException if the interpolator throws one because + * the number of functions evaluations is exceeded + */ + public void handleStep(final FieldStepInterpolator<T> interpolator, final boolean isLast) + throws MaxCountExceededException { + // The first time, update the last state with the start information. + if (last == null) { + + first = interpolator.getPreviousState(); + last = first; + + // Take the integration direction into account. + forward = interpolator.isForward(); + if (!forward) { + h = -h; + } + } + + // Calculate next normalized step time. + T nextTime = (mode == StepNormalizerMode.INCREMENT) ? + last.getTime().add(h) : + last.getTime().getField().getZero().add((FastMath.floor(last.getTime().getReal() / h) + 1) * h); + if (mode == StepNormalizerMode.MULTIPLES && + Precision.equals(nextTime.getReal(), last.getTime().getReal(), 1)) { + nextTime = nextTime.add(h); + } + + // Process normalized steps as long as they are in the current step. + boolean nextInStep = isNextInStep(nextTime, interpolator); + while (nextInStep) { + // Output the stored previous step. + doNormalizedStep(false); + + // Store the next step as last step. + last = interpolator.getInterpolatedState(nextTime); + + // Move on to the next step. + nextTime = nextTime.add(h); + nextInStep = isNextInStep(nextTime, interpolator); + } + + if (isLast) { + // There will be no more steps. The stored one should be given to + // the handler. We may have to output one more step. Only the last + // one of those should be flagged as being the last. + final boolean addLast = bounds.lastIncluded() && + last.getTime().getReal() != interpolator.getCurrentState().getTime().getReal(); + doNormalizedStep(!addLast); + if (addLast) { + last = interpolator.getCurrentState(); + doNormalizedStep(true); + } + } + } + + /** + * Returns a value indicating whether the next normalized time is in the + * current step. + * @param nextTime the next normalized time + * @param interpolator interpolator for the last accepted step, to use to + * get the end time of the current step + * @return value indicating whether the next normalized time is in the + * current step + */ + private boolean isNextInStep(final T nextTime, final FieldStepInterpolator<T> interpolator) { + return forward ? + nextTime.getReal() <= interpolator.getCurrentState().getTime().getReal() : + nextTime.getReal() >= interpolator.getCurrentState().getTime().getReal(); + } + + /** + * Invokes the underlying step handler for the current normalized step. + * @param isLast true if the step is the last one + */ + private void doNormalizedStep(final boolean isLast) { + if (!bounds.firstIncluded() && first.getTime().getReal() == last.getTime().getReal()) { + return; + } + handler.handleStep(last, isLast); + } + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/FixedStepHandler.java b/src/main/java/org/apache/commons/math3/ode/sampling/FixedStepHandler.java new file mode 100644 index 0000000..bbe6ac5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/FixedStepHandler.java @@ -0,0 +1,71 @@ +/* + * 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.sampling; + + +/** + * This interface represents a handler that should be called after + * each successful fixed step. + + * <p>This interface should be implemented by anyone who is interested + * in getting the solution of an ordinary differential equation at + * fixed time steps. Objects implementing this interface should be + * wrapped within an instance of {@link StepNormalizer} that itself + * is used as the general {@link StepHandler} by the integrator. The + * {@link StepNormalizer} object is called according to the integrator + * internal algorithms and it calls objects implementing this + * interface as necessary at fixed time steps.</p> + * + * @see StepHandler + * @see StepNormalizer + * @since 1.2 + */ + +public interface FixedStepHandler { + + /** Initialize step handler at the start of an ODE integration. + * <p> + * This method is called once at the start of the integration. It + * may be used by the step handler to initialize some internal data + * if needed. + * </p> + * @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 + */ + void init(double t0, double[] y0, double t); + + /** + * Handle the last accepted step + * @param t time of the current step + * @param y state vector at t. For efficiency purposes, the {@link + * StepNormalizer} class reuses the same array on each call, so if + * the instance wants to keep it across all calls (for example to + * provide at the end of the integration a complete array of all + * steps), it should build a local copy store this copy. + * @param yDot derivatives of the state vector state vector at t. + * For efficiency purposes, the {@link StepNormalizer} class reuses + * the same array on each call, so if + * the instance wants to keep it across all calls (for example to + * provide at the end of the integration a complete array of all + * steps), it should build a local copy store this copy. + * @param isLast true if the step is the last one + */ + void handleStep(double t, double[] y, double[] yDot, boolean isLast); + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolator.java new file mode 100644 index 0000000..39b05ab --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolator.java @@ -0,0 +1,293 @@ +/* + * 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.sampling; + +import java.io.IOException; +import java.io.ObjectInput; +import java.io.ObjectOutput; +import java.util.Arrays; + +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.ode.EquationsMapper; +import org.apache.commons.math3.util.FastMath; + +/** + * This class implements an interpolator for integrators using Nordsieck representation. + * + * <p>This interpolator computes dense output around the current point. + * The interpolation equation is based on Taylor series formulas. + * + * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator + * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator + * @since 2.0 + */ + +public class NordsieckStepInterpolator extends AbstractStepInterpolator { + + /** Serializable version identifier */ + private static final long serialVersionUID = -7179861704951334960L; + + /** State variation. */ + protected double[] stateVariation; + + /** Step size used in the first scaled derivative and Nordsieck vector. */ + private double scalingH; + + /** Reference time for all arrays. + * <p>Sometimes, the reference time is the same as previousTime, + * sometimes it is the same as currentTime, so we use a separate + * field to avoid any confusion. + * </p> + */ + private double referenceTime; + + /** First scaled derivative. */ + private double[] scaled; + + /** Nordsieck vector. */ + private Array2DRowRealMatrix nordsieck; + + /** Simple constructor. + * This constructor builds an instance that is not usable yet, the + * {@link AbstractStepInterpolator#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. + */ + public NordsieckStepInterpolator() { + } + + /** Copy constructor. + * @param interpolator interpolator to copy from. The copy is a deep + * copy: its arrays are separated from the original arrays of the + * instance + */ + public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) { + super(interpolator); + scalingH = interpolator.scalingH; + referenceTime = interpolator.referenceTime; + if (interpolator.scaled != null) { + scaled = interpolator.scaled.clone(); + } + if (interpolator.nordsieck != null) { + nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true); + } + if (interpolator.stateVariation != null) { + stateVariation = interpolator.stateVariation.clone(); + } + } + + /** {@inheritDoc} */ + @Override + protected StepInterpolator doCopy() { + return new NordsieckStepInterpolator(this); + } + + /** Reinitialize the instance. + * <p>Beware that all arrays <em>must</em> be references to integrator + * arrays, in order to ensure proper update without copy.</p> + * @param y reference to the integrator array holding the state at + * the end of the step + * @param forward integration direction indicator + * @param primaryMapper equations mapper for the primary equations set + * @param secondaryMappers equations mappers for the secondary equations sets + */ + @Override + public void reinitialize(final double[] y, final boolean forward, + final EquationsMapper primaryMapper, + final EquationsMapper[] secondaryMappers) { + super.reinitialize(y, forward, primaryMapper, secondaryMappers); + stateVariation = new double[y.length]; + } + + /** Reinitialize the instance. + * <p>Beware that all arrays <em>must</em> be references to integrator + * arrays, in order to ensure proper update without copy.</p> + * @param time time at which all arrays are defined + * @param stepSize step size used in the scaled and Nordsieck arrays + * @param scaledDerivative reference to the integrator array holding the first + * scaled derivative + * @param nordsieckVector reference to the integrator matrix holding the + * Nordsieck vector + */ + public void reinitialize(final double time, final double stepSize, + final double[] scaledDerivative, + final Array2DRowRealMatrix nordsieckVector) { + this.referenceTime = time; + this.scalingH = stepSize; + this.scaled = scaledDerivative; + this.nordsieck = nordsieckVector; + + // make sure the state and derivatives will depend on the new arrays + setInterpolatedTime(getInterpolatedTime()); + + } + + /** Rescale the instance. + * <p>Since the scaled and Nordsieck arrays are shared with the caller, + * this method has the side effect of rescaling this arrays in the caller too.</p> + * @param stepSize new step size to use in the scaled and Nordsieck arrays + */ + public void rescale(final double stepSize) { + + final double ratio = stepSize / scalingH; + for (int i = 0; i < scaled.length; ++i) { + scaled[i] *= ratio; + } + + final double[][] nData = nordsieck.getDataRef(); + double power = ratio; + for (int i = 0; i < nData.length; ++i) { + power *= ratio; + final double[] nDataI = nData[i]; + for (int j = 0; j < nDataI.length; ++j) { + nDataI[j] *= power; + } + } + + scalingH = stepSize; + + } + + /** + * Get the state vector variation from current to interpolated state. + * <p>This method is aimed at computing y(t<sub>interpolation</sub>) + * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors + * that would occur if the subtraction were performed explicitly.</p> + * <p>The returned vector is a reference to a reused array, so + * it should not be modified and it should be copied if it needs + * to be preserved across several calls.</p> + * @return state vector at time {@link #getInterpolatedTime} + * @see #getInterpolatedDerivatives() + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + public double[] getInterpolatedStateVariation() throws MaxCountExceededException { + // compute and ignore interpolated state + // to make sure state variation is computed as a side effect + getInterpolatedState(); + return stateVariation; + } + + /** {@inheritDoc} */ + @Override + protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) { + + final double x = interpolatedTime - referenceTime; + final double normalizedAbscissa = x / scalingH; + + Arrays.fill(stateVariation, 0.0); + Arrays.fill(interpolatedDerivatives, 0.0); + + // apply Taylor formula from high order to low order, + // for the sake of numerical accuracy + final double[][] nData = nordsieck.getDataRef(); + for (int i = nData.length - 1; i >= 0; --i) { + final int order = i + 2; + final double[] nDataI = nData[i]; + final double power = FastMath.pow(normalizedAbscissa, order); + for (int j = 0; j < nDataI.length; ++j) { + final double d = nDataI[j] * power; + stateVariation[j] += d; + interpolatedDerivatives[j] += order * d; + } + } + + for (int j = 0; j < currentState.length; ++j) { + stateVariation[j] += scaled[j] * normalizedAbscissa; + interpolatedState[j] = currentState[j] + stateVariation[j]; + interpolatedDerivatives[j] = + (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x; + } + + } + + /** {@inheritDoc} */ + @Override + public void writeExternal(final ObjectOutput out) + throws IOException { + + // save the state of the base class + writeBaseExternal(out); + + // save the local attributes + out.writeDouble(scalingH); + out.writeDouble(referenceTime); + + final int n = (currentState == null) ? -1 : currentState.length; + if (scaled == null) { + out.writeBoolean(false); + } else { + out.writeBoolean(true); + for (int j = 0; j < n; ++j) { + out.writeDouble(scaled[j]); + } + } + + if (nordsieck == null) { + out.writeBoolean(false); + } else { + out.writeBoolean(true); + out.writeObject(nordsieck); + } + + // we don't save state variation, it will be recomputed + + } + + /** {@inheritDoc} */ + @Override + public void readExternal(final ObjectInput in) + throws IOException, ClassNotFoundException { + + // read the base class + final double t = readBaseExternal(in); + + // read the local attributes + scalingH = in.readDouble(); + referenceTime = in.readDouble(); + + final int n = (currentState == null) ? -1 : currentState.length; + final boolean hasScaled = in.readBoolean(); + if (hasScaled) { + scaled = new double[n]; + for (int j = 0; j < n; ++j) { + scaled[j] = in.readDouble(); + } + } else { + scaled = null; + } + + final boolean hasNordsieck = in.readBoolean(); + if (hasNordsieck) { + nordsieck = (Array2DRowRealMatrix) in.readObject(); + } else { + nordsieck = null; + } + + if (hasScaled && hasNordsieck) { + // we can now set the interpolated time and state + stateVariation = new double[n]; + setInterpolatedTime(t); + } else { + stateVariation = null; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/StepHandler.java b/src/main/java/org/apache/commons/math3/ode/sampling/StepHandler.java new file mode 100644 index 0000000..f4c63df --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/StepHandler.java @@ -0,0 +1,75 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.exception.MaxCountExceededException; + + +/** + * This interface represents a handler that should be called after + * each successful step. + * + * <p>The ODE integrators compute the evolution of the state vector at + * some grid points that depend on their own internal algorithm. Once + * they have found a new grid point (possibly after having computed + * several evaluation of the derivative at intermediate points), they + * provide it to objects implementing this interface. These objects + * typically either ignore the intermediate steps and wait for the + * last one, store the points in an ephemeris, or forward them to + * specialized processing or output methods.</p> + * + * @see org.apache.commons.math3.ode.FirstOrderIntegrator + * @see org.apache.commons.math3.ode.SecondOrderIntegrator + * @see StepInterpolator + * @since 1.2 + */ + +public interface StepHandler { + + /** Initialize step handler at the start of an ODE integration. + * <p> + * This method is called once at the start of the integration. It + * may be used by the step handler to initialize some internal data + * if needed. + * </p> + * @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 + */ + void init(double t0, double[] y0, double t); + + /** + * Handle the last accepted step + * @param interpolator interpolator for the last accepted step. For + * efficiency purposes, the various integrators reuse the same + * object on each call, so if the instance wants to keep it across + * all calls (for example to provide at the end of the integration a + * continuous model valid throughout the integration range, as the + * {@link org.apache.commons.math3.ode.ContinuousOutputModel + * ContinuousOutputModel} class does), it should build a local copy + * using the clone method of the interpolator and store this copy. + * Keeping only a reference to the interpolator and reusing it will + * result in unpredictable behavior (potentially crashing the application). + * @param isLast true if the step is the last one + * @exception MaxCountExceededException if the interpolator throws one because + * the number of functions evaluations is exceeded + */ + void handleStep(StepInterpolator interpolator, boolean isLast) + throws MaxCountExceededException; + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java new file mode 100644 index 0000000..5d27bf2 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java @@ -0,0 +1,181 @@ +/* + * 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.sampling; + +import java.io.Externalizable; + +import org.apache.commons.math3.exception.MaxCountExceededException; + +/** This interface represents an interpolator over the last step + * during an ODE integration. + * + * <p>The various ODE integrators provide objects implementing this + * interface to the step handlers. These objects are often custom + * objects tightly bound to the integrator internal algorithms. The + * handlers can use these objects to retrieve the state vector at + * intermediate times between the previous and the current grid points + * (this feature is often called dense output).</p> + * <p>One important thing to note is that the step handlers may be so + * tightly bound to the integrators that they often share some internal + * state arrays. This imply that one should <em>never</em> use a direct + * reference to a step interpolator outside of the step handler, either + * for future use or for use in another thread. If such a need arise, the + * step interpolator <em>must</em> be copied using the dedicated + * {@link #copy()} method. + * </p> + * + * @see org.apache.commons.math3.ode.FirstOrderIntegrator + * @see org.apache.commons.math3.ode.SecondOrderIntegrator + * @see StepHandler + * @since 1.2 + */ + +public interface StepInterpolator extends Externalizable { + + /** + * Get the previous grid point time. + * @return previous grid point time + */ + double getPreviousTime(); + + /** + * Get the current grid point time. + * @return current grid point time + */ + double getCurrentTime(); + + /** + * Get the time of the interpolated point. + * If {@link #setInterpolatedTime} has not been called, it returns + * the current grid point time. + * @return interpolation point time + */ + double getInterpolatedTime(); + + /** + * Set the time of the interpolated point. + * <p>Setting the time outside of the current step is now allowed, but + * should be used with care since the accuracy of the interpolator will + * probably be very poor far from this step. This allowance has been + * added to simplify implementation of search algorithms near the + * step endpoints.</p> + * <p>Setting the time changes the instance internal state. This includes + * the internal arrays returned in {@link #getInterpolatedState()}, + * {@link #getInterpolatedDerivatives()}, {@link + * #getInterpolatedSecondaryState(int)} and {@link + * #getInterpolatedSecondaryDerivatives(int)}. So if their content must be preserved + * across several calls, user must copy them.</p> + * @param time time of the interpolated point + * @see #getInterpolatedState() + * @see #getInterpolatedDerivatives() + * @see #getInterpolatedSecondaryState(int) + * @see #getInterpolatedSecondaryDerivatives(int) + */ + void setInterpolatedTime(double time); + + /** + * Get the state vector of the interpolated point. + * <p>The returned vector is a reference to a reused array, so + * it should not be modified and it should be copied if it needs + * to be preserved across several calls to the associated + * {@link #setInterpolatedTime(double)} method.</p> + * @return state vector at time {@link #getInterpolatedTime} + * @see #getInterpolatedDerivatives() + * @see #getInterpolatedSecondaryState(int) + * @see #getInterpolatedSecondaryDerivatives(int) + * @see #setInterpolatedTime(double) + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + double[] getInterpolatedState() throws MaxCountExceededException; + + /** + * Get the derivatives of the state vector of the interpolated point. + * <p>The returned vector is a reference to a reused array, so + * it should not be modified and it should be copied if it needs + * to be preserved across several calls to the associated + * {@link #setInterpolatedTime(double)} method.</p> + * @return derivatives of the state vector at time {@link #getInterpolatedTime} + * @see #getInterpolatedState() + * @see #getInterpolatedSecondaryState(int) + * @see #getInterpolatedSecondaryDerivatives(int) + * @see #setInterpolatedTime(double) + * @since 2.0 + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + double[] getInterpolatedDerivatives() throws MaxCountExceededException; + + /** Get the interpolated secondary state corresponding to the secondary equations. + * <p>The returned vector is a reference to a reused array, so + * it should not be modified and it should be copied if it needs + * to be preserved across several calls to the associated + * {@link #setInterpolatedTime(double)} method.</p> + * @param index index of the secondary set, as returned by {@link + * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( + * org.apache.commons.math3.ode.SecondaryEquations) + * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)} + * @return interpolated secondary state at the current interpolation date + * @see #getInterpolatedState() + * @see #getInterpolatedDerivatives() + * @see #getInterpolatedSecondaryDerivatives(int) + * @see #setInterpolatedTime(double) + * @since 3.0 + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + double[] getInterpolatedSecondaryState(int index) throws MaxCountExceededException; + + /** Get the interpolated secondary derivatives corresponding to the secondary equations. + * <p>The returned vector is a reference to a reused array, so + * it should not be modified and it should be copied if it needs + * to be preserved across several calls.</p> + * @param index index of the secondary set, as returned by {@link + * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( + * org.apache.commons.math3.ode.SecondaryEquations) + * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)} + * @return interpolated secondary derivatives at the current interpolation date + * @see #getInterpolatedState() + * @see #getInterpolatedDerivatives() + * @see #getInterpolatedSecondaryState(int) + * @see #setInterpolatedTime(double) + * @since 3.0 + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + double[] getInterpolatedSecondaryDerivatives(int index) throws MaxCountExceededException; + + /** Check if the natural integration direction is forward. + * <p>This method provides the integration direction as specified by + * the integrator itself, it avoid some nasty problems in + * degenerated cases like null steps due to cancellation at step + * initialization, step control or discrete events + * triggering.</p> + * @return true if the integration variable (time) increases during + * integration + */ + boolean isForward(); + + /** Copy the instance. + * <p>The copied instance is guaranteed to be independent from the + * original one. Both can be used with different settings for + * interpolated time without any side effect.</p> + * @return a deep copy of the instance, which can be used independently. + * @see #setInterpolatedTime(double) + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + * during step finalization + */ + StepInterpolator copy() throws MaxCountExceededException; + +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizer.java b/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizer.java new file mode 100644 index 0000000..307c0d9 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizer.java @@ -0,0 +1,300 @@ +/* + * 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.sampling; + +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; + +/** + * This class wraps an object implementing {@link FixedStepHandler} + * into a {@link StepHandler}. + + * <p>This wrapper allows to use fixed step handlers with general + * integrators which cannot guaranty their integration steps will + * remain constant and therefore only accept general step + * handlers.</p> + * + * <p>The stepsize used is selected at construction time. The {@link + * FixedStepHandler#handleStep handleStep} method of the underlying + * {@link FixedStepHandler} object is called at normalized times. The + * normalized times can be influenced by the {@link StepNormalizerMode} and + * {@link StepNormalizerBounds}.</p> + * + * <p>There is no constraint on the integrator, it can use any time step + * it needs (time steps longer or shorter than the fixed time step and + * non-integer ratios are all allowed).</p> + * + * <p> + * <table border="1" align="center"> + * <tr BGCOLOR="#CCCCFF"><td colspan=6><font size="+2">Examples (step size = 0.5)</font></td></tr> + * <tr BGCOLOR="#EEEEFF"><font size="+1"><td>Start time</td><td>End time</td> + * <td>Direction</td><td>{@link StepNormalizerMode Mode}</td> + * <td>{@link StepNormalizerBounds Bounds}</td><td>Output</td></font></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.8, 1.3, 1.8, 2.3, 2.8</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.3, 0.8, 1.3, 1.8, 2.3, 2.8</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.8, 1.3, 1.8, 2.3, 2.8, 3.1</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.3, 0.8, 1.3, 1.8, 2.3, 2.8, 3.1</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.1</td></tr> + * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.1</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.6, 2.1, 1.6, 1.1, 0.6</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.1, 2.6, 2.1, 1.6, 1.1, 0.6</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.6, 2.1, 1.6, 1.1, 0.6, 0.3</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.1, 2.6, 2.1, 1.6, 1.1, 0.6, 0.3</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.1, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.3</td></tr> + * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.1, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.3</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr> + * </table> + * </p> + * + * @see StepHandler + * @see FixedStepHandler + * @see StepNormalizerMode + * @see StepNormalizerBounds + * @since 1.2 + */ + +public class StepNormalizer implements StepHandler { + /** Fixed time step. */ + private double h; + + /** Underlying step handler. */ + private final FixedStepHandler handler; + + /** First step time. */ + private double firstTime; + + /** Last step time. */ + private double lastTime; + + /** Last state vector. */ + private double[] lastState; + + /** Last derivatives vector. */ + private double[] lastDerivatives; + + /** Integration direction indicator. */ + private boolean forward; + + /** The step normalizer bounds settings to use. */ + private final StepNormalizerBounds bounds; + + /** The step normalizer mode to use. */ + private final StepNormalizerMode mode; + + /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT} + * mode, and {@link StepNormalizerBounds#FIRST FIRST} bounds setting, for + * backwards compatibility. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + */ + public StepNormalizer(final double h, final FixedStepHandler handler) { + this(h, handler, StepNormalizerMode.INCREMENT, + StepNormalizerBounds.FIRST); + } + + /** Simple constructor. Uses {@link StepNormalizerBounds#FIRST FIRST} + * bounds setting. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + * @param mode step normalizer mode to use + * @since 3.0 + */ + public StepNormalizer(final double h, final FixedStepHandler handler, + final StepNormalizerMode mode) { + this(h, handler, mode, StepNormalizerBounds.FIRST); + } + + /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT} + * mode. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + * @param bounds step normalizer bounds setting to use + * @since 3.0 + */ + public StepNormalizer(final double h, final FixedStepHandler handler, + final StepNormalizerBounds bounds) { + this(h, handler, StepNormalizerMode.INCREMENT, bounds); + } + + /** Simple constructor. + * @param h fixed time step (sign is not used) + * @param handler fixed time step handler to wrap + * @param mode step normalizer mode to use + * @param bounds step normalizer bounds setting to use + * @since 3.0 + */ + public StepNormalizer(final double h, final FixedStepHandler handler, + final StepNormalizerMode mode, + final StepNormalizerBounds bounds) { + this.h = FastMath.abs(h); + this.handler = handler; + this.mode = mode; + this.bounds = bounds; + firstTime = Double.NaN; + lastTime = Double.NaN; + lastState = null; + lastDerivatives = null; + forward = true; + } + + /** {@inheritDoc} */ + public void init(double t0, double[] y0, double t) { + + firstTime = Double.NaN; + lastTime = Double.NaN; + lastState = null; + lastDerivatives = null; + forward = true; + + // initialize the underlying handler + handler.init(t0, y0, t); + + } + + /** + * Handle the last accepted step + * @param interpolator interpolator for the last accepted step. For + * efficiency purposes, the various integrators reuse the same + * object on each call, so if the instance wants to keep it across + * all calls (for example to provide at the end of the integration a + * continuous model valid throughout the integration range), it + * should build a local copy using the clone method and store this + * copy. + * @param isLast true if the step is the last one + * @exception MaxCountExceededException if the interpolator throws one because + * the number of functions evaluations is exceeded + */ + public void handleStep(final StepInterpolator interpolator, final boolean isLast) + throws MaxCountExceededException { + // The first time, update the last state with the start information. + if (lastState == null) { + firstTime = interpolator.getPreviousTime(); + lastTime = interpolator.getPreviousTime(); + interpolator.setInterpolatedTime(lastTime); + lastState = interpolator.getInterpolatedState().clone(); + lastDerivatives = interpolator.getInterpolatedDerivatives().clone(); + + // Take the integration direction into account. + forward = interpolator.getCurrentTime() >= lastTime; + if (!forward) { + h = -h; + } + } + + // Calculate next normalized step time. + double nextTime = (mode == StepNormalizerMode.INCREMENT) ? + lastTime + h : + (FastMath.floor(lastTime / h) + 1) * h; + if (mode == StepNormalizerMode.MULTIPLES && + Precision.equals(nextTime, lastTime, 1)) { + nextTime += h; + } + + // Process normalized steps as long as they are in the current step. + boolean nextInStep = isNextInStep(nextTime, interpolator); + while (nextInStep) { + // Output the stored previous step. + doNormalizedStep(false); + + // Store the next step as last step. + storeStep(interpolator, nextTime); + + // Move on to the next step. + nextTime += h; + nextInStep = isNextInStep(nextTime, interpolator); + } + + if (isLast) { + // There will be no more steps. The stored one should be given to + // the handler. We may have to output one more step. Only the last + // one of those should be flagged as being the last. + boolean addLast = bounds.lastIncluded() && + lastTime != interpolator.getCurrentTime(); + doNormalizedStep(!addLast); + if (addLast) { + storeStep(interpolator, interpolator.getCurrentTime()); + doNormalizedStep(true); + } + } + } + + /** + * Returns a value indicating whether the next normalized time is in the + * current step. + * @param nextTime the next normalized time + * @param interpolator interpolator for the last accepted step, to use to + * get the end time of the current step + * @return value indicating whether the next normalized time is in the + * current step + */ + private boolean isNextInStep(double nextTime, + StepInterpolator interpolator) { + return forward ? + nextTime <= interpolator.getCurrentTime() : + nextTime >= interpolator.getCurrentTime(); + } + + /** + * Invokes the underlying step handler for the current normalized step. + * @param isLast true if the step is the last one + */ + private void doNormalizedStep(boolean isLast) { + if (!bounds.firstIncluded() && firstTime == lastTime) { + return; + } + handler.handleStep(lastTime, lastState, lastDerivatives, isLast); + } + + /** Stores the interpolated information for the given time in the current + * state. + * @param interpolator interpolator for the last accepted step, to use to + * get the interpolated information + * @param t the time for which to store the interpolated information + * @exception MaxCountExceededException if the interpolator throws one because + * the number of functions evaluations is exceeded + */ + private void storeStep(StepInterpolator interpolator, double t) + throws MaxCountExceededException { + lastTime = t; + interpolator.setInterpolatedTime(lastTime); + System.arraycopy(interpolator.getInterpolatedState(), 0, + lastState, 0, lastState.length); + System.arraycopy(interpolator.getInterpolatedDerivatives(), 0, + lastDerivatives, 0, lastDerivatives.length); + } +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizerBounds.java b/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizerBounds.java new file mode 100644 index 0000000..ca35e82 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizerBounds.java @@ -0,0 +1,84 @@ +/* + * 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.sampling; + +/** {@link StepNormalizer Step normalizer} bounds settings. They influence + * whether the underlying fixed step size step handler is called for the first + * and last points. Note that if the last point coincides with a normalized + * point, then the underlying fixed step size step handler is always called, + * regardless of these settings. + * @see FieldStepNormalizer + * @see StepNormalizer + * @see StepNormalizerMode + * @since 3.0 + */ +public enum StepNormalizerBounds { + /** Do not include the first and last points. */ + NEITHER(false, false), + + /** Include the first point, but not the last point. */ + FIRST(true, false), + + /** Include the last point, but not the first point. */ + LAST(false, true), + + /** Include both the first and last points. */ + BOTH(true, true); + + /** Whether the first point should be passed to the underlying fixed + * step size step handler. + */ + private final boolean first; + + /** Whether the last point should be passed to the underlying fixed + * step size step handler. + */ + private final boolean last; + + /** + * Simple constructor. + * @param first Whether the first point should be passed to the + * underlying fixed step size step handler. + * @param last Whether the last point should be passed to the + * underlying fixed step size step handler. + */ + StepNormalizerBounds(final boolean first, final boolean last) { + this.first = first; + this.last = last; + } + + /** + * Returns a value indicating whether the first point should be passed + * to the underlying fixed step size step handler. + * @return value indicating whether the first point should be passed + * to the underlying fixed step size step handler. + */ + public boolean firstIncluded() { + return first; + } + + /** + * Returns a value indicating whether the last point should be passed + * to the underlying fixed step size step handler. + * @return value indicating whether the last point should be passed + * to the underlying fixed step size step handler. + */ + public boolean lastIncluded() { + return last; + } +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizerMode.java b/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizerMode.java new file mode 100644 index 0000000..c6f4c69 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/StepNormalizerMode.java @@ -0,0 +1,70 @@ +/* + * 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.sampling; + + +/** {@link StepNormalizer Step normalizer} modes. Determines how the step size + * is interpreted. + * @see FieldStepNormalizer + * @see StepNormalizer + * @see StepNormalizerBounds + * @since 3.0 + */ +public enum StepNormalizerMode { + /** + * Steps are fixed increments of the start value. In other words, they + * are relative to the start value. + * + * <p>If the integration start time is t0, then the points handled by + * the underlying fixed step size step handler are t0 (depending on + * the {@link StepNormalizerBounds bounds settings}), t0+h, t0+2h, ...</p> + * + * <p>If the integration range is an integer multiple of the step size + * (h), then the last point handled will be the end point of the + * integration (tend). If not, the last point may be the end point + * tend, or it may be a point belonging to the interval [tend - h ; + * tend], depending on the {@link StepNormalizerBounds bounds settings}. + * </p> + * + * @see StepNormalizer + * @see StepNormalizerBounds + */ + INCREMENT, + + /** Steps are multiples of a fixed value. In other words, they are + * relative to the first multiple of the step size that is encountered + * after the start value. + * + * <p>If the integration start time is t0, and the first multiple of + * the fixed step size that is encountered is t1, then the points + * handled by the underlying fixed step size step handler are t0 + * (depending on the {@link StepNormalizerBounds bounds settings}), t1, + * t1+h, t1+2h, ...</p> + * + * <p>If the end point of the integration range (tend) is an integer + * multiple of the step size (h) added to t1, then the last point + * handled will be the end point of the integration (tend). If not, + * the last point may be the end point tend, or it may be a point + * belonging to the interval [tend - h ; tend], depending on the + * {@link StepNormalizerBounds bounds settings}.</p> + * + * @see StepNormalizer + * @see StepNormalizerBounds + */ + MULTIPLES; +} diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/package-info.java b/src/main/java/org/apache/commons/math3/ode/sampling/package-info.java new file mode 100644 index 0000000..29c65e0 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/sampling/package-info.java @@ -0,0 +1,60 @@ +/* + * 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. + */ +/** + * + * <p> + * This package provides classes to handle sampling steps during + * Ordinary Differential Equations integration. + * </p> + * + * <p> + * In addition to computing the evolution of the state vector at some grid points, all + * ODE integrators also build up interpolation models of this evolution <em>inside</em> the + * last computed step. If users are interested in these interpolators, they can register a + * {@link org.apache.commons.math3.ode.sampling.StepHandler StepHandler} instance using the + * {@link org.apache.commons.math3.ode.FirstOrderIntegrator#addStepHandler addStepHandler} + * method which is supported by all integrators. The integrator will call this instance + * at the end of each accepted step and provide it the interpolator. The user can do + * whatever he wants with this interpolator, which computes both the state and its + * time-derivative. A typical use of step handler is to provide some output to monitor + * the integration process. + * </p> + * + * <p> + * In a sense, this is a kind of Inversion Of Control: rather than having the master + * application driving the slave integrator by providing the target end value for + * the free variable, we get a master integrator scheduling the free variable + * evolution and calling the slave application callbacks that were registered at + * configuration time. + * </p> + * + * <p> + * Since some integrators may use variable step size, the generic {@link + * org.apache.commons.math3.ode.sampling.StepHandler StepHandler} interface can be called + * either at regular or irregular rate. This interface allows to navigate to any location + * within the last computed step, thanks to the provided {@link + * org.apache.commons.math3.ode.sampling.StepInterpolator StepInterpolator} object. + * If regular output is desired (for example in order to write an ephemeris file), then + * the simpler {@link org.apache.commons.math3.ode.sampling.FixedStepHandler FixedStepHandler} + * interface can be used. Objects implementing this interface should be wrapped within a + * {@link org.apache.commons.math3.ode.sampling.StepNormalizer StepNormalizer} instance + * in order to be registered to the integrator. + * </p> + * + * + */ +package org.apache.commons.math3.ode.sampling; |