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    *      https://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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.ode;
24  
25  import org.hipparchus.Field;
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.linear.Array2DRowFieldMatrix;
30  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
31  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
32  import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
33  import org.hipparchus.ode.sampling.FieldODEStepHandler;
34  import org.hipparchus.util.FastMath;
35  import org.hipparchus.util.MathArrays;
36  
37  /**
38   * This class is the base class for multistep integrators for Ordinary
39   * Differential Equations.
40   * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
41   * \[
42   *   \left\{\begin{align}
43   *   s_1(n) &amp;= h y'_n \text{ for first derivative}\\
44   *   s_2(n) &amp;= \frac{h^2}{2} y_n'' \text{ for second derivative}\\
45   *   s_3(n) &amp;= \frac{h^3}{6} y_n''' \text{ for third derivative}\\
46   *   &amp;\cdots\\
47   *   s_k(n) &amp;= \frac{h^k}{k!} y_n^{(k)} \text{ for } k^\mathrm{th} \text{ derivative}
48   *   \end{align}\right.
49   * \]</p>
50   * <p>Rather than storing several previous steps separately, this implementation uses
51   * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
52   * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
53   * \[
54   * r_n = [ s_2(n), s_3(n) \ldots s_k(n) ]^T
55   * \]
56   * (we omit the k index in the notation for clarity)</p>
57   * <p>
58   * Multistep integrators with Nordsieck representation are highly sensitive to
59   * large step changes because when the step is multiplied by factor a, the
60   * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
61   * and the last components are the least accurate ones. The default max growth
62   * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
63   * </p>
64   *
65   * @see org.hipparchus.ode.nonstiff.AdamsBashforthFieldIntegrator
66   * @see org.hipparchus.ode.nonstiff.AdamsMoultonFieldIntegrator
67   * @param <T> the type of the field elements
68   */
69  public abstract class MultistepFieldIntegrator<T extends CalculusFieldElement<T>>
70      extends AdaptiveStepsizeFieldIntegrator<T> {
71  
72      /** First scaled derivative (h y'). */
73      protected T[] scaled;
74  
75      /** Nordsieck matrix of the higher scaled derivatives.
76       * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
77       */
78      protected Array2DRowFieldMatrix<T> nordsieck;
79  
80      /** Starter integrator. */
81      private FieldODEIntegrator<T> starter;
82  
83      /** Number of steps of the multistep method (excluding the one being computed). */
84      private final int nSteps;
85  
86      /** Stepsize control exponent. */
87      private double exp;
88  
89      /** Safety factor for stepsize control. */
90      private double safety;
91  
92      /** Minimal reduction factor for stepsize control. */
93      private double minReduction;
94  
95      /** Maximal growth factor for stepsize control. */
96      private double maxGrowth;
97  
98      /**
99       * Build a multistep integrator with the given stepsize bounds.
100      * <p>The default starter integrator is set to the {@link
101      * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
102      * some defaults settings.</p>
103      * <p>
104      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
105      * </p>
106      * @param field field to which the time and state vector elements belong
107      * @param name name of the method
108      * @param nSteps number of steps of the multistep method
109      * (excluding the one being computed)
110      * @param order order of the method
111      * @param minStep minimal step (must be positive even for backward
112      * integration), the last step can be smaller than this
113      * @param maxStep maximal step (must be positive even for backward
114      * integration)
115      * @param scalAbsoluteTolerance allowed absolute error
116      * @param scalRelativeTolerance allowed relative error
117      * @exception MathIllegalArgumentException if number of steps is smaller than 2
118      */
119     protected MultistepFieldIntegrator(final Field<T> field, final String name,
120                                        final int nSteps, final int order,
121                                        final double minStep, final double maxStep,
122                                        final double scalAbsoluteTolerance,
123                                        final double scalRelativeTolerance)
124         throws MathIllegalArgumentException {
125 
126         super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
127 
128         if (nSteps < 2) {
129             throw new MathIllegalArgumentException(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
130                                                    nSteps, 2, true);
131         }
132 
133         starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
134                                                         scalAbsoluteTolerance,
135                                                         scalRelativeTolerance);
136         this.nSteps = nSteps;
137 
138         exp = -1.0 / order;
139 
140         // set the default values of the algorithm control parameters
141         setSafety(0.9);
142         setMinReduction(0.2);
143         setMaxGrowth(FastMath.pow(2.0, -exp));
144 
145     }
146 
147     /**
148      * Build a multistep integrator with the given stepsize bounds.
149      * <p>The default starter integrator is set to the {@link
150      * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
151      * some defaults settings.</p>
152      * <p>
153      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
154      * </p>
155      * @param field field to which the time and state vector elements belong
156      * @param name name of the method
157      * @param nSteps number of steps of the multistep method
158      * (excluding the one being computed)
159      * @param order order of the method
160      * @param minStep minimal step (must be positive even for backward
161      * integration), the last step can be smaller than this
162      * @param maxStep maximal step (must be positive even for backward
163      * integration)
164      * @param vecAbsoluteTolerance allowed absolute error
165      * @param vecRelativeTolerance allowed relative error
166      */
167     protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
168                                        final int order,
169                                        final double minStep, final double maxStep,
170                                        final double[] vecAbsoluteTolerance,
171                                        final double[] vecRelativeTolerance) {
172         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
173 
174         if (nSteps < 2) {
175             throw new MathIllegalArgumentException(LocalizedODEFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
176                                                    nSteps, 2, true);
177         }
178 
179         starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
180                                                         vecAbsoluteTolerance,
181                                                         vecRelativeTolerance);
182         this.nSteps = nSteps;
183 
184         exp = -1.0 / order;
185 
186         // set the default values of the algorithm control parameters
187         setSafety(0.9);
188         setMinReduction(0.2);
189         setMaxGrowth(FastMath.pow(2.0, -exp));
190 
191     }
192 
193     /**
194      * Get the starter integrator.
195      * @return starter integrator
196      */
197     public FieldODEIntegrator<T> getStarterIntegrator() {
198         return starter;
199     }
200 
201     /**
202      * Set the starter integrator.
203      * <p>The various step and event handlers for this starter integrator
204      * will be managed automatically by the multi-step integrator. Any
205      * user configuration for these elements will be cleared before use.</p>
206      * @param starterIntegrator starter integrator
207      */
208     public void setStarterIntegrator(FieldODEIntegrator<T> starterIntegrator) {
209         this.starter = starterIntegrator;
210     }
211 
212     /** Start the integration.
213      * <p>This method computes one step using the underlying starter integrator,
214      * and initializes the Nordsieck vector at step start. The starter integrator
215      * purpose is only to establish initial conditions, it does not really change
216      * time by itself. The top level multistep integrator remains in charge of
217      * handling time propagation and events handling as it will starts its own
218      * computation right from the beginning. In a sense, the starter integrator
219      * can be seen as a dummy one and so it will never trigger any user event nor
220      * call any user step handler.</p>
221      * @param equations complete set of differential equations to integrate
222      * @param initialState initial state (time, primary and secondary state vectors)
223      * @param t target time for the integration
224      * (can be set to a value smaller than <code>t0</code> for backward integration)
225      * @exception MathIllegalArgumentException if arrays dimension do not match equations settings
226      * @exception MathIllegalArgumentException if integration step is too small
227      * @exception MathIllegalStateException if the number of functions evaluations is exceeded
228      * @exception MathIllegalArgumentException if the location of an event cannot be bracketed
229      */
230     protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t)
231         throws MathIllegalArgumentException, MathIllegalStateException {
232 
233         // make sure NO user events nor user step handlers are triggered,
234         // this is the task of the top level integrator, not the task of the starter integrator
235         starter.clearEventDetectors();
236         starter.clearStepHandlers();
237 
238         // set up one specific step handler to extract initial Nordsieck vector
239         starter.addStepHandler(new FieldNordsieckInitializer((nSteps + 3) / 2));
240 
241         // start integration, expecting a InitializationCompletedMarkerException
242         try {
243 
244             starter.integrate(equations, initialState, t);
245 
246             // we should not reach this step
247             throw new MathIllegalStateException(LocalizedODEFormats.MULTISTEP_STARTER_STOPPED_EARLY);
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().increment(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 Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(T h, T[] t, T[][] y, T[][] yDot);
271 
272     /** Get the minimal reduction factor for stepsize control.
273      * @return minimal reduction factor
274      */
275     public double getMinReduction() {
276         return minReduction;
277     }
278 
279     /** Set the minimal reduction factor for stepsize control.
280      * @param minReduction minimal reduction factor
281      */
282     public void setMinReduction(final double minReduction) {
283         this.minReduction = minReduction;
284     }
285 
286     /** Get the maximal growth factor for stepsize control.
287      * @return maximal growth factor
288      */
289     public double getMaxGrowth() {
290         return maxGrowth;
291     }
292 
293     /** Set the maximal growth factor for stepsize control.
294      * @param maxGrowth maximal growth factor
295      */
296     public void setMaxGrowth(final double maxGrowth) {
297         this.maxGrowth = maxGrowth;
298     }
299 
300     /** Get the safety factor for stepsize control.
301      * @return safety factor
302      */
303     public double getSafety() {
304       return safety;
305     }
306 
307     /** Set the safety factor for stepsize control.
308      * @param safety safety factor
309      */
310     public void setSafety(final double safety) {
311       this.safety = safety;
312     }
313 
314     /** Get the number of steps of the multistep method (excluding the one being computed).
315      * @return number of steps of the multistep method (excluding the one being computed)
316      */
317     public int getNSteps() {
318       return nSteps;
319     }
320 
321     /** Rescale the instance.
322      * <p>Since the scaled and Nordsieck arrays are shared with the caller,
323      * this method has the side effect of rescaling this arrays in the caller too.</p>
324      * @param newStepSize new step size to use in the scaled and Nordsieck arrays
325      */
326     protected void rescale(final T newStepSize) {
327 
328         final T ratio = newStepSize.divide(getStepSize());
329         for (int i = 0; i < scaled.length; ++i) {
330             scaled[i] = scaled[i].multiply(ratio);
331         }
332 
333         final T[][] nData = nordsieck.getDataRef();
334         T power = ratio;
335         for (int i = 0; i < nData.length; ++i) {
336             power = power.multiply(ratio);
337             final T[] nDataI = nData[i];
338             for (int j = 0; j < nDataI.length; ++j) {
339                 nDataI[j] = nDataI[j].multiply(power);
340             }
341         }
342 
343         setStepSize(newStepSize);
344 
345     }
346 
347     /** Compute step grow/shrink factor according to normalized error.
348      * @param error normalized error of the current step
349      * @return grow/shrink factor for next step
350      */
351     protected double computeStepGrowShrinkFactor(final double error) {
352         return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
353     }
354 
355     /** Specialized step handler storing the first step.
356      */
357     private class FieldNordsieckInitializer implements FieldODEStepHandler<T> {
358 
359         /** Steps counter. */
360         private int count;
361 
362         /** Saved start. */
363         private FieldODEStateAndDerivative<T> savedStart;
364 
365         /** First steps times. */
366         private final T[] t;
367 
368         /** First steps states. */
369         private final T[][] y;
370 
371         /** First steps derivatives. */
372         private final T[][] yDot;
373 
374         /** Simple constructor.
375          * @param nbStartPoints number of start points (including the initial point)
376          */
377         FieldNordsieckInitializer(final int nbStartPoints) {
378             this.count  = 0;
379             this.t      = MathArrays.buildArray(getField(), nbStartPoints);
380             this.y      = MathArrays.buildArray(getField(), nbStartPoints, -1);
381             this.yDot   = MathArrays.buildArray(getField(), nbStartPoints, -1);
382         }
383 
384         /** {@inheritDoc} */
385         @Override
386         public void handleStep(FieldODEStateInterpolator<T> interpolator) {
387 
388 
389             if (count == 0) {
390                 // first step, we need to store also the point at the beginning of the step
391                 final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
392                 savedStart  = prev;
393                 t[count]    = prev.getTime();
394                 y[count]    = prev.getCompleteState();
395                 yDot[count] = prev.getCompleteDerivative();
396             }
397 
398             // store the point at the end of the step
399             ++count;
400             final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
401             t[count]    = curr.getTime();
402             y[count]    = curr.getCompleteState();
403             yDot[count] = curr.getCompleteDerivative();
404 
405             if (count == t.length - 1) {
406 
407                 // this was the last point we needed, we can compute the derivatives
408                 setStepStart(savedStart);
409                 final T rawStep = t[t.length - 1].subtract(t[0]).divide(t.length - 1);
410                 setStepSize(getStepSizeHelper().filterStep(rawStep, rawStep.getReal() >= 0, true));
411 
412                 // first scaled derivative
413                 scaled = MathArrays.buildArray(getField(), yDot[0].length);
414                 for (int j = 0; j < scaled.length; ++j) {
415                     scaled[j] = yDot[0][j].multiply(getStepSize());
416                 }
417 
418                 // higher order derivatives
419                 nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);
420 
421                 // stop the integrator now that all needed steps have been handled
422                 throw new InitializationCompletedMarkerException();
423 
424             }
425 
426         }
427 
428         /** {@inheritDoc} */
429         @Override
430         public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
431             // nothing to do
432         }
433 
434     }
435 
436     /** Marker exception used ONLY to stop the starter integrator after first step. */
437     private static class InitializationCompletedMarkerException
438         extends RuntimeException {
439 
440         /** Serializable version identifier. */
441         private static final long serialVersionUID = -1914085471038046418L;
442 
443         /** Simple constructor. */
444         InitializationCompletedMarkerException() {
445             super((Throwable) null);
446         }
447 
448     }
449 
450 }