1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package org.apache.commons.math.ode;
19
20import org.apache.commons.math.MathRuntimeException;
21import org.apache.commons.math.ode.DerivativeException;
22import org.apache.commons.math.exception.util.LocalizedFormats;
23import org.apache.commons.math.linear.Array2DRowRealMatrix;
24import org.apache.commons.math.linear.RealMatrix;
25import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
26import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
27import org.apache.commons.math.ode.sampling.StepHandler;
28import org.apache.commons.math.ode.sampling.StepInterpolator;
29import org.apache.commons.math.util.FastMath;
30
31/**
32 * This class is the base class for multistep integrators for Ordinary
33 * Differential Equations.
34 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
35 * <pre>
36 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
37 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
38 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
39 * ...
40 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
41 * </pre></p>
42 * <p>Rather than storing several previous steps separately, this implementation uses
43 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
44 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
45 * <pre>
46 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
47 * </pre>
48 * (we omit the k index in the notation for clarity)</p>
49 * <p>
50 * Multistep integrators with Nordsieck representation are highly sensitive to
51 * large step changes because when the step is multiplied by a factor a, the
52 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
53 * and the last components are the least accurate ones. The default max growth
54 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
55 * </p>
56 *
57 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
58 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
59 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
60 * @since 2.0
61 */
62public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
63
64    /** First scaled derivative (h y'). */
65    protected double[] scaled;
66
67    /** Nordsieck matrix of the higher scaled derivatives.
68     * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p>
69     */
70    protected Array2DRowRealMatrix nordsieck;
71
72    /** Starter integrator. */
73    private FirstOrderIntegrator starter;
74
75    /** Number of steps of the multistep method (excluding the one being computed). */
76    private final int nSteps;
77
78    /** Stepsize control exponent. */
79    private double exp;
80
81    /** Safety factor for stepsize control. */
82    private double safety;
83
84    /** Minimal reduction factor for stepsize control. */
85    private double minReduction;
86
87    /** Maximal growth factor for stepsize control. */
88    private double maxGrowth;
89
90    /**
91     * Build a multistep integrator with the given stepsize bounds.
92     * <p>The default starter integrator is set to the {@link
93     * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
94     * some defaults settings.</p>
95     * <p>
96     * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
97     * </p>
98     * @param name name of the method
99     * @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}
413