    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.math.ode;
    019
    020    import org.apache.commons.math.MathRuntimeException;
    021    import org.apache.commons.math.ode.DerivativeException;
    022    import org.apache.commons.math.exception.util.LocalizedFormats;
    023    import org.apache.commons.math.linear.Array2DRowRealMatrix;
    024    import org.apache.commons.math.linear.RealMatrix;
    025    import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
    026    import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
    027    import org.apache.commons.math.ode.sampling.StepHandler;
    028    import org.apache.commons.math.ode.sampling.StepInterpolator;
    029    import org.apache.commons.math.util.FastMath;
    030
    031    /**
    032     * This class is the base class for multistep integrators for Ordinary
    033     * Differential Equations.
    034     * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
    035     * <pre>
    036     * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
    037     * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
    038     * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
    039     * ...
    040     * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
    041     * </pre></p>
    042     * <p>Rather than storing several previous steps separately, this implementation uses
    043     * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
    044     * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
    045     * <pre>
    046     * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
    047     * </pre>
    048     * (we omit the k index in the notation for clarity)</p>
    049     * <p>
    050     * Multistep integrators with Nordsieck representation are highly sensitive to
    051     * large step changes because when the step is multiplied by a factor a, the
    052     * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
    053     * and the last components are the least accurate ones. The default max growth
    054     * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
    055     * </p>
    056     *
    057     * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
    058     * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
    059     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    060     * @since 2.0
    061     */
    062    public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
    063
    064        /** First scaled derivative (h y'). */
    065        protected double[] scaled;
    066
    067        /** Nordsieck matrix of the higher scaled derivatives.
    068         * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p>
    069         */
    070        protected Array2DRowRealMatrix nordsieck;
    071
    072        /** Starter integrator. */
    073        private FirstOrderIntegrator starter;
    074
    075        /** Number of steps of the multistep method (excluding the one being computed). */
    076        private final int nSteps;
    077
    078        /** Stepsize control exponent. */
    079        private double exp;
    080
    081        /** Safety factor for stepsize control. */
    082        private double safety;
    083
    084        /** Minimal reduction factor for stepsize control. */
    085        private double minReduction;
    086
    087        /** Maximal growth factor for stepsize control. */
    088        private double maxGrowth;
    089
    090        /**
    091         * Build a multistep integrator with the given stepsize bounds.
    092         * <p>The default starter integrator is set to the {@link
    093         * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
    094         * some defaults settings.</p>
    095         * <p>
    096         * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
    097         * </p>
    098         * @param name name of the method
    099         * @param nSteps number of steps of the multistep method
    100         * (excluding the one being computed)
    101         * @param order order of the method
    102         * @param minStep minimal step (must be positive even for backward
    103         * integration), the last step can be smaller than this
    104         * @param maxStep maximal step (must be positive even for backward
    105         * integration)
    106         * @param scalAbsoluteTolerance allowed absolute error
    107         * @param scalRelativeTolerance allowed relative error
    108         */
    109        protected MultistepIntegrator(final String name, final int nSteps,
    110                                      final int order,
    111                                      final double minStep, final double maxStep,
    112                                      final double scalAbsoluteTolerance,
    113                                      final double scalRelativeTolerance) {
    114
    115            super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
    116
    117            if (nSteps <= 0) {
    118                throw MathRuntimeException.createIllegalArgumentException(
    119                      LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT,
    120                      name);
    121            }
    122
    123            starter = new DormandPrince853Integrator(minStep, maxStep,
    124                                                     scalAbsoluteTolerance,
    125                                                     scalRelativeTolerance);
    126            this.nSteps = nSteps;
    127
    128            exp = -1.0 / order;
    129
    130            // set the default values of the algorithm control parameters
    131            setSafety(0.9);
    132            setMinReduction(0.2);
    133            setMaxGrowth(FastMath.pow(2.0, -exp));
    134
    135        }
    136
    137        /**
    138         * Build a multistep integrator with the given stepsize bounds.
    139         * <p>The default starter integrator is set to the {@link
    140         * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
    141         * some defaults settings.</p>
    142         * <p>
    143         * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
    144         * </p>
    145         * @param name name of the method
    146         * @param nSteps number of steps of the multistep method
    147         * (excluding the one being computed)
    148         * @param order order of the method
    149         * @param minStep minimal step (must be positive even for backward
    150         * integration), the last step can be smaller than this
    151         * @param maxStep maximal step (must be positive even for backward
    152         * integration)
    153         * @param vecAbsoluteTolerance allowed absolute error
    154         * @param vecRelativeTolerance allowed relative error
    155         */
    156        protected MultistepIntegrator(final String name, final int nSteps,
    157                                      final int order,
    158                                      final double minStep, final double maxStep,
    159                                      final double[] vecAbsoluteTolerance,
    160                                      final double[] vecRelativeTolerance) {
    161            super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
    162            starter = new DormandPrince853Integrator(minStep, maxStep,
    163                                                     vecAbsoluteTolerance,
    164                                                     vecRelativeTolerance);
    165            this.nSteps = nSteps;
    166
    167            exp = -1.0 / order;
    168
    169            // set the default values of the algorithm control parameters
    170            setSafety(0.9);
    171            setMinReduction(0.2);
    172            setMaxGrowth(FastMath.pow(2.0, -exp));
    173
    174        }
    175
    176        /**
    177         * Get the starter integrator.
    178         * @return starter integrator
    179         */
    180        public ODEIntegrator getStarterIntegrator() {
    181            return starter;
    182        }
    183
    184        /**
    185         * Set the starter integrator.
    186         * <p>The various step and event handlers for this starter integrator
    187         * will be managed automatically by the multi-step integrator. Any
    188         * user configuration for these elements will be cleared before use.</p>
    189         * @param starterIntegrator starter integrator
    190         */
    191        public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
    192            this.starter = starterIntegrator;
    193        }
    194
    195        /** Start the integration.
    196         * <p>This method computes one step using the underlying starter integrator,
    197         * and initializes the Nordsieck vector at step start. The starter integrator
    198         * purpose is only to establish initial conditions, it does not really change
    199         * time by itself. The top level multistep integrator remains in charge of
    200         * handling time propagation and events handling as it will starts its own
    201         * computation right from the beginning. In a sense, the starter integrator
    202         * can be seen as a dummy one and so it will never trigger any user event nor
    203         * call any user step handler.</p>
    204         * @param t0 initial time
    205         * @param y0 initial value of the state vector at t0
    206         * @param t target time for the integration
    207         * (can be set to a value smaller than <code>t0</code> for backward integration)
    208         * @throws IntegratorException if the integrator cannot perform integration
    209         * @throws DerivativeException this exception is propagated to the caller if
    210         * the underlying user function triggers one
    211         */
    212        protected void start(final double t0, final double[] y0, final double t)
    213            throws DerivativeException, IntegratorException {
    214
    215            // make sure NO user event nor user step handler is triggered,
    216            // this is the task of the top level integrator, not the task
    217            // of the starter integrator
    218            starter.clearEventHandlers();
    219            starter.clearStepHandlers();
    220
    221            // set up one specific step handler to extract initial Nordsieck vector
    222            starter.addStepHandler(new NordsieckInitializer(y0.length));
    223
    224            // start integration, expecting a InitializationCompletedMarkerException
    225            try {
    226                starter.integrate(new CountingDifferentialEquations(y0.length),
    227                                  t0, y0, t, new double[y0.length]);
    228            } catch (DerivativeException mue) {
    229                if (!(mue instanceof InitializationCompletedMarkerException)) {
    230                    // this is not the expected nominal interruption of the start integrator
    231                    throw mue;
    232                }
    233            }
    234
    235            // remove the specific step handler
    236            starter.clearStepHandlers();
    237
    238        }
    239
    240        /** Initialize the high order scaled derivatives at step start.
    241         * @param first first scaled derivative at step start
    242         * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
    243         * will be modified
    244         * @return high order scaled derivatives at step start
    245         */
    246        protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
    247                                                                               final double[][] multistep);
    248
    249        /** Get the minimal reduction factor for stepsize control.
    250         * @return minimal reduction factor
    251         */
    252        public double getMinReduction() {
    253            return minReduction;
    254        }
    255
    256        /** Set the minimal reduction factor for stepsize control.
    257         * @param minReduction minimal reduction factor
    258         */
    259        public void setMinReduction(final double minReduction) {
    260            this.minReduction = minReduction;
    261        }
    262
    263        /** Get the maximal growth factor for stepsize control.
    264         * @return maximal growth factor
    265         */
    266        public double getMaxGrowth() {
    267            return maxGrowth;
    268        }
    269
    270        /** Set the maximal growth factor for stepsize control.
    271         * @param maxGrowth maximal growth factor
    272         */
    273        public void setMaxGrowth(final double maxGrowth) {
    274            this.maxGrowth = maxGrowth;
    275        }
    276
    277        /** Get the safety factor for stepsize control.
    278         * @return safety factor
    279         */
    280        public double getSafety() {
    281          return safety;
    282        }
    283
    284        /** Set the safety factor for stepsize control.
    285         * @param safety safety factor
    286         */
    287        public void setSafety(final double safety) {
    288          this.safety = safety;
    289        }
    290
    291        /** Compute step grow/shrink factor according to normalized error.
    292         * @param error normalized error of the current step
    293         * @return grow/shrink factor for next step
    294         */
    295        protected double computeStepGrowShrinkFactor(final double error) {
    296            return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
    297        }
    298
    299        /** Transformer used to convert the first step to Nordsieck representation. */
    300        public static interface NordsieckTransformer {
    301            /** Initialize the high order scaled derivatives at step start.
    302             * @param first first scaled derivative at step start
    303             * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
    304             * will be modified
    305             * @return high order derivatives at step start
    306             */
    307            RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
    308        }
    309
    310        /** Specialized step handler storing the first step. */
    311        private class NordsieckInitializer implements StepHandler {
    312
    313            /** Problem dimension. */
    314            private final int n;
    315
    316            /** Simple constructor.
    317             * @param n problem dimension
    318             */
    319            public NordsieckInitializer(final int n) {
    320                this.n = n;
    321            }
    322
    323            /** {@inheritDoc} */
    324            public void handleStep(StepInterpolator interpolator, boolean isLast)
    325                throws DerivativeException {
    326
    327                final double prev = interpolator.getPreviousTime();
    328                final double curr = interpolator.getCurrentTime();
    329                stepStart = prev;
    330                stepSize  = (curr - prev) / (nSteps + 1);
    331
    332                // compute the first scaled derivative
    333                interpolator.setInterpolatedTime(prev);
    334                scaled = interpolator.getInterpolatedDerivatives().clone();
    335                for (int j = 0; j < n; ++j) {
    336                    scaled[j] *= stepSize;
    337                }
    338
    339                // compute the high order scaled derivatives
    340                final double[][] multistep = new double[nSteps][];
    341                for (int i = 1; i <= nSteps; ++i) {
    342                    interpolator.setInterpolatedTime(prev + stepSize * i);
    343                    final double[] msI = interpolator.getInterpolatedDerivatives().clone();
    344                    for (int j = 0; j < n; ++j) {
    345                        msI[j] *= stepSize;
    346                    }
    347                    multistep[i - 1] = msI;
    348                }
    349                nordsieck = initializeHighOrderDerivatives(scaled, multistep);
    350
    351                // stop the integrator after the first step has been handled
    352                throw new InitializationCompletedMarkerException();
    353
    354            }
    355
    356            /** {@inheritDoc} */
    357            public boolean requiresDenseOutput() {
    358                return true;
    359            }
    360
    361            /** {@inheritDoc} */
    362            public void reset() {
    363                // nothing to do
    364            }
    365
    366        }
    367
    368        /** Marker exception used ONLY to stop the starter integrator after first step. */
    369        private static class InitializationCompletedMarkerException
    370            extends DerivativeException {
    371
    372            /** Serializable version identifier. */
    373            private static final long serialVersionUID = -4105805787353488365L;
    374
    375            /** Simple constructor. */
    376            public InitializationCompletedMarkerException() {
    377                super((Throwable) null);
    378            }
    379
    380        }
    381
    382        /** Wrapper for differential equations, ensuring start evaluations are counted. */
    383        private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations {
    384
    385            /** Dimension of the problem. */
    386            private final int dimension;
    387
    388            /** Simple constructor.
    389             * @param dimension dimension of the problem
    390             */
    391            public CountingDifferentialEquations(final int dimension) {
    392                this.dimension = dimension;
    393            }
    394
    395            /** {@inheritDoc} */
    396            public void computeDerivatives(double t, double[] y, double[] dot)
    397                    throws DerivativeException {
    398                MultistepIntegrator.this.computeDerivatives(t, y, dot);
    399            }
    400
    401            /** {@inheritDoc} */
    402            public int getDimension() {
    403                return dimension;
    404            }
    405
    406            /** {@inheritDoc} */
    407            public int getMainSetDimension() {
    408                return mainSetDimension;
    409            }
    410        }
    411
    412    }