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
018package org.apache.commons.math3.ode;
019
020import org.apache.commons.math3.exception.DimensionMismatchException;
021import org.apache.commons.math3.exception.MathIllegalStateException;
022import org.apache.commons.math3.exception.MaxCountExceededException;
023import org.apache.commons.math3.exception.NoBracketingException;
024import org.apache.commons.math3.exception.NumberIsTooSmallException;
025import org.apache.commons.math3.exception.util.LocalizedFormats;
026import org.apache.commons.math3.linear.Array2DRowRealMatrix;
027import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
028import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
029import org.apache.commons.math3.ode.sampling.StepHandler;
030import org.apache.commons.math3.ode.sampling.StepInterpolator;
031import org.apache.commons.math3.util.FastMath;
032
033/**
034 * This class is the base class for multistep integrators for Ordinary
035 * Differential Equations.
036 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
037 * <pre>
038 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
039 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
040 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
041 * ...
042 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
043 * </pre></p>
044 * <p>Rather than storing several previous steps separately, this implementation uses
045 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
046 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
047 * <pre>
048 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
049 * </pre>
050 * (we omit the k index in the notation for clarity)</p>
051 * <p>
052 * Multistep integrators with Nordsieck representation are highly sensitive to
053 * large step changes because when the step is multiplied by factor a, the
054 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
055 * and the last components are the least accurate ones. The default max growth
056 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
057 * </p>
058 *
059 * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator
060 * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator
061 * @since 2.0
062 */
063public 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 + 3) / 2, y0.length));
228
229        // start integration, expecting a InitializationCompletedMarkerException
230        try {
231
232            if (starter instanceof AbstractIntegrator) {
233                ((AbstractIntegrator) starter).integrate(getExpandable(), t);
234            } else {
235                starter.integrate(new FirstOrderDifferentialEquations() {
236
237                    /** {@inheritDoc} */
238                    public int getDimension() {
239                        return getExpandable().getTotalDimension();
240                    }
241
242                    /** {@inheritDoc} */
243                    public void computeDerivatives(double t, double[] y, double[] yDot) {
244                        getExpandable().computeDerivatives(t, y, yDot);
245                    }
246
247                }, t0, y0, t, new double[y0.length]);
248            }
249
250            // we should not reach this step
251            throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
252
253        } catch (InitializationCompletedMarkerException icme) { // NOPMD
254            // this is the expected nominal interruption of the start integrator
255
256            // count the evaluations used by the starter
257            getCounter().increment(starter.getEvaluations());
258
259        }
260
261        // remove the specific step handler
262        starter.clearStepHandlers();
263
264    }
265
266    /** Initialize the high order scaled derivatives at step start.
267     * @param h step size to use for scaling
268     * @param t first steps times
269     * @param y first steps states
270     * @param yDot first steps derivatives
271     * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
272     * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
273     */
274    protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
275                                                                           final double[][] y,
276                                                                           final double[][] yDot);
277
278    /** Get the minimal reduction factor for stepsize control.
279     * @return minimal reduction factor
280     */
281    public double getMinReduction() {
282        return minReduction;
283    }
284
285    /** Set the minimal reduction factor for stepsize control.
286     * @param minReduction minimal reduction factor
287     */
288    public void setMinReduction(final double minReduction) {
289        this.minReduction = minReduction;
290    }
291
292    /** Get the maximal growth factor for stepsize control.
293     * @return maximal growth factor
294     */
295    public double getMaxGrowth() {
296        return maxGrowth;
297    }
298
299    /** Set the maximal growth factor for stepsize control.
300     * @param maxGrowth maximal growth factor
301     */
302    public void setMaxGrowth(final double maxGrowth) {
303        this.maxGrowth = maxGrowth;
304    }
305
306    /** Get the safety factor for stepsize control.
307     * @return safety factor
308     */
309    public double getSafety() {
310      return safety;
311    }
312
313    /** Set the safety factor for stepsize control.
314     * @param safety safety factor
315     */
316    public void setSafety(final double safety) {
317      this.safety = safety;
318    }
319
320    /** Get the number of steps of the multistep method (excluding the one being computed).
321     * @return number of steps of the multistep method (excluding the one being computed)
322     */
323    public int getNSteps() {
324      return nSteps;
325    }
326
327    /** Compute step grow/shrink factor according to normalized error.
328     * @param error normalized error of the current step
329     * @return grow/shrink factor for next step
330     */
331    protected double computeStepGrowShrinkFactor(final double error) {
332        return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
333    }
334
335    /** Transformer used to convert the first step to Nordsieck representation.
336     * @deprecated as of 3.6 this unused interface is deprecated
337     */
338    @Deprecated
339    public interface NordsieckTransformer {
340        /** Initialize the high order scaled derivatives at step start.
341         * @param h step size to use for scaling
342         * @param t first steps times
343         * @param y first steps states
344         * @param yDot first steps derivatives
345         * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
346         * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
347         */
348        Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
349                                                            final double[][] y,
350                                                            final double[][] yDot);
351    }
352
353    /** Specialized step handler storing the first step. */
354    private class NordsieckInitializer implements StepHandler {
355
356        /** Steps counter. */
357        private int count;
358
359        /** First steps times. */
360        private final double[] t;
361
362        /** First steps states. */
363        private final double[][] y;
364
365        /** First steps derivatives. */
366        private final double[][] yDot;
367
368        /** Simple constructor.
369         * @param nbStartPoints number of start points (including the initial point)
370         * @param n problem dimension
371         */
372        NordsieckInitializer(final int nbStartPoints, final int n) {
373            this.count = 0;
374            this.t     = new double[nbStartPoints];
375            this.y     = new double[nbStartPoints][n];
376            this.yDot  = new double[nbStartPoints][n];
377        }
378
379        /** {@inheritDoc} */
380        public void handleStep(StepInterpolator interpolator, boolean isLast)
381            throws MaxCountExceededException {
382
383            final double prev = interpolator.getPreviousTime();
384            final double curr = interpolator.getCurrentTime();
385
386            if (count == 0) {
387                // first step, we need to store also the point at the beginning of the step
388                interpolator.setInterpolatedTime(prev);
389                t[0] = prev;
390                final ExpandableStatefulODE expandable = getExpandable();
391                final EquationsMapper primary = expandable.getPrimaryMapper();
392                primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
393                primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
394                int index = 0;
395                for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
396                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
397                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
398                    ++index;
399                }
400            }
401
402            // store the point at the end of the step
403            ++count;
404            interpolator.setInterpolatedTime(curr);
405            t[count] = curr;
406
407            final ExpandableStatefulODE expandable = getExpandable();
408            final EquationsMapper primary = expandable.getPrimaryMapper();
409            primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
410            primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
411            int index = 0;
412            for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
413                secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
414                secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
415                ++index;
416            }
417
418            if (count == t.length - 1) {
419
420                // this was the last point we needed, we can compute the derivatives
421                stepStart = t[0];
422                stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
423
424                // first scaled derivative
425                scaled = yDot[0].clone();
426                for (int j = 0; j < scaled.length; ++j) {
427                    scaled[j] *= stepSize;
428                }
429
430                // higher order derivatives
431                nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
432
433                // stop the integrator now that all needed steps have been handled
434                throw new InitializationCompletedMarkerException();
435
436            }
437
438        }
439
440        /** {@inheritDoc} */
441        public void init(double t0, double[] y0, double time) {
442            // nothing to do
443        }
444
445    }
446
447    /** Marker exception used ONLY to stop the starter integrator after first step. */
448    private static class InitializationCompletedMarkerException
449        extends RuntimeException {
450
451        /** Serializable version identifier. */
452        private static final long serialVersionUID = -1914085471038046418L;
453
454        /** Simple constructor. */
455        InitializationCompletedMarkerException() {
456            super((Throwable) null);
457        }
458
459    }
460
461}