View Javadoc

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  
18  package org.apache.commons.math.ode;
19  
20  import org.apache.commons.math.MathRuntimeException;
21  import org.apache.commons.math.linear.Array2DRowRealMatrix;
22  import org.apache.commons.math.linear.RealMatrix;
23  import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
24  import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
25  import org.apache.commons.math.ode.sampling.StepHandler;
26  import org.apache.commons.math.ode.sampling.StepInterpolator;
27  
28  /**
29   * This class is the base class for multistep integrators for Ordinary
30   * Differential Equations.
31   * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
32   * <pre>
33   * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
34   * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
35   * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
36   * ...
37   * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
38   * </pre></p>
39   * <p>Rather than storing several previous steps separately, this implementation uses
40   * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
41   * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
42   * <pre>
43   * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
44   * </pre>
45   * (we omit the k index in the notation for clarity)</p>
46   * <p>
47   * Multistep integrators with Nordsieck representation are highly sensitive to
48   * large step changes because when the step is multiplied by a factor a, the
49   * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
50   * and the last components are the least accurate ones. The default max growth
51   * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
52   * </p>
53   *
54   * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
55   * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
56   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
57   * @since 2.0
58   */
59  public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
60  
61      /** Starter integrator. */
62      private FirstOrderIntegrator starter;
63  
64      /** Number of steps of the multistep method (excluding the one being computed). */
65      private final int nSteps;
66  
67      /** First scaled derivative (h y'). */
68      protected double[] scaled;
69  
70      /** Nordsieck matrix of the higher scaled derivatives.
71       * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p>
72       */
73      protected Array2DRowRealMatrix nordsieck;
74  
75      /** Stepsize control exponent. */
76      private double exp;
77  
78      /** Safety factor for stepsize control. */
79      private double safety;
80  
81      /** Minimal reduction factor for stepsize control. */
82      private double minReduction;
83  
84      /** Maximal growth factor for stepsize control. */
85      private double maxGrowth;
86  
87      /**
88       * Build a multistep integrator with the given stepsize bounds.
89       * <p>The default starter integrator is set to the {@link
90       * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
91       * some defaults settings.</p>
92       * <p>
93       * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
94       * </p>
95       * @param name name of the method
96       * @param nSteps number of steps of the multistep method
97       * (excluding the one being computed)
98       * @param order order of the method
99       * @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 starter starter integrator
187      */
188     public void setStarterIntegrator(FirstOrderIntegrator starter) {
189         this.starter = starter;
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 }