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.MathIllegalStateException;
22  import org.apache.commons.math3.exception.MaxCountExceededException;
23  import org.apache.commons.math3.exception.NoBracketingException;
24  import org.apache.commons.math3.exception.NumberIsTooSmallException;
25  import org.apache.commons.math3.exception.util.LocalizedFormats;
26  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
27  import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator;
28  import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
29  import org.apache.commons.math3.ode.sampling.StepHandler;
30  import org.apache.commons.math3.ode.sampling.StepInterpolator;
31  import org.apache.commons.math3.util.FastMath;
32  
33  /**
34   * This class is the base class for multistep integrators for Ordinary
35   * Differential Equations.
36   * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
37   * <pre>
38   * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
39   * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
40   * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
41   * ...
42   * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
43   * </pre></p>
44   * <p>Rather than storing several previous steps separately, this implementation uses
45   * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
46   * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
47   * <pre>
48   * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
49   * </pre>
50   * (we omit the k index in the notation for clarity)</p>
51   * <p>
52   * Multistep integrators with Nordsieck representation are highly sensitive to
53   * large step changes because when the step is multiplied by factor a, the
54   * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
55   * and the last components are the least accurate ones. The default max growth
56   * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
57   * </p>
58   *
59   * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator
60   * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator
61   * @since 2.0
62   */
63  public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
64  
65      /** First scaled derivative (h y'). */
66      protected double[] scaled;
67  
68      /** Nordsieck matrix of the higher scaled derivatives.
69       * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
70       */
71      protected Array2DRowRealMatrix nordsieck;
72  
73      /** Starter integrator. */
74      private FirstOrderIntegrator starter;
75  
76      /** Number of steps of the multistep method (excluding the one being computed). */
77      private final int nSteps;
78  
79      /** Stepsize control exponent. */
80      private double exp;
81  
82      /** Safety factor for stepsize control. */
83      private double safety;
84  
85      /** Minimal reduction factor for stepsize control. */
86      private double minReduction;
87  
88      /** Maximal growth factor for stepsize control. */
89      private double maxGrowth;
90  
91      /**
92       * Build a multistep integrator with the given stepsize bounds.
93       * <p>The default starter integrator is set to the {@link
94       * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
95       * some defaults settings.</p>
96       * <p>
97       * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
98       * </p>
99       * @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 + 3) / 2, 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             // we should not reach this step
251             throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
252 
253         } catch (InitializationCompletedMarkerException icme) { // NOPMD
254             // this is the expected nominal interruption of the start integrator
255 
256             // count the evaluations used by the starter
257             getCounter().increment(starter.getEvaluations());
258 
259         }
260 
261         // remove the specific step handler
262         starter.clearStepHandlers();
263 
264     }
265 
266     /** Initialize the high order scaled derivatives at step start.
267      * @param h step size to use for scaling
268      * @param t first steps times
269      * @param y first steps states
270      * @param yDot first steps derivatives
271      * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
272      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
273      */
274     protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
275                                                                            final double[][] y,
276                                                                            final double[][] yDot);
277 
278     /** Get the minimal reduction factor for stepsize control.
279      * @return minimal reduction factor
280      */
281     public double getMinReduction() {
282         return minReduction;
283     }
284 
285     /** Set the minimal reduction factor for stepsize control.
286      * @param minReduction minimal reduction factor
287      */
288     public void setMinReduction(final double minReduction) {
289         this.minReduction = minReduction;
290     }
291 
292     /** Get the maximal growth factor for stepsize control.
293      * @return maximal growth factor
294      */
295     public double getMaxGrowth() {
296         return maxGrowth;
297     }
298 
299     /** Set the maximal growth factor for stepsize control.
300      * @param maxGrowth maximal growth factor
301      */
302     public void setMaxGrowth(final double maxGrowth) {
303         this.maxGrowth = maxGrowth;
304     }
305 
306     /** Get the safety factor for stepsize control.
307      * @return safety factor
308      */
309     public double getSafety() {
310       return safety;
311     }
312 
313     /** Set the safety factor for stepsize control.
314      * @param safety safety factor
315      */
316     public void setSafety(final double safety) {
317       this.safety = safety;
318     }
319 
320     /** Get the number of steps of the multistep method (excluding the one being computed).
321      * @return number of steps of the multistep method (excluding the one being computed)
322      */
323     public int getNSteps() {
324       return nSteps;
325     }
326 
327     /** Compute step grow/shrink factor according to normalized error.
328      * @param error normalized error of the current step
329      * @return grow/shrink factor for next step
330      */
331     protected double computeStepGrowShrinkFactor(final double error) {
332         return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
333     }
334 
335     /** Transformer used to convert the first step to Nordsieck representation.
336      * @deprecated as of 3.6 this unused interface is deprecated
337      */
338     @Deprecated
339     public interface NordsieckTransformer {
340         /** Initialize the high order scaled derivatives at step start.
341          * @param h step size to use for scaling
342          * @param t first steps times
343          * @param y first steps states
344          * @param yDot first steps derivatives
345          * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
346          * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
347          */
348         Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
349                                                             final double[][] y,
350                                                             final double[][] yDot);
351     }
352 
353     /** Specialized step handler storing the first step. */
354     private class NordsieckInitializer implements StepHandler {
355 
356         /** Steps counter. */
357         private int count;
358 
359         /** First steps times. */
360         private final double[] t;
361 
362         /** First steps states. */
363         private final double[][] y;
364 
365         /** First steps derivatives. */
366         private final double[][] yDot;
367 
368         /** Simple constructor.
369          * @param nbStartPoints number of start points (including the initial point)
370          * @param n problem dimension
371          */
372         NordsieckInitializer(final int nbStartPoints, final int n) {
373             this.count = 0;
374             this.t     = new double[nbStartPoints];
375             this.y     = new double[nbStartPoints][n];
376             this.yDot  = new double[nbStartPoints][n];
377         }
378 
379         /** {@inheritDoc} */
380         public void handleStep(StepInterpolator interpolator, boolean isLast)
381             throws MaxCountExceededException {
382 
383             final double prev = interpolator.getPreviousTime();
384             final double curr = interpolator.getCurrentTime();
385 
386             if (count == 0) {
387                 // first step, we need to store also the point at the beginning of the step
388                 interpolator.setInterpolatedTime(prev);
389                 t[0] = prev;
390                 final ExpandableStatefulODE expandable = getExpandable();
391                 final EquationsMapper primary = expandable.getPrimaryMapper();
392                 primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
393                 primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
394                 int index = 0;
395                 for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
396                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
397                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
398                     ++index;
399                 }
400             }
401 
402             // store the point at the end of the step
403             ++count;
404             interpolator.setInterpolatedTime(curr);
405             t[count] = curr;
406 
407             final ExpandableStatefulODE expandable = getExpandable();
408             final EquationsMapper primary = expandable.getPrimaryMapper();
409             primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
410             primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
411             int index = 0;
412             for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
413                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
414                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
415                 ++index;
416             }
417 
418             if (count == t.length - 1) {
419 
420                 // this was the last point we needed, we can compute the derivatives
421                 stepStart = t[0];
422                 stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
423 
424                 // first scaled derivative
425                 scaled = yDot[0].clone();
426                 for (int j = 0; j < scaled.length; ++j) {
427                     scaled[j] *= stepSize;
428                 }
429 
430                 // higher order derivatives
431                 nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
432 
433                 // stop the integrator now that all needed steps have been handled
434                 throw new InitializationCompletedMarkerException();
435 
436             }
437 
438         }
439 
440         /** {@inheritDoc} */
441         public void init(double t0, double[] y0, double time) {
442             // nothing to do
443         }
444 
445     }
446 
447     /** Marker exception used ONLY to stop the starter integrator after first step. */
448     private static class InitializationCompletedMarkerException
449         extends RuntimeException {
450 
451         /** Serializable version identifier. */
452         private static final long serialVersionUID = -1914085471038046418L;
453 
454         /** Simple constructor. */
455         InitializationCompletedMarkerException() {
456             super((Throwable) null);
457         }
458 
459     }
460 
461 }