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.math3.ode;
19  
20  import org.apache.commons.math3.exception.DimensionMismatchException;
21  import org.apache.commons.math3.exception.MaxCountExceededException;
22  import org.apache.commons.math3.exception.NoBracketingException;
23  import org.apache.commons.math3.exception.NumberIsTooSmallException;
24  import org.apache.commons.math3.exception.util.LocalizedFormats;
25  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
26  import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
27  import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
28  import org.apache.commons.math3.ode.sampling.StepHandler;
29  import org.apache.commons.math3.ode.sampling.StepInterpolator;
30  import org.apache.commons.math3.util.FastMath;
31  
32  /**
33   * This class is the base class for multistep integrators for Ordinary
34   * Differential Equations.
35   * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
36   * <pre>
37   * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
38   * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
39   * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
40   * ...
41   * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
42   * </pre></p>
43   * <p>Rather than storing several previous steps separately, this implementation uses
44   * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
45   * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
46   * <pre>
47   * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
48   * </pre>
49   * (we omit the k index in the notation for clarity)</p>
50   * <p>
51   * Multistep integrators with Nordsieck representation are highly sensitive to
52   * large step changes because when the step is multiplied by factor a, the
53   * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
54   * and the last components are the least accurate ones. The default max growth
55   * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
56   * </p>
57   *
58   * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator
59   * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator
60   * @since 2.0
61   */
62  public 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<sup>(k)</sup>)</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      * @exception NumberIsTooSmallException if number of steps is smaller than 2
109      */
110     protected MultistepIntegrator(final String name, final int nSteps,
111                                   final int order,
112                                   final double minStep, final double maxStep,
113                                   final double scalAbsoluteTolerance,
114                                   final double scalRelativeTolerance)
115         throws NumberIsTooSmallException {
116 
117         super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
118 
119         if (nSteps < 2) {
120             throw new NumberIsTooSmallException(
121                   LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
122                   nSteps, 2, true);
123         }
124 
125         starter = new DormandPrince853Integrator(minStep, maxStep,
126                                                  scalAbsoluteTolerance,
127                                                  scalRelativeTolerance);
128         this.nSteps = nSteps;
129 
130         exp = -1.0 / order;
131 
132         // set the default values of the algorithm control parameters
133         setSafety(0.9);
134         setMinReduction(0.2);
135         setMaxGrowth(FastMath.pow(2.0, -exp));
136 
137     }
138 
139     /**
140      * Build a multistep integrator with the given stepsize bounds.
141      * <p>The default starter integrator is set to the {@link
142      * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
143      * some defaults settings.</p>
144      * <p>
145      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
146      * </p>
147      * @param name name of the method
148      * @param nSteps number of steps of the multistep method
149      * (excluding the one being computed)
150      * @param order order of the method
151      * @param minStep minimal step (must be positive even for backward
152      * integration), the last step can be smaller than this
153      * @param maxStep maximal step (must be positive even for backward
154      * integration)
155      * @param vecAbsoluteTolerance allowed absolute error
156      * @param vecRelativeTolerance allowed relative error
157      */
158     protected MultistepIntegrator(final String name, final int nSteps,
159                                   final int order,
160                                   final double minStep, final double maxStep,
161                                   final double[] vecAbsoluteTolerance,
162                                   final double[] vecRelativeTolerance) {
163         super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
164         starter = new DormandPrince853Integrator(minStep, maxStep,
165                                                  vecAbsoluteTolerance,
166                                                  vecRelativeTolerance);
167         this.nSteps = nSteps;
168 
169         exp = -1.0 / order;
170 
171         // set the default values of the algorithm control parameters
172         setSafety(0.9);
173         setMinReduction(0.2);
174         setMaxGrowth(FastMath.pow(2.0, -exp));
175 
176     }
177 
178     /**
179      * Get the starter integrator.
180      * @return starter integrator
181      */
182     public ODEIntegrator getStarterIntegrator() {
183         return starter;
184     }
185 
186     /**
187      * Set the starter integrator.
188      * <p>The various step and event handlers for this starter integrator
189      * will be managed automatically by the multi-step integrator. Any
190      * user configuration for these elements will be cleared before use.</p>
191      * @param starterIntegrator starter integrator
192      */
193     public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
194         this.starter = starterIntegrator;
195     }
196 
197     /** Start the integration.
198      * <p>This method computes one step using the underlying starter integrator,
199      * and initializes the Nordsieck vector at step start. The starter integrator
200      * purpose is only to establish initial conditions, it does not really change
201      * time by itself. The top level multistep integrator remains in charge of
202      * handling time propagation and events handling as it will starts its own
203      * computation right from the beginning. In a sense, the starter integrator
204      * can be seen as a dummy one and so it will never trigger any user event nor
205      * call any user step handler.</p>
206      * @param t0 initial time
207      * @param y0 initial value of the state vector at t0
208      * @param t target time for the integration
209      * (can be set to a value smaller than <code>t0</code> for backward integration)
210      * @exception DimensionMismatchException if arrays dimension do not match equations settings
211      * @exception NumberIsTooSmallException if integration step is too small
212      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
213      * @exception NoBracketingException if the location of an event cannot be bracketed
214      */
215     protected void start(final double t0, final double[] y0, final double t)
216         throws DimensionMismatchException, NumberIsTooSmallException,
217                MaxCountExceededException, NoBracketingException {
218 
219         // make sure NO user event nor user step handler is triggered,
220         // this is the task of the top level integrator, not the task
221         // of the starter integrator
222         starter.clearEventHandlers();
223         starter.clearStepHandlers();
224 
225         // set up one specific step handler to extract initial Nordsieck vector
226         starter.addStepHandler(new NordsieckInitializer(nSteps, y0.length));
227 
228         // start integration, expecting a InitializationCompletedMarkerException
229         try {
230 
231             if (starter instanceof AbstractIntegrator) {
232                 ((AbstractIntegrator) starter).integrate(getExpandable(), t);
233             } else {
234                 starter.integrate(new FirstOrderDifferentialEquations() {
235 
236                     /** {@inheritDoc} */
237                     public int getDimension() {
238                         return getExpandable().getTotalDimension();
239                     }
240 
241                     /** {@inheritDoc} */
242                     public void computeDerivatives(double t, double[] y, double[] yDot) {
243                         getExpandable().computeDerivatives(t, y, yDot);
244                     }
245 
246                 }, t0, y0, t, new double[y0.length]);
247             }
248 
249         } catch (InitializationCompletedMarkerException icme) { // NOPMD
250             // this is the expected nominal interruption of the start integrator
251 
252             // count the evaluations used by the starter
253             getEvaluationsCounter().incrementCount(starter.getEvaluations());
254 
255         }
256 
257         // remove the specific step handler
258         starter.clearStepHandlers();
259 
260     }
261 
262     /** Initialize the high order scaled derivatives at step start.
263      * @param h step size to use for scaling
264      * @param t first steps times
265      * @param y first steps states
266      * @param yDot first steps derivatives
267      * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
268      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
269      */
270     protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
271                                                                            final double[][] y,
272                                                                            final double[][] yDot);
273 
274     /** Get the minimal reduction factor for stepsize control.
275      * @return minimal reduction factor
276      */
277     public double getMinReduction() {
278         return minReduction;
279     }
280 
281     /** Set the minimal reduction factor for stepsize control.
282      * @param minReduction minimal reduction factor
283      */
284     public void setMinReduction(final double minReduction) {
285         this.minReduction = minReduction;
286     }
287 
288     /** Get the maximal growth factor for stepsize control.
289      * @return maximal growth factor
290      */
291     public double getMaxGrowth() {
292         return maxGrowth;
293     }
294 
295     /** Set the maximal growth factor for stepsize control.
296      * @param maxGrowth maximal growth factor
297      */
298     public void setMaxGrowth(final double maxGrowth) {
299         this.maxGrowth = maxGrowth;
300     }
301 
302     /** Get the safety factor for stepsize control.
303      * @return safety factor
304      */
305     public double getSafety() {
306       return safety;
307     }
308 
309     /** Set the safety factor for stepsize control.
310      * @param safety safety factor
311      */
312     public void setSafety(final double safety) {
313       this.safety = safety;
314     }
315 
316     /** Compute step grow/shrink factor according to normalized error.
317      * @param error normalized error of the current step
318      * @return grow/shrink factor for next step
319      */
320     protected double computeStepGrowShrinkFactor(final double error) {
321         return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
322     }
323 
324     /** Transformer used to convert the first step to Nordsieck representation. */
325     public interface NordsieckTransformer {
326         /** Initialize the high order scaled derivatives at step start.
327          * @param h step size to use for scaling
328          * @param t first steps times
329          * @param y first steps states
330          * @param yDot first steps derivatives
331          * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
332          * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
333          */
334         Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
335                                                             final double[][] y,
336                                                             final double[][] yDot);
337     }
338 
339     /** Specialized step handler storing the first step. */
340     private class NordsieckInitializer implements StepHandler {
341 
342         /** Steps counter. */
343         private int count;
344 
345         /** First steps times. */
346         private final double[] t;
347 
348         /** First steps states. */
349         private final double[][] y;
350 
351         /** First steps derivatives. */
352         private final double[][] yDot;
353 
354         /** Simple constructor.
355          * @param nSteps number of steps of the multistep method (excluding the one being computed)
356          * @param n problem dimension
357          */
358         public NordsieckInitializer(final int nSteps, final int n) {
359             this.count = 0;
360             this.t     = new double[nSteps];
361             this.y     = new double[nSteps][n];
362             this.yDot  = new double[nSteps][n];
363         }
364 
365         /** {@inheritDoc} */
366         public void handleStep(StepInterpolator interpolator, boolean isLast)
367             throws MaxCountExceededException {
368 
369             final double prev = interpolator.getPreviousTime();
370             final double curr = interpolator.getCurrentTime();
371 
372             if (count == 0) {
373                 // first step, we need to store also the beginning of the step
374                 interpolator.setInterpolatedTime(prev);
375                 t[0] = prev;
376                 final ExpandableStatefulODE expandable = getExpandable();
377                 final EquationsMapper primary = expandable.getPrimaryMapper();
378                 primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
379                 primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
380                 int index = 0;
381                 for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
382                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
383                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
384                     ++index;
385                 }
386             }
387 
388             // store the end of the step
389             ++count;
390             interpolator.setInterpolatedTime(curr);
391             t[count] = curr;
392 
393             final ExpandableStatefulODE expandable = getExpandable();
394             final EquationsMapper primary = expandable.getPrimaryMapper();
395             primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
396             primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
397             int index = 0;
398             for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
399                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
400                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
401                 ++index;
402             }
403 
404             if (count == t.length - 1) {
405 
406                 // this was the last step we needed, we can compute the derivatives
407                 stepStart = t[0];
408                 stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
409 
410                 // first scaled derivative
411                 scaled = yDot[0].clone();
412                 for (int j = 0; j < scaled.length; ++j) {
413                     scaled[j] *= stepSize;
414                 }
415 
416                 // higher order derivatives
417                 nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
418 
419                 // stop the integrator now that all needed steps have been handled
420                 throw new InitializationCompletedMarkerException();
421 
422             }
423 
424         }
425 
426         /** {@inheritDoc} */
427         public void init(double t0, double[] y0, double time) {
428             // nothing to do
429         }
430 
431     }
432 
433     /** Marker exception used ONLY to stop the starter integrator after first step. */
434     private static class InitializationCompletedMarkerException
435         extends RuntimeException {
436 
437         /** Serializable version identifier. */
438         private static final long serialVersionUID = -1914085471038046418L;
439 
440         /** Simple constructor. */
441         public InitializationCompletedMarkerException() {
442             super((Throwable) null);
443         }
444 
445     }
446 
447 }