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