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.ode.DerivativeException;
022    import org.apache.commons.math.exception.util.LocalizedFormats;
023    import org.apache.commons.math.linear.Array2DRowRealMatrix;
024    import org.apache.commons.math.linear.RealMatrix;
025    import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
026    import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
027    import org.apache.commons.math.ode.sampling.StepHandler;
028    import org.apache.commons.math.ode.sampling.StepInterpolator;
029    import org.apache.commons.math.util.FastMath;
030    
031    /**
032     * This class is the base class for multistep integrators for Ordinary
033     * Differential Equations.
034     * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
035     * <pre>
036     * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
037     * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
038     * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
039     * ...
040     * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
041     * </pre></p>
042     * <p>Rather than storing several previous steps separately, this implementation uses
043     * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
044     * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
045     * <pre>
046     * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
047     * </pre>
048     * (we omit the k index in the notation for clarity)</p>
049     * <p>
050     * Multistep integrators with Nordsieck representation are highly sensitive to
051     * large step changes because when the step is multiplied by a factor a, the
052     * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
053     * and the last components are the least accurate ones. The default max growth
054     * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
055     * </p>
056     *
057     * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
058     * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
059     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
060     * @since 2.0
061     */
062    public 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(k))</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         */
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    }