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.MaxCountExceededException;
022import org.apache.commons.math3.exception.NoBracketingException;
023import org.apache.commons.math3.exception.NumberIsTooSmallException;
024import org.apache.commons.math3.exception.util.LocalizedFormats;
025import org.apache.commons.math3.linear.Array2DRowRealMatrix;
026import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
027import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
028import org.apache.commons.math3.ode.sampling.StepHandler;
029import org.apache.commons.math3.ode.sampling.StepInterpolator;
030import 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 * @since 2.0
061 */
062public 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<sup>(k)</sup>)</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     * @exception NumberIsTooSmallException if number of steps is smaller than 2
109     */
110    protected MultistepIntegrator(final String name, final int nSteps,
111                                  final int order,
112                                  final double minStep, final double maxStep,
113                                  final double scalAbsoluteTolerance,
114                                  final double scalRelativeTolerance)
115        throws NumberIsTooSmallException {
116
117        super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
118
119        if (nSteps < 2) {
120            throw new NumberIsTooSmallException(
121                  LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
122                  nSteps, 2, true);
123        }
124
125        starter = new DormandPrince853Integrator(minStep, maxStep,
126                                                 scalAbsoluteTolerance,
127                                                 scalRelativeTolerance);
128        this.nSteps = nSteps;
129
130        exp = -1.0 / order;
131
132        // set the default values of the algorithm control parameters
133        setSafety(0.9);
134        setMinReduction(0.2);
135        setMaxGrowth(FastMath.pow(2.0, -exp));
136
137    }
138
139    /**
140     * Build a multistep integrator with the given stepsize bounds.
141     * <p>The default starter integrator is set to the {@link
142     * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
143     * some defaults settings.</p>
144     * <p>
145     * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
146     * </p>
147     * @param name name of the method
148     * @param nSteps number of steps of the multistep method
149     * (excluding the one being computed)
150     * @param order order of the method
151     * @param minStep minimal step (must be positive even for backward
152     * integration), the last step can be smaller than this
153     * @param maxStep maximal step (must be positive even for backward
154     * integration)
155     * @param vecAbsoluteTolerance allowed absolute error
156     * @param vecRelativeTolerance allowed relative error
157     */
158    protected MultistepIntegrator(final String name, final int nSteps,
159                                  final int order,
160                                  final double minStep, final double maxStep,
161                                  final double[] vecAbsoluteTolerance,
162                                  final double[] vecRelativeTolerance) {
163        super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
164        starter = new DormandPrince853Integrator(minStep, maxStep,
165                                                 vecAbsoluteTolerance,
166                                                 vecRelativeTolerance);
167        this.nSteps = nSteps;
168
169        exp = -1.0 / order;
170
171        // set the default values of the algorithm control parameters
172        setSafety(0.9);
173        setMinReduction(0.2);
174        setMaxGrowth(FastMath.pow(2.0, -exp));
175
176    }
177
178    /**
179     * Get the starter integrator.
180     * @return starter integrator
181     */
182    public ODEIntegrator getStarterIntegrator() {
183        return starter;
184    }
185
186    /**
187     * Set the starter integrator.
188     * <p>The various step and event handlers for this starter integrator
189     * will be managed automatically by the multi-step integrator. Any
190     * user configuration for these elements will be cleared before use.</p>
191     * @param starterIntegrator starter integrator
192     */
193    public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
194        this.starter = starterIntegrator;
195    }
196
197    /** Start the integration.
198     * <p>This method computes one step using the underlying starter integrator,
199     * and initializes the Nordsieck vector at step start. The starter integrator
200     * purpose is only to establish initial conditions, it does not really change
201     * time by itself. The top level multistep integrator remains in charge of
202     * handling time propagation and events handling as it will starts its own
203     * computation right from the beginning. In a sense, the starter integrator
204     * can be seen as a dummy one and so it will never trigger any user event nor
205     * call any user step handler.</p>
206     * @param t0 initial time
207     * @param y0 initial value of the state vector at t0
208     * @param t target time for the integration
209     * (can be set to a value smaller than <code>t0</code> for backward integration)
210     * @exception DimensionMismatchException if arrays dimension do not match equations settings
211     * @exception NumberIsTooSmallException if integration step is too small
212     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
213     * @exception NoBracketingException if the location of an event cannot be bracketed
214     */
215    protected void start(final double t0, final double[] y0, final double t)
216        throws DimensionMismatchException, NumberIsTooSmallException,
217               MaxCountExceededException, NoBracketingException {
218
219        // make sure NO user event nor user step handler is triggered,
220        // this is the task of the top level integrator, not the task
221        // of the starter integrator
222        starter.clearEventHandlers();
223        starter.clearStepHandlers();
224
225        // set up one specific step handler to extract initial Nordsieck vector
226        starter.addStepHandler(new NordsieckInitializer(nSteps, y0.length));
227
228        // start integration, expecting a InitializationCompletedMarkerException
229        try {
230
231            if (starter instanceof AbstractIntegrator) {
232                ((AbstractIntegrator) starter).integrate(getExpandable(), t);
233            } else {
234                starter.integrate(new FirstOrderDifferentialEquations() {
235
236                    /** {@inheritDoc} */
237                    public int getDimension() {
238                        return getExpandable().getTotalDimension();
239                    }
240
241                    /** {@inheritDoc} */
242                    public void computeDerivatives(double t, double[] y, double[] yDot) {
243                        getExpandable().computeDerivatives(t, y, yDot);
244                    }
245
246                }, t0, y0, t, new double[y0.length]);
247            }
248
249        } catch (InitializationCompletedMarkerException icme) { // NOPMD
250            // this is the expected nominal interruption of the start integrator
251
252            // count the evaluations used by the starter
253            getEvaluationsCounter().incrementCount(starter.getEvaluations());
254
255        }
256
257        // remove the specific step handler
258        starter.clearStepHandlers();
259
260    }
261
262    /** Initialize the high order scaled derivatives at step start.
263     * @param h step size to use for scaling
264     * @param t first steps times
265     * @param y first steps states
266     * @param yDot first steps derivatives
267     * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
268     * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
269     */
270    protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
271                                                                           final double[][] y,
272                                                                           final double[][] yDot);
273
274    /** Get the minimal reduction factor for stepsize control.
275     * @return minimal reduction factor
276     */
277    public double getMinReduction() {
278        return minReduction;
279    }
280
281    /** Set the minimal reduction factor for stepsize control.
282     * @param minReduction minimal reduction factor
283     */
284    public void setMinReduction(final double minReduction) {
285        this.minReduction = minReduction;
286    }
287
288    /** Get the maximal growth factor for stepsize control.
289     * @return maximal growth factor
290     */
291    public double getMaxGrowth() {
292        return maxGrowth;
293    }
294
295    /** Set the maximal growth factor for stepsize control.
296     * @param maxGrowth maximal growth factor
297     */
298    public void setMaxGrowth(final double maxGrowth) {
299        this.maxGrowth = maxGrowth;
300    }
301
302    /** Get the safety factor for stepsize control.
303     * @return safety factor
304     */
305    public double getSafety() {
306      return safety;
307    }
308
309    /** Set the safety factor for stepsize control.
310     * @param safety safety factor
311     */
312    public void setSafety(final double safety) {
313      this.safety = safety;
314    }
315
316    /** Compute step grow/shrink factor according to normalized error.
317     * @param error normalized error of the current step
318     * @return grow/shrink factor for next step
319     */
320    protected double computeStepGrowShrinkFactor(final double error) {
321        return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
322    }
323
324    /** Transformer used to convert the first step to Nordsieck representation. */
325    public interface NordsieckTransformer {
326        /** Initialize the high order scaled derivatives at step start.
327         * @param h step size to use for scaling
328         * @param t first steps times
329         * @param y first steps states
330         * @param yDot first steps derivatives
331         * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
332         * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
333         */
334        Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
335                                                            final double[][] y,
336                                                            final double[][] yDot);
337    }
338
339    /** Specialized step handler storing the first step. */
340    private class NordsieckInitializer implements StepHandler {
341
342        /** Steps counter. */
343        private int count;
344
345        /** First steps times. */
346        private final double[] t;
347
348        /** First steps states. */
349        private final double[][] y;
350
351        /** First steps derivatives. */
352        private final double[][] yDot;
353
354        /** Simple constructor.
355         * @param nSteps number of steps of the multistep method (excluding the one being computed)
356         * @param n problem dimension
357         */
358        public NordsieckInitializer(final int nSteps, final int n) {
359            this.count = 0;
360            this.t     = new double[nSteps];
361            this.y     = new double[nSteps][n];
362            this.yDot  = new double[nSteps][n];
363        }
364
365        /** {@inheritDoc} */
366        public void handleStep(StepInterpolator interpolator, boolean isLast)
367            throws MaxCountExceededException {
368
369            final double prev = interpolator.getPreviousTime();
370            final double curr = interpolator.getCurrentTime();
371
372            if (count == 0) {
373                // first step, we need to store also the beginning of the step
374                interpolator.setInterpolatedTime(prev);
375                t[0] = prev;
376                final ExpandableStatefulODE expandable = getExpandable();
377                final EquationsMapper primary = expandable.getPrimaryMapper();
378                primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
379                primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
380                int index = 0;
381                for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
382                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
383                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
384                    ++index;
385                }
386            }
387
388            // store the end of the step
389            ++count;
390            interpolator.setInterpolatedTime(curr);
391            t[count] = curr;
392
393            final ExpandableStatefulODE expandable = getExpandable();
394            final EquationsMapper primary = expandable.getPrimaryMapper();
395            primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
396            primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
397            int index = 0;
398            for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
399                secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
400                secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
401                ++index;
402            }
403
404            if (count == t.length - 1) {
405
406                // this was the last step we needed, we can compute the derivatives
407                stepStart = t[0];
408                stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
409
410                // first scaled derivative
411                scaled = yDot[0].clone();
412                for (int j = 0; j < scaled.length; ++j) {
413                    scaled[j] *= stepSize;
414                }
415
416                // higher order derivatives
417                nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
418
419                // stop the integrator now that all needed steps have been handled
420                throw new InitializationCompletedMarkerException();
421
422            }
423
424        }
425
426        /** {@inheritDoc} */
427        public void init(double t0, double[] y0, double time) {
428            // nothing to do
429        }
430
431    }
432
433    /** Marker exception used ONLY to stop the starter integrator after first step. */
434    private static class InitializationCompletedMarkerException
435        extends RuntimeException {
436
437        /** Serializable version identifier. */
438        private static final long serialVersionUID = -1914085471038046418L;
439
440        /** Simple constructor. */
441        public InitializationCompletedMarkerException() {
442            super((Throwable) null);
443        }
444
445    }
446
447}