diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java | 482 |
1 files changed, 482 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java new file mode 100644 index 0000000..44a8b81 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java @@ -0,0 +1,482 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math3.ode; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NoBracketingException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.linear.Array2DRowFieldMatrix; +import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeFieldIntegrator; +import org.apache.commons.math3.ode.nonstiff.DormandPrince853FieldIntegrator; +import org.apache.commons.math3.ode.sampling.FieldStepHandler; +import org.apache.commons.math3.ode.sampling.FieldStepInterpolator; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.MathUtils; + +/** + * This class is the base class for multistep integrators for Ordinary Differential Equations. + * + * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as: + * + * <pre> + * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative + * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative + * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative + * ... + * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative + * </pre> + * + * <p>Rather than storing several previous steps separately, this implementation uses the Nordsieck + * vector with higher degrees scaled derivatives all taken at the same step (y<sub>n</sub>, + * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as: + * + * <pre> + * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> + * </pre> + * + * (we omit the k index in the notation for clarity) + * + * <p>Multistep integrators with Nordsieck representation are highly sensitive to large step changes + * because when the step is multiplied by factor a, the k<sup>th</sup> component of the Nordsieck + * vector is multiplied by a<sup>k</sup> and the last components are the least accurate ones. The + * default max growth factor is therefore set to a quite low value: 2<sup>1/order</sup>. + * + * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthFieldIntegrator + * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonFieldIntegrator + * @param <T> the type of the field elements + * @since 3.6 + */ +public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>> + extends AdaptiveStepsizeFieldIntegrator<T> { + + /** First scaled derivative (h y'). */ + protected T[] scaled; + + /** + * Nordsieck matrix of the higher scaled derivatives. + * + * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>) + */ + protected Array2DRowFieldMatrix<T> nordsieck; + + /** Starter integrator. */ + private FirstOrderFieldIntegrator<T> starter; + + /** Number of steps of the multistep method (excluding the one being computed). */ + private final int nSteps; + + /** Stepsize control exponent. */ + private double exp; + + /** Safety factor for stepsize control. */ + private double safety; + + /** Minimal reduction factor for stepsize control. */ + private double minReduction; + + /** Maximal growth factor for stepsize control. */ + private double maxGrowth; + + /** + * Build a multistep integrator with the given stepsize bounds. + * + * <p>The default starter integrator is set to the {@link DormandPrince853FieldIntegrator + * Dormand-Prince 8(5,3)} integrator with some defaults settings. + * + * <p>The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. + * + * @param field field to which the time and state vector elements belong + * @param name name of the method + * @param nSteps number of steps of the multistep method (excluding the one being computed) + * @param order order of the method + * @param minStep minimal step (must be positive even for backward integration), the last step + * can be smaller than this + * @param maxStep maximal step (must be positive even for backward integration) + * @param scalAbsoluteTolerance allowed absolute error + * @param scalRelativeTolerance allowed relative error + * @exception NumberIsTooSmallException if number of steps is smaller than 2 + */ + protected MultistepFieldIntegrator( + final Field<T> field, + final String name, + final int nSteps, + final int order, + final double minStep, + final double maxStep, + final double scalAbsoluteTolerance, + final double scalRelativeTolerance) + throws NumberIsTooSmallException { + + super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); + + if (nSteps < 2) { + throw new NumberIsTooSmallException( + LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS, + nSteps, + 2, + true); + } + + starter = + new DormandPrince853FieldIntegrator<T>( + field, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); + this.nSteps = nSteps; + + exp = -1.0 / order; + + // set the default values of the algorithm control parameters + setSafety(0.9); + setMinReduction(0.2); + setMaxGrowth(FastMath.pow(2.0, -exp)); + } + + /** + * Build a multistep integrator with the given stepsize bounds. + * + * <p>The default starter integrator is set to the {@link DormandPrince853FieldIntegrator + * Dormand-Prince 8(5,3)} integrator with some defaults settings. + * + * <p>The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. + * + * @param field field to which the time and state vector elements belong + * @param name name of the method + * @param nSteps number of steps of the multistep method (excluding the one being computed) + * @param order order of the method + * @param minStep minimal step (must be positive even for backward integration), the last step + * can be smaller than this + * @param maxStep maximal step (must be positive even for backward integration) + * @param vecAbsoluteTolerance allowed absolute error + * @param vecRelativeTolerance allowed relative error + */ + protected MultistepFieldIntegrator( + final Field<T> field, + final String name, + final int nSteps, + final int order, + final double minStep, + final double maxStep, + final double[] vecAbsoluteTolerance, + final double[] vecRelativeTolerance) { + super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); + starter = + new DormandPrince853FieldIntegrator<T>( + field, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); + this.nSteps = nSteps; + + exp = -1.0 / order; + + // set the default values of the algorithm control parameters + setSafety(0.9); + setMinReduction(0.2); + setMaxGrowth(FastMath.pow(2.0, -exp)); + } + + /** + * Get the starter integrator. + * + * @return starter integrator + */ + public FirstOrderFieldIntegrator<T> getStarterIntegrator() { + return starter; + } + + /** + * Set the starter integrator. + * + * <p>The various step and event handlers for this starter integrator will be managed + * automatically by the multi-step integrator. Any user configuration for these elements will be + * cleared before use. + * + * @param starterIntegrator starter integrator + */ + public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) { + this.starter = starterIntegrator; + } + + /** + * Start the integration. + * + * <p>This method computes one step using the underlying starter integrator, and initializes the + * Nordsieck vector at step start. The starter integrator purpose is only to establish initial + * conditions, it does not really change time by itself. The top level multistep integrator + * remains in charge of handling time propagation and events handling as it will starts its own + * computation right from the beginning. In a sense, the starter integrator can be seen as a + * dummy one and so it will never trigger any user event nor call any user step handler. + * + * @param equations complete set of differential equations to integrate + * @param initialState initial state (time, primary and secondary state vectors) + * @param t target time for the integration (can be set to a value smaller than <code>t0</code> + * for backward integration) + * @exception DimensionMismatchException if arrays dimension do not match equations settings + * @exception NumberIsTooSmallException if integration step is too small + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + * @exception NoBracketingException if the location of an event cannot be bracketed + */ + protected void start( + final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t) + throws DimensionMismatchException, + NumberIsTooSmallException, + MaxCountExceededException, + NoBracketingException { + + // make sure NO user event nor user step handler is triggered, + // this is the task of the top level integrator, not the task + // of the starter integrator + starter.clearEventHandlers(); + starter.clearStepHandlers(); + + // set up one specific step handler to extract initial Nordsieck vector + starter.addStepHandler( + new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2)); + + // start integration, expecting a InitializationCompletedMarkerException + try { + + starter.integrate(equations, initialState, t); + + // we should not reach this step + throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY); + + } catch (InitializationCompletedMarkerException icme) { // NOPMD + // this is the expected nominal interruption of the start integrator + + // count the evaluations used by the starter + getEvaluationsCounter().increment(starter.getEvaluations()); + } + + // remove the specific step handler + starter.clearStepHandlers(); + } + + /** + * Initialize the high order scaled derivatives at step start. + * + * @param h step size to use for scaling + * @param t first steps times + * @param y first steps states + * @param yDot first steps derivatives + * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>, h<sup>3</sup>/6 + * y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>) + */ + protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives( + final T h, final T[] t, final T[][] y, final T[][] yDot); + + /** + * Get the minimal reduction factor for stepsize control. + * + * @return minimal reduction factor + */ + public double getMinReduction() { + return minReduction; + } + + /** + * Set the minimal reduction factor for stepsize control. + * + * @param minReduction minimal reduction factor + */ + public void setMinReduction(final double minReduction) { + this.minReduction = minReduction; + } + + /** + * Get the maximal growth factor for stepsize control. + * + * @return maximal growth factor + */ + public double getMaxGrowth() { + return maxGrowth; + } + + /** + * Set the maximal growth factor for stepsize control. + * + * @param maxGrowth maximal growth factor + */ + public void setMaxGrowth(final double maxGrowth) { + this.maxGrowth = maxGrowth; + } + + /** + * Get the safety factor for stepsize control. + * + * @return safety factor + */ + public double getSafety() { + return safety; + } + + /** + * Set the safety factor for stepsize control. + * + * @param safety safety factor + */ + public void setSafety(final double safety) { + this.safety = safety; + } + + /** + * Get the number of steps of the multistep method (excluding the one being computed). + * + * @return number of steps of the multistep method (excluding the one being computed) + */ + public int getNSteps() { + return nSteps; + } + + /** + * 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. + * + * @param newStepSize new step size to use in the scaled and Nordsieck arrays + */ + protected void rescale(final T newStepSize) { + + final T ratio = newStepSize.divide(getStepSize()); + for (int i = 0; i < scaled.length; ++i) { + scaled[i] = scaled[i].multiply(ratio); + } + + final T[][] nData = nordsieck.getDataRef(); + T power = ratio; + for (int i = 0; i < nData.length; ++i) { + power = power.multiply(ratio); + final T[] nDataI = nData[i]; + for (int j = 0; j < nDataI.length; ++j) { + nDataI[j] = nDataI[j].multiply(power); + } + } + + setStepSize(newStepSize); + } + + /** + * Compute step grow/shrink factor according to normalized error. + * + * @param error normalized error of the current step + * @return grow/shrink factor for next step + */ + protected T computeStepGrowShrinkFactor(final T error) { + return MathUtils.min( + error.getField().getZero().add(maxGrowth), + MathUtils.max( + error.getField().getZero().add(minReduction), + error.pow(exp).multiply(safety))); + } + + /** Specialized step handler storing the first step. */ + private class FieldNordsieckInitializer implements FieldStepHandler<T> { + + /** Equation mapper. */ + private final FieldEquationsMapper<T> mapper; + + /** Steps counter. */ + private int count; + + /** Saved start. */ + private FieldODEStateAndDerivative<T> savedStart; + + /** First steps times. */ + private final T[] t; + + /** First steps states. */ + private final T[][] y; + + /** First steps derivatives. */ + private final T[][] yDot; + + /** + * Simple constructor. + * + * @param mapper equation mapper + * @param nbStartPoints number of start points (including the initial point) + */ + FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) { + this.mapper = mapper; + this.count = 0; + this.t = MathArrays.buildArray(getField(), nbStartPoints); + this.y = MathArrays.buildArray(getField(), nbStartPoints, -1); + this.yDot = MathArrays.buildArray(getField(), nbStartPoints, -1); + } + + /** {@inheritDoc} */ + public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) + throws MaxCountExceededException { + + if (count == 0) { + // first step, we need to store also the point at the beginning of the step + final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState(); + savedStart = prev; + t[count] = prev.getTime(); + y[count] = mapper.mapState(prev); + yDot[count] = mapper.mapDerivative(prev); + } + + // store the point at the end of the step + ++count; + final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState(); + t[count] = curr.getTime(); + y[count] = mapper.mapState(curr); + yDot[count] = mapper.mapDerivative(curr); + + if (count == t.length - 1) { + + // this was the last point we needed, we can compute the derivatives + setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1)); + + // first scaled derivative + scaled = MathArrays.buildArray(getField(), yDot[0].length); + for (int j = 0; j < scaled.length; ++j) { + scaled[j] = yDot[0][j].multiply(getStepSize()); + } + + // higher order derivatives + nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot); + + // stop the integrator now that all needed steps have been handled + setStepStart(savedStart); + throw new InitializationCompletedMarkerException(); + } + } + + /** {@inheritDoc} */ + public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) { + // nothing to do + } + } + + /** Marker exception used ONLY to stop the starter integrator after first step. */ + private static class InitializationCompletedMarkerException extends RuntimeException { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -1914085471038046418L; + + /** Simple constructor. */ + InitializationCompletedMarkerException() { + super((Throwable) null); + } + } +} |