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 FieldButcherArrayProvider<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      /** Stepsize control exponent. */
89      private final double exp;
90  
91      /** Safety factor for stepsize control. */
92      private T safety;
93  
94      /** Minimal reduction factor for stepsize control. */
95      private T minReduction;
96  
97      /** Maximal growth factor for stepsize control. */
98      private T maxGrowth;
99  
100     /** Build a Runge-Kutta integrator with the given Butcher array.
101      * @param field field to which the time and state vector elements belong
102      * @param name name of the method
103      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
104      * or -1 if method is not <i>fsal</i>
105      * @param minStep minimal step (sign is irrelevant, regardless of
106      * integration direction, forward or backward), the last step can
107      * be smaller than this
108      * @param maxStep maximal step (sign is irrelevant, regardless of
109      * integration direction, forward or backward), the last step can
110      * be smaller than this
111      * @param scalAbsoluteTolerance allowed absolute error
112      * @param scalRelativeTolerance allowed relative error
113      */
114     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
115                                                 final double minStep, final double maxStep,
116                                                 final double scalAbsoluteTolerance,
117                                                 final double scalRelativeTolerance) {
118 
119         super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
120 
121         this.fsal = fsal;
122         this.c    = getC();
123         this.a    = getA();
124         this.b    = getB();
125 
126         exp = -1.0 / getOrder();
127 
128         // set the default values of the algorithm control parameters
129         setSafety(field.getZero().add(0.9));
130         setMinReduction(field.getZero().add(0.2));
131         setMaxGrowth(field.getZero().add(10.0));
132 
133     }
134 
135     /** Build a Runge-Kutta integrator with the given Butcher array.
136      * @param field field to which the time and state vector elements belong
137      * @param name name of the method
138      * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
139      * or -1 if method is not <i>fsal</i>
140      * @param minStep minimal step (must be positive even for backward
141      * integration), the last step can be smaller than this
142      * @param maxStep maximal step (must be positive even for backward
143      * integration)
144      * @param vecAbsoluteTolerance allowed absolute error
145      * @param vecRelativeTolerance allowed relative error
146      */
147     protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
148                                                 final double   minStep, final double maxStep,
149                                                 final double[] vecAbsoluteTolerance,
150                                                 final double[] vecRelativeTolerance) {
151 
152         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
153 
154         this.fsal = fsal;
155         this.c    = getC();
156         this.a    = getA();
157         this.b    = getB();
158 
159         exp = -1.0 / getOrder();
160 
161         // set the default values of the algorithm control parameters
162         setSafety(field.getZero().add(0.9));
163         setMinReduction(field.getZero().add(0.2));
164         setMaxGrowth(field.getZero().add(10.0));
165 
166     }
167 
168     /** Create a fraction.
169      * @param p numerator
170      * @param q denominator
171      * @return p/q computed in the instance field
172      */
173     protected T fraction(final int p, final int q) {
174         return getField().getOne().multiply(p).divide(q);
175     }
176 
177     /** Create a fraction.
178      * @param p numerator
179      * @param q denominator
180      * @return p/q computed in the instance field
181      */
182     protected T fraction(final double p, final double q) {
183         return getField().getOne().multiply(p).divide(q);
184     }
185 
186     /** Create an interpolator.
187      * @param forward integration direction indicator
188      * @param yDotK slopes at the intermediate points
189      * @param globalPreviousState start of the global step
190      * @param globalCurrentState end of the global step
191      * @param mapper equations mapper for the all equations
192      * @return external weights for the high order method from Butcher array
193      */
194     protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
195                                                                               FieldODEStateAndDerivative<T> globalPreviousState,
196                                                                               FieldODEStateAndDerivative<T> globalCurrentState,
197                                                                               FieldEquationsMapper<T> mapper);
198     /** Get the order of the method.
199      * @return order of the method
200      */
201     public abstract int getOrder();
202 
203     /** Get the safety factor for stepsize control.
204      * @return safety factor
205      */
206     public T getSafety() {
207         return safety;
208     }
209 
210     /** Set the safety factor for stepsize control.
211      * @param safety safety factor
212      */
213     public void setSafety(final T safety) {
214         this.safety = safety;
215     }
216 
217     /** {@inheritDoc} */
218     @Override
219     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
220                                                    final FieldODEState<T> initialState, final T finalTime)
221         throws MathIllegalArgumentException, MathIllegalStateException {
222 
223         sanityChecks(initialState, finalTime);
224         setStepStart(initIntegration(equations, initialState, finalTime));
225         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
226 
227         // create some internal working arrays
228         final int   stages = c.length + 1;
229         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
230         final T[]   yTmp   = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
231 
232         // set up integration control objects
233         T  hNew           = getField().getZero();
234         boolean firstTime = true;
235 
236         // main integration loop
237         setIsLastStep(false);
238         do {
239 
240             // iterate over step size, ensuring local normalized error is smaller than 1
241             double error = 10.0;
242             while (error >= 1.0) {
243 
244                 // first stage
245                 final T[] y = getStepStart().getCompleteState();
246                 yDotK[0] = getStepStart().getCompleteDerivative();
247 
248                 if (firstTime) {
249                     final StepsizeHelper helper = getStepSizeHelper();
250                     final T[] scale = MathArrays.buildArray(getField(), helper.getMainSetDimension());
251                     for (int i = 0; i < scale.length; ++i) {
252                         scale[i] = helper.getTolerance(i, y[i].abs());
253                     }
254                     hNew = getField().getZero().add(initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper()));
255                     firstTime = false;
256                 }
257 
258                 setStepSize(hNew);
259                 if (forward) {
260                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
261                         setStepSize(finalTime.subtract(getStepStart().getTime()));
262                     }
263                 } else {
264                     if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
265                         setStepSize(finalTime.subtract(getStepStart().getTime()));
266                     }
267                 }
268 
269                 // next stages
270                 for (int k = 1; k < stages; ++k) {
271 
272                     for (int j = 0; j < y.length; ++j) {
273                         T sum = yDotK[0][j].multiply(a[k-1][0]);
274                         for (int l = 1; l < k; ++l) {
275                             sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
276                         }
277                         yTmp[j] = y[j].add(getStepSize().multiply(sum));
278                     }
279 
280                     yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp);
281 
282                 }
283 
284                 // estimate the state at the end of the step
285                 for (int j = 0; j < y.length; ++j) {
286                     T sum    = yDotK[0][j].multiply(b[0]);
287                     for (int l = 1; l < stages; ++l) {
288                         sum = sum.add(yDotK[l][j].multiply(b[l]));
289                     }
290                     yTmp[j] = y[j].add(getStepSize().multiply(sum));
291                 }
292 
293                 // estimate the error at the end of the step
294                 error = estimateError(yDotK, y, yTmp, getStepSize());
295                 if (error >= 1.0) {
296                     // reject the step and attempt to reduce error by stepsize control
297                     final T factor = MathUtils.min(maxGrowth,
298                                                    MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
299                     hNew = getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false);
300                 }
301 
302             }
303             final T   stepEnd = getStepStart().getTime().add(getStepSize());
304             final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
305             final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
306 
307             // local error is small enough: accept the step, trigger events and step handlers
308             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
309                                     finalTime));
310 
311             if (!isLastStep()) {
312 
313                 // stepsize control for next step
314                 final T factor = MathUtils.min(maxGrowth,
315                                                MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
316                 final T  scaledH    = getStepSize().multiply(factor);
317                 final T  nextT      = getStepStart().getTime().add(scaledH);
318                 final boolean nextIsLast = forward ?
319                                            nextT.subtract(finalTime).getReal() >= 0 :
320                                            nextT.subtract(finalTime).getReal() <= 0;
321                 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
322 
323                 final T  filteredNextT      = getStepStart().getTime().add(hNew);
324                 final boolean filteredNextIsLast = forward ?
325                                                    filteredNextT.subtract(finalTime).getReal() >= 0 :
326                                                    filteredNextT.subtract(finalTime).getReal() <= 0;
327                 if (filteredNextIsLast) {
328                     hNew = finalTime.subtract(getStepStart().getTime());
329                 }
330 
331             }
332 
333         } while (!isLastStep());
334 
335         final FieldODEStateAndDerivative<T> finalState = getStepStart();
336         resetInternalState();
337         return finalState;
338 
339     }
340 
341     /** Get the minimal reduction factor for stepsize control.
342      * @return minimal reduction factor
343      */
344     public T getMinReduction() {
345         return minReduction;
346     }
347 
348     /** Set the minimal reduction factor for stepsize control.
349      * @param minReduction minimal reduction factor
350      */
351     public void setMinReduction(final T minReduction) {
352         this.minReduction = minReduction;
353     }
354 
355     /** Get the maximal growth factor for stepsize control.
356      * @return maximal growth factor
357      */
358     public T getMaxGrowth() {
359         return maxGrowth;
360     }
361 
362     /** Set the maximal growth factor for stepsize control.
363      * @param maxGrowth maximal growth factor
364      */
365     public void setMaxGrowth(final T maxGrowth) {
366         this.maxGrowth = maxGrowth;
367     }
368 
369     /** Compute the error ratio.
370      * @param yDotK derivatives computed during the first stages
371      * @param y0 estimate of the step at the start of the step
372      * @param y1 estimate of the step at the end of the step
373      * @param h  current step
374      * @return error ratio, greater than 1 if step should be rejected
375      */
376     protected abstract double estimateError(T[][] yDotK, T[] y0, T[] y1, T h);
377 
378 }