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 1463684 2013-04-02 19:04:13Z luc $
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    
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            } catch (InitializationCompletedMarkerException icme) { // NOPMD
251                // this is the expected nominal interruption of the start integrator
252    
253                // count the evaluations used by the starter
254                getEvaluationsCounter().incrementCount(starter.getEvaluations());
255    
256            }
257    
258            // remove the specific step handler
259            starter.clearStepHandlers();
260    
261        }
262    
263        /** Initialize the high order scaled derivatives at step start.
264         * @param h step size to use for scaling
265         * @param t first steps times
266         * @param y first steps states
267         * @param yDot first steps derivatives
268         * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
269         * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
270         */
271        protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
272                                                                               final double[][] y,
273                                                                               final double[][] yDot);
274    
275        /** Get the minimal reduction factor for stepsize control.
276         * @return minimal reduction factor
277         */
278        public double getMinReduction() {
279            return minReduction;
280        }
281    
282        /** Set the minimal reduction factor for stepsize control.
283         * @param minReduction minimal reduction factor
284         */
285        public void setMinReduction(final double minReduction) {
286            this.minReduction = minReduction;
287        }
288    
289        /** Get the maximal growth factor for stepsize control.
290         * @return maximal growth factor
291         */
292        public double getMaxGrowth() {
293            return maxGrowth;
294        }
295    
296        /** Set the maximal growth factor for stepsize control.
297         * @param maxGrowth maximal growth factor
298         */
299        public void setMaxGrowth(final double maxGrowth) {
300            this.maxGrowth = maxGrowth;
301        }
302    
303        /** Get the safety factor for stepsize control.
304         * @return safety factor
305         */
306        public double getSafety() {
307          return safety;
308        }
309    
310        /** Set the safety factor for stepsize control.
311         * @param safety safety factor
312         */
313        public void setSafety(final double safety) {
314          this.safety = safety;
315        }
316    
317        /** Compute step grow/shrink factor according to normalized error.
318         * @param error normalized error of the current step
319         * @return grow/shrink factor for next step
320         */
321        protected double computeStepGrowShrinkFactor(final double error) {
322            return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
323        }
324    
325        /** Transformer used to convert the first step to Nordsieck representation. */
326        public interface NordsieckTransformer {
327            /** Initialize the high order scaled derivatives at step start.
328             * @param h step size to use for scaling
329             * @param t first steps times
330             * @param y first steps states
331             * @param yDot first steps derivatives
332             * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
333             * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
334             */
335            Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
336                                                                final double[][] y,
337                                                                final double[][] yDot);
338        }
339    
340        /** Specialized step handler storing the first step. */
341        private class NordsieckInitializer implements StepHandler {
342    
343            /** Steps counter. */
344            private int count;
345    
346            /** First steps times. */
347            private final double[] t;
348    
349            /** First steps states. */
350            private final double[][] y;
351    
352            /** First steps derivatives. */
353            private final double[][] yDot;
354    
355            /** Simple constructor.
356             * @param nSteps number of steps of the multistep method (excluding the one being computed)
357             * @param n problem dimension
358             */
359            public NordsieckInitializer(final int nSteps, final int n) {
360                this.count = 0;
361                this.t     = new double[nSteps];
362                this.y     = new double[nSteps][n];
363                this.yDot  = new double[nSteps][n];
364            }
365    
366            /** {@inheritDoc} */
367            public void handleStep(StepInterpolator interpolator, boolean isLast)
368                throws MaxCountExceededException {
369    
370                final double prev = interpolator.getPreviousTime();
371                final double curr = interpolator.getCurrentTime();
372    
373                if (count == 0) {
374                    // first step, we need to store also the beginning of the step
375                    interpolator.setInterpolatedTime(prev);
376                    t[0] = prev;
377                    final ExpandableStatefulODE expandable = getExpandable();
378                    final EquationsMapper primary = expandable.getPrimaryMapper();
379                    primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
380                    primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
381                    int index = 0;
382                    for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
383                        secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
384                        secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
385                        ++index;
386                    }
387                }
388    
389                // store the end of the step
390                ++count;
391                interpolator.setInterpolatedTime(curr);
392                t[count] = curr;
393    
394                final ExpandableStatefulODE expandable = getExpandable();
395                final EquationsMapper primary = expandable.getPrimaryMapper();
396                primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
397                primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
398                int index = 0;
399                for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
400                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
401                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
402                    ++index;
403                }
404    
405                if (count == t.length - 1) {
406    
407                    // this was the last step we needed, we can compute the derivatives
408                    stepStart = t[0];
409                    stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
410    
411                    // first scaled derivative
412                    scaled = yDot[0].clone();
413                    for (int j = 0; j < scaled.length; ++j) {
414                        scaled[j] *= stepSize;
415                    }
416    
417                    // higher order derivatives
418                    nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
419    
420                    // stop the integrator now that all needed steps have been handled
421                    throw new InitializationCompletedMarkerException();
422    
423                }
424    
425            }
426    
427            /** {@inheritDoc} */
428            public void init(double t0, double[] y0, double time) {
429                // nothing to do
430            }
431    
432        }
433    
434        /** Marker exception used ONLY to stop the starter integrator after first step. */
435        private static class InitializationCompletedMarkerException
436            extends RuntimeException {
437    
438            /** Serializable version identifier. */
439            private static final long serialVersionUID = -1914085471038046418L;
440    
441            /** Simple constructor. */
442            public InitializationCompletedMarkerException() {
443                super((Throwable) null);
444            }
445    
446        }
447    
448    }