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.nonstiff;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.ode.FieldEquationsMapper;
30  import org.hipparchus.ode.FieldExpandableODE;
31  import org.hipparchus.ode.FieldODEState;
32  import org.hipparchus.ode.FieldODEStateAndDerivative;
33  import org.hipparchus.util.FastMath;
34  import org.hipparchus.util.MathArrays;
35  import org.hipparchus.util.MathUtils;
36  
37  /**
38   * This class implements the common part of all embedded Runge-Kutta
39   * integrators for Ordinary Differential Equations.
40   *
41   * <p>These methods are embedded explicit Runge-Kutta methods with two
42   * sets of coefficients allowing to estimate the error, their Butcher
43   * arrays are as follows :</p>
44   * <pre>
45   *    0  |
46   *   c2  | a21
47   *   c3  | a31  a32
48   *   ... |        ...
49   *   cs  | as1  as2  ...  ass-1
50   *       |--------------------------
51   *       |  b1   b2  ...   bs-1  bs
52   *       |  b'1  b'2 ...   b's-1 b's
53   * </pre>
54   *
55   * <p>In fact, we rather use the array defined by ej = bj - b'j to
56   * compute directly the error rather than computing two estimates and
57   * then comparing them.</p>
58   *
59   * <p>Some methods are qualified as <i>fsal</i> (first same as last)
60   * methods. This means the last evaluation of the derivatives in one
61   * step is the same as the first in the next step. Then, this
62   * evaluation can be reused from one step to the next one and the cost
63   * of such a method is really s-1 evaluations despite the method still
64   * has s stages. This behaviour is true only for successful steps, if
65   * the step is rejected after the error estimation phase, no
66   * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
67   * asi = bi for all i.</p>
68   *
69   * @param <T> the type of the field elements
70   */
71  
72  public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
73      extends AdaptiveStepsizeFieldIntegrator<T>
74      implements FieldExplicitRungeKuttaIntegrator<T> {
75  
76      /** Index of the pre-computed derivative for <i>fsal</i> methods. */
77      private final int fsal;
78  
79      /** Time steps from Butcher array (without the first zero). */
80      private final T[] c;
81  
82      /** Internal weights from Butcher array (without the first empty row). */
83      private final T[][] a;
84  
85      /** External weights for the high order method from Butcher array. */
86      private final T[] b;
87  
88      /** Time steps from Butcher array (without the first zero). */
89      private double[] realC = new double[0];
90  
91      /** Internal weights from Butcher array (without the first empty row). Real version, optional. */
92      private double[][] realA = new double[0][];
93  
94      /** External weights for the high order method from Butcher array. Real version, optional. */
95      private double[] realB = new double[0];
96  
97      /** Stepsize control exponent. */
98      private final double exp;
99  
100     /** Safety factor for stepsize control. */
101     private T safety;
102 
103     /** Minimal reduction factor for stepsize control. */
104     private T minReduction;
105 
106     /** Maximal growth factor for stepsize control. */
107     private T maxGrowth;
108 
109     /** Flag setting whether coefficients in Butcher array are interpreted as Field or real numbers. */
110     private boolean usingFieldCoefficients;
111 
112     /** Build a Runge-Kutta integrator with the given Butcher array.
113      * @param field field to which the time and state vector elements belong
114      * @param name name of the method
115      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
116      * or -1 if method is not <i>fsal</i>
117      * @param minStep minimal step (sign is irrelevant, regardless of
118      * integration direction, forward or backward), the last step can
119      * be smaller than this
120      * @param maxStep maximal step (sign is irrelevant, regardless of
121      * integration direction, forward or backward), the last step can
122      * be smaller than this
123      * @param scalAbsoluteTolerance allowed absolute error
124      * @param scalRelativeTolerance allowed relative error
125      */
126     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
127                                                 final double minStep, final double maxStep,
128                                                 final double scalAbsoluteTolerance,
129                                                 final double scalRelativeTolerance) {
130 
131         super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
132 
133         this.fsal = fsal;
134         this.c    = getC();
135         this.a    = getA();
136         this.b    = getB();
137 
138         exp = -1.0 / getOrder();
139 
140         // set the default values of the algorithm control parameters
141         setSafety(field.getZero().add(0.9));
142         setMinReduction(field.getZero().add(0.2));
143         setMaxGrowth(field.getZero().add(10.0));
144 
145     }
146 
147     /** Build a Runge-Kutta integrator with the given Butcher array.
148      * @param field field to which the time and state vector elements belong
149      * @param name name of the method
150      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
151      * or -1 if method is not <i>fsal</i>
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 EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
160                                                 final double   minStep, final double maxStep,
161                                                 final double[] vecAbsoluteTolerance,
162                                                 final double[] vecRelativeTolerance) {
163 
164         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
165         this.usingFieldCoefficients = false;
166 
167         this.fsal = fsal;
168         this.c    = getC();
169         this.a    = getA();
170         this.b    = getB();
171 
172         exp = -1.0 / getOrder();
173 
174         // set the default values of the algorithm control parameters
175         setSafety(field.getZero().add(0.9));
176         setMinReduction(field.getZero().add(0.2));
177         setMaxGrowth(field.getZero().add(10.0));
178 
179     }
180 
181     /** Create an interpolator.
182      * @param forward integration direction indicator
183      * @param yDotK slopes at the intermediate points
184      * @param globalPreviousState start of the global step
185      * @param globalCurrentState end of the global step
186      * @param mapper equations mapper for the all equations
187      * @return external weights for the high order method from Butcher array
188      */
189     protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
190                                                                               FieldODEStateAndDerivative<T> globalPreviousState,
191                                                                               FieldODEStateAndDerivative<T> globalCurrentState,
192                                                                               FieldEquationsMapper<T> mapper);
193 
194     /** Get the order of the method.
195      * @return order of the method
196      */
197     public abstract int getOrder();
198 
199     /** Get the safety factor for stepsize control.
200      * @return safety factor
201      */
202     public T getSafety() {
203         return safety;
204     }
205 
206     /** Set the safety factor for stepsize control.
207      * @param safety safety factor
208      */
209     public void setSafety(final T safety) {
210         this.safety = safety;
211     }
212 
213     /**
214      * Setter for the flag between real or Field coefficients in the Butcher array.
215      *
216      * @param usingFieldCoefficients new value for flag
217      */
218     public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
219         this.usingFieldCoefficients = usingFieldCoefficients;
220     }
221 
222     /** {@inheritDoc} */
223     @Override
224     public boolean isUsingFieldCoefficients() {
225         return usingFieldCoefficients;
226     }
227 
228     /** {@inheritDoc} */
229     @Override
230     public int getNumberOfStages() {
231         return b.length;
232     }
233 
234     /** {@inheritDoc} */
235     @Override
236     protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
237         if (!isUsingFieldCoefficients()) {
238             realA = getRealA();
239             realB = getRealB();
240             realC = getRealC();
241         }
242         return super.initIntegration(eqn, s0, t);
243     }
244 
245     /** {@inheritDoc} */
246     @Override
247     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
248                                                    final FieldODEState<T> initialState, final T finalTime)
249         throws MathIllegalArgumentException, MathIllegalStateException {
250 
251         sanityChecks(initialState, finalTime);
252         setStepStart(initIntegration(equations, initialState, finalTime));
253         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
254 
255         // create some internal working arrays
256         final int   stages = getNumberOfStages();
257         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
258         T[]   yTmp   = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
259 
260         // set up integration control objects
261         T  hNew           = getField().getZero();
262         boolean firstTime = true;
263 
264         // main integration loop
265         setIsLastStep(false);
266         do {
267 
268             // iterate over step size, ensuring local normalized error is smaller than 1
269             double error = 10.0;
270             while (error >= 1.0) {
271 
272                 // first stage
273                 final T[] y = getStepStart().getCompleteState();
274                 yDotK[0] = getStepStart().getCompleteDerivative();
275 
276                 if (firstTime) {
277                     final StepsizeHelper helper = getStepSizeHelper();
278                     final T[] scale = MathArrays.buildArray(getField(), helper.getMainSetDimension());
279                     for (int i = 0; i < scale.length; ++i) {
280                         scale[i] = helper.getTolerance(i, y[i].abs());
281                     }
282                     hNew = getField().getZero().add(initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper()));
283                     firstTime = false;
284                 }
285 
286                 setStepSize(hNew);
287                 if (forward) {
288                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
289                         setStepSize(finalTime.subtract(getStepStart().getTime()));
290                     }
291                 } else {
292                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
293                         setStepSize(finalTime.subtract(getStepStart().getTime()));
294                     }
295                 }
296 
297                 // next stages
298                 if (isUsingFieldCoefficients()) {
299                     FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
300                             getStepStart().getTime(), y, getStepSize(), a, c, yDotK);
301                     yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
302                 } else {
303                     FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
304                             getStepStart().getTime(), y, getStepSize(), realA, realC, yDotK);
305                     yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
306                 }
307 
308                 incrementEvaluations(stages - 1);
309 
310                 // estimate the error at the end of the step
311                 error = estimateError(yDotK, y, yTmp, getStepSize());
312                 if (error >= 1.0) {
313                     // reject the step and attempt to reduce error by stepsize control
314                     final T factor = MathUtils.min(maxGrowth,
315                                                    MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
316                     hNew = getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false);
317                 }
318 
319             }
320             final T   stepEnd = getStepStart().getTime().add(getStepSize());
321             final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
322             final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
323 
324             // local error is small enough: accept the step, trigger events and step handlers
325             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
326                                     finalTime));
327 
328             if (!isLastStep()) {
329 
330                 // stepsize control for next step
331                 final T factor = MathUtils.min(maxGrowth,
332                                                MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
333                 final T  scaledH    = getStepSize().multiply(factor);
334                 final T  nextT      = getStepStart().getTime().add(scaledH);
335                 final boolean nextIsLast = forward ?
336                                            nextT.subtract(finalTime).getReal() >= 0 :
337                                            nextT.subtract(finalTime).getReal() <= 0;
338                 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
339 
340                 final T  filteredNextT      = getStepStart().getTime().add(hNew);
341                 final boolean filteredNextIsLast = forward ?
342                                                    filteredNextT.subtract(finalTime).getReal() >= 0 :
343                                                    filteredNextT.subtract(finalTime).getReal() <= 0;
344                 if (filteredNextIsLast) {
345                     hNew = finalTime.subtract(getStepStart().getTime());
346                 }
347 
348             }
349 
350         } while (!isLastStep());
351 
352         final FieldODEStateAndDerivative<T> finalState = getStepStart();
353         resetInternalState();
354         return finalState;
355 
356     }
357 
358     /** Get the minimal reduction factor for stepsize control.
359      * @return minimal reduction factor
360      */
361     public T getMinReduction() {
362         return minReduction;
363     }
364 
365     /** Set the minimal reduction factor for stepsize control.
366      * @param minReduction minimal reduction factor
367      */
368     public void setMinReduction(final T minReduction) {
369         this.minReduction = minReduction;
370     }
371 
372     /** Get the maximal growth factor for stepsize control.
373      * @return maximal growth factor
374      */
375     public T getMaxGrowth() {
376         return maxGrowth;
377     }
378 
379     /** Set the maximal growth factor for stepsize control.
380      * @param maxGrowth maximal growth factor
381      */
382     public void setMaxGrowth(final T maxGrowth) {
383         this.maxGrowth = maxGrowth;
384     }
385 
386     /** Compute the error ratio.
387      * @param yDotK derivatives computed during the first stages
388      * @param y0 estimate of the step at the start of the step
389      * @param y1 estimate of the step at the end of the step
390      * @param h  current step
391      * @return error ratio, greater than 1 if step should be rejected
392      */
393     protected abstract double estimateError(T[][] yDotK, T[] y0, T[] y1, T h);
394 
395 }