    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.ode;
    019
    020    import org.apache.commons.math3.exception.DimensionMismatchException;
    021    import org.apache.commons.math3.exception.MaxCountExceededException;
    022    import org.apache.commons.math3.exception.NoBracketingException;
    023    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
    026    import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
    027    import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
    028    import org.apache.commons.math3.ode.sampling.StepHandler;
    029    import org.apache.commons.math3.ode.sampling.StepInterpolator;
    030    import org.apache.commons.math3.util.FastMath;
    031
    032    /**
    033     * This class is the base class for multistep integrators for Ordinary
    034     * Differential Equations.
    035     * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
    036     * <pre>
    037     * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
    038     * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
    039     * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
    040     * ...
    041     * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
    042     * </pre></p>
    043     * <p>Rather than storing several previous steps separately, this implementation uses
    044     * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
    045     * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
    046     * <pre>
    047     * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
    048     * </pre>
    049     * (we omit the k index in the notation for clarity)</p>
    050     * <p>
    051     * Multistep integrators with Nordsieck representation are highly sensitive to
    052     * large step changes because when the step is multiplied by factor a, the
    053     * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
    054     * and the last components are the least accurate ones. The default max growth
    055     * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
    056     * </p>
    057     *
    058     * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator
    059     * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator
    060     * @version $Id: MultistepIntegrator.java 1421448 2012-12-13 19:45:57Z tn $
    061     * @since 2.0
    062     */
    063    public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
    064
    065        /** First scaled derivative (h y'). */
    066        protected double[] scaled;
    067
    068        /** Nordsieck matrix of the higher scaled derivatives.
    069         * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
    070         */
    071        protected Array2DRowRealMatrix nordsieck;
    072
    073        /** Starter integrator. */
    074        private FirstOrderIntegrator starter;
    075
    076        /** Number of steps of the multistep method (excluding the one being computed). */
    077        private final int nSteps;
    078
    079        /** Stepsize control exponent. */
    080        private double exp;
    081
    082        /** Safety factor for stepsize control. */
    083        private double safety;
    084
    085        /** Minimal reduction factor for stepsize control. */
    086        private double minReduction;
    087
    088        /** Maximal growth factor for stepsize control. */
    089        private double maxGrowth;
    090
    091        /**
    092         * Build a multistep integrator with the given stepsize bounds.
    093         * <p>The default starter integrator is set to the {@link
    094         * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
    095         * some defaults settings.</p>
    096         * <p>
    097         * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
    098         * </p>
    099         * @param name name of the method
    100         * @param nSteps number of steps of the multistep method
    101         * (excluding the one being computed)
    102         * @param order order of the method
    103         * @param minStep minimal step (must be positive even for backward
    104         * integration), the last step can be smaller than this
    105         * @param maxStep maximal step (must be positive even for backward
    106         * integration)
    107         * @param scalAbsoluteTolerance allowed absolute error
    108         * @param scalRelativeTolerance allowed relative error
    109         * @exception NumberIsTooSmallException if number of steps is smaller than 2
    110         */
    111        protected MultistepIntegrator(final String name, final int nSteps,
    112                                      final int order,
    113                                      final double minStep, final double maxStep,
    114                                      final double scalAbsoluteTolerance,
    115                                      final double scalRelativeTolerance)
    116            throws NumberIsTooSmallException {
    117
    118            super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
    119
    120            if (nSteps < 2) {
    121                throw new NumberIsTooSmallException(
    122                      LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
    123                      nSteps, 2, true);
    124            }
    125
    126            starter = new DormandPrince853Integrator(minStep, maxStep,
    127                                                     scalAbsoluteTolerance,
    128                                                     scalRelativeTolerance);
    129            this.nSteps = nSteps;
    130
    131            exp = -1.0 / order;
    132
    133            // set the default values of the algorithm control parameters
    134            setSafety(0.9);
    135            setMinReduction(0.2);
    136            setMaxGrowth(FastMath.pow(2.0, -exp));
    137
    138        }
    139
    140        /**
    141         * Build a multistep integrator with the given stepsize bounds.
    142         * <p>The default starter integrator is set to the {@link
    143         * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
    144         * some defaults settings.</p>
    145         * <p>
    146         * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
    147         * </p>
    148         * @param name name of the method
    149         * @param nSteps number of steps of the multistep method
    150         * (excluding the one being computed)
    151         * @param order order of the method
    152         * @param minStep minimal step (must be positive even for backward
    153         * integration), the last step can be smaller than this
    154         * @param maxStep maximal step (must be positive even for backward
    155         * integration)
    156         * @param vecAbsoluteTolerance allowed absolute error
    157         * @param vecRelativeTolerance allowed relative error
    158         */
    159        protected MultistepIntegrator(final String name, final int nSteps,
    160                                      final int order,
    161                                      final double minStep, final double maxStep,
    162                                      final double[] vecAbsoluteTolerance,
    163                                      final double[] vecRelativeTolerance) {
    164            super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
    165            starter = new DormandPrince853Integrator(minStep, maxStep,
    166                                                     vecAbsoluteTolerance,
    167                                                     vecRelativeTolerance);
    168            this.nSteps = nSteps;
    169
    170            exp = -1.0 / order;
    171
    172            // set the default values of the algorithm control parameters
    173            setSafety(0.9);
    174            setMinReduction(0.2);
    175            setMaxGrowth(FastMath.pow(2.0, -exp));
    176
    177        }
    178
    179        /**
    180         * Get the starter integrator.
    181         * @return starter integrator
    182         */
    183        public ODEIntegrator getStarterIntegrator() {
    184            return starter;
    185        }
    186
    187        /**
    188         * Set the starter integrator.
    189         * <p>The various step and event handlers for this starter integrator
    190         * will be managed automatically by the multi-step integrator. Any
    191         * user configuration for these elements will be cleared before use.</p>
    192         * @param starterIntegrator starter integrator
    193         */
    194        public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
    195            this.starter = starterIntegrator;
    196        }
    197
    198        /** Start the integration.
    199         * <p>This method computes one step using the underlying starter integrator,
    200         * and initializes the Nordsieck vector at step start. The starter integrator
    201         * purpose is only to establish initial conditions, it does not really change
    202         * time by itself. The top level multistep integrator remains in charge of
    203         * handling time propagation and events handling as it will starts its own
    204         * computation right from the beginning. In a sense, the starter integrator
    205         * can be seen as a dummy one and so it will never trigger any user event nor
    206         * call any user step handler.</p>
    207         * @param t0 initial time
    208         * @param y0 initial value of the state vector at t0
    209         * @param t target time for the integration
    210         * (can be set to a value smaller than <code>t0</code> for backward integration)
    211         * @exception DimensionMismatchException if arrays dimension do not match equations settings
    212         * @exception NumberIsTooSmallException if integration step is too small
    213         * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    214         * @exception NoBracketingException if the location of an event cannot be bracketed
    215         */
    216        protected void start(final double t0, final double[] y0, final double t)
    217            throws DimensionMismatchException, NumberIsTooSmallException,
    218                   MaxCountExceededException, NoBracketingException {
    219
    220            // make sure NO user event nor user step handler is triggered,
    221            // this is the task of the top level integrator, not the task
    222            // of the starter integrator
    223            starter.clearEventHandlers();
    224            starter.clearStepHandlers();
    225
    226            // set up one specific step handler to extract initial Nordsieck vector
    227            starter.addStepHandler(new NordsieckInitializer(nSteps, y0.length));
    228
    229            // start integration, expecting a InitializationCompletedMarkerException
    230            try {
    231                starter.integrate(new CountingDifferentialEquations(y0.length),
    232                                  t0, y0, t, new double[y0.length]);
    233            } catch (InitializationCompletedMarkerException icme) { // NOPMD
    234                // this is the expected nominal interruption of the start integrator
    235            }
    236
    237            // remove the specific step handler
    238            starter.clearStepHandlers();
    239
    240        }
    241
    242        /** Initialize the high order scaled derivatives at step start.
    243         * @param h step size to use for scaling
    244         * @param t first steps times
    245         * @param y first steps states
    246         * @param yDot first steps derivatives
    247         * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
    248         * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
    249         */
    250        protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
    251                                                                               final double[][] y,
    252                                                                               final double[][] yDot);
    253
    254        /** Get the minimal reduction factor for stepsize control.
    255         * @return minimal reduction factor
    256         */
    257        public double getMinReduction() {
    258            return minReduction;
    259        }
    260
    261        /** Set the minimal reduction factor for stepsize control.
    262         * @param minReduction minimal reduction factor
    263         */
    264        public void setMinReduction(final double minReduction) {
    265            this.minReduction = minReduction;
    266        }
    267
    268        /** Get the maximal growth factor for stepsize control.
    269         * @return maximal growth factor
    270         */
    271        public double getMaxGrowth() {
    272            return maxGrowth;
    273        }
    274
    275        /** Set the maximal growth factor for stepsize control.
    276         * @param maxGrowth maximal growth factor
    277         */
    278        public void setMaxGrowth(final double maxGrowth) {
    279            this.maxGrowth = maxGrowth;
    280        }
    281
    282        /** Get the safety factor for stepsize control.
    283         * @return safety factor
    284         */
    285        public double getSafety() {
    286          return safety;
    287        }
    288
    289        /** Set the safety factor for stepsize control.
    290         * @param safety safety factor
    291         */
    292        public void setSafety(final double safety) {
    293          this.safety = safety;
    294        }
    295
    296        /** Compute step grow/shrink factor according to normalized error.
    297         * @param error normalized error of the current step
    298         * @return grow/shrink factor for next step
    299         */
    300        protected double computeStepGrowShrinkFactor(final double error) {
    301            return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
    302        }
    303
    304        /** Transformer used to convert the first step to Nordsieck representation. */
    305        public interface NordsieckTransformer {
    306            /** Initialize the high order scaled derivatives at step start.
    307             * @param h step size to use for scaling
    308             * @param t first steps times
    309             * @param y first steps states
    310             * @param yDot first steps derivatives
    311             * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
    312             * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
    313             */
    314            Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
    315                                                                final double[][] y,
    316                                                                final double[][] yDot);
    317        }
    318
    319        /** Specialized step handler storing the first step. */
    320        private class NordsieckInitializer implements StepHandler {
    321
    322            /** Steps counter. */
    323            private int count;
    324
    325            /** First steps times. */
    326            private final double[] t;
    327
    328            /** First steps states. */
    329            private final double[][] y;
    330
    331            /** First steps derivatives. */
    332            private final double[][] yDot;
    333
    334            /** Simple constructor.
    335             * @param nSteps number of steps of the multistep method (excluding the one being computed)
    336             * @param n problem dimension
    337             */
    338            public NordsieckInitializer(final int nSteps, final int n) {
    339                this.count = 0;
    340                this.t     = new double[nSteps];
    341                this.y     = new double[nSteps][n];
    342                this.yDot  = new double[nSteps][n];
    343            }
    344
    345            /** {@inheritDoc} */
    346            public void handleStep(StepInterpolator interpolator, boolean isLast)
    347                throws MaxCountExceededException {
    348
    349                final double prev = interpolator.getPreviousTime();
    350                final double curr = interpolator.getCurrentTime();
    351
    352                if (count == 0) {
    353                    // first step, we need to store also the beginning of the step
    354                    interpolator.setInterpolatedTime(prev);
    355                    t[0] = prev;
    356                    System.arraycopy(interpolator.getInterpolatedState(), 0,
    357                                     y[0],    0, y[0].length);
    358                    System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
    359                                     yDot[0], 0, yDot[0].length);
    360                }
    361
    362                // store the end of the step
    363                ++count;
    364                interpolator.setInterpolatedTime(curr);
    365                t[count] = curr;
    366                System.arraycopy(interpolator.getInterpolatedState(), 0,
    367                                 y[count],    0, y[count].length);
    368                System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
    369                                 yDot[count], 0, yDot[count].length);
    370
    371                if (count == t.length - 1) {
    372
    373                    // this was the last step we needed, we can compute the derivatives
    374                    stepStart = t[0];
    375                    stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
    376
    377                    // first scaled derivative
    378                    scaled = yDot[0].clone();
    379                    for (int j = 0; j < scaled.length; ++j) {
    380                        scaled[j] *= stepSize;
    381                    }
    382
    383                    // higher order derivatives
    384                    nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
    385
    386                    // stop the integrator now that all needed steps have been handled
    387                    throw new InitializationCompletedMarkerException();
    388
    389                }
    390
    391            }
    392
    393            /** {@inheritDoc} */
    394            public void init(double t0, double[] y0, double time) {
    395                // nothing to do
    396            }
    397
    398        }
    399
    400        /** Marker exception used ONLY to stop the starter integrator after first step. */
    401        private static class InitializationCompletedMarkerException
    402            extends RuntimeException {
    403
    404            /** Serializable version identifier. */
    405            private static final long serialVersionUID = -1914085471038046418L;
    406
    407            /** Simple constructor. */
    408            public InitializationCompletedMarkerException() {
    409                super((Throwable) null);
    410            }
    411
    412        }
    413
    414        /** Wrapper for differential equations, ensuring start evaluations are counted. */
    415        private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {
    416
    417            /** Dimension of the problem. */
    418            private final int dimension;
    419
    420            /** Simple constructor.
    421             * @param dimension dimension of the problem
    422             */
    423            public CountingDifferentialEquations(final int dimension) {
    424                this.dimension = dimension;
    425            }
    426
    427            /** {@inheritDoc} */
    428            public void computeDerivatives(double t, double[] y, double[] dot)
    429                throws MaxCountExceededException, DimensionMismatchException {
    430                MultistepIntegrator.this.computeDerivatives(t, y, dot);
    431            }
    432
    433            /** {@inheritDoc} */
    434            public int getDimension() {
    435                return dimension;
    436            }
    437
    438        }
    439
    440    }