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.nonstiff;
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.ode.ExpandableStatefulODE;
25  import org.apache.commons.math3.util.FastMath;
26  
27  /**
28   * This class implements the common part of all embedded Runge-Kutta
29   * integrators for Ordinary Differential Equations.
30   *
31   * <p>These methods are embedded explicit Runge-Kutta methods with two
32   * sets of coefficients allowing to estimate the error, their Butcher
33   * arrays are as follows :
34   * <pre>
35   *    0  |
36   *   c2  | a21
37   *   c3  | a31  a32
38   *   ... |        ...
39   *   cs  | as1  as2  ...  ass-1
40   *       |--------------------------
41   *       |  b1   b2  ...   bs-1  bs
42   *       |  b'1  b'2 ...   b's-1 b's
43   * </pre>
44   * </p>
45   *
46   * <p>In fact, we rather use the array defined by ej = bj - b'j to
47   * compute directly the error rather than computing two estimates and
48   * then comparing them.</p>
49   *
50   * <p>Some methods are qualified as <i>fsal</i> (first same as last)
51   * methods. This means the last evaluation of the derivatives in one
52   * step is the same as the first in the next step. Then, this
53   * evaluation can be reused from one step to the next one and the cost
54   * of such a method is really s-1 evaluations despite the method still
55   * has s stages. This behaviour is true only for successful steps, if
56   * the step is rejected after the error estimation phase, no
57   * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
58   * asi = bi for all i.</p>
59   *
60   * @since 1.2
61   */
62  
63  public abstract class EmbeddedRungeKuttaIntegrator
64    extends AdaptiveStepsizeIntegrator {
65  
66      /** Indicator for <i>fsal</i> methods. */
67      private final boolean fsal;
68  
69      /** Time steps from Butcher array (without the first zero). */
70      private final double[] c;
71  
72      /** Internal weights from Butcher array (without the first empty row). */
73      private final double[][] a;
74  
75      /** External weights for the high order method from Butcher array. */
76      private final double[] b;
77  
78      /** Prototype of the step interpolator. */
79      private final RungeKuttaStepInterpolator prototype;
80  
81      /** Stepsize control exponent. */
82      private final double exp;
83  
84      /** Safety factor for stepsize control. */
85      private double safety;
86  
87      /** Minimal reduction factor for stepsize control. */
88      private double minReduction;
89  
90      /** Maximal growth factor for stepsize control. */
91      private double maxGrowth;
92  
93    /** Build a Runge-Kutta integrator with the given Butcher array.
94     * @param name name of the method
95     * @param fsal indicate that the method is an <i>fsal</i>
96     * @param c time steps from Butcher array (without the first zero)
97     * @param a internal weights from Butcher array (without the first empty row)
98     * @param b propagation weights for the high order method from Butcher array
99     * @param prototype prototype of the step interpolator to use
100    * @param minStep minimal step (sign is irrelevant, regardless of
101    * integration direction, forward or backward), the last step can
102    * be smaller than this
103    * @param maxStep maximal step (sign is irrelevant, regardless of
104    * integration direction, forward or backward), the last step can
105    * be smaller than this
106    * @param scalAbsoluteTolerance allowed absolute error
107    * @param scalRelativeTolerance allowed relative error
108    */
109   protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
110                                          final double[] c, final double[][] a, final double[] b,
111                                          final RungeKuttaStepInterpolator prototype,
112                                          final double minStep, final double maxStep,
113                                          final double scalAbsoluteTolerance,
114                                          final double scalRelativeTolerance) {
115 
116     super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
117 
118     this.fsal      = fsal;
119     this.c         = c;
120     this.a         = a;
121     this.b         = b;
122     this.prototype = prototype;
123 
124     exp = -1.0 / getOrder();
125 
126     // set the default values of the algorithm control parameters
127     setSafety(0.9);
128     setMinReduction(0.2);
129     setMaxGrowth(10.0);
130 
131   }
132 
133   /** Build a Runge-Kutta integrator with the given Butcher array.
134    * @param name name of the method
135    * @param fsal indicate that the method is an <i>fsal</i>
136    * @param c time steps from Butcher array (without the first zero)
137    * @param a internal weights from Butcher array (without the first empty row)
138    * @param b propagation weights for the high order method from Butcher array
139    * @param prototype prototype of the step interpolator to use
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 EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
148                                          final double[] c, final double[][] a, final double[] b,
149                                          final RungeKuttaStepInterpolator prototype,
150                                          final double   minStep, final double maxStep,
151                                          final double[] vecAbsoluteTolerance,
152                                          final double[] vecRelativeTolerance) {
153 
154     super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
155 
156     this.fsal      = fsal;
157     this.c         = c;
158     this.a         = a;
159     this.b         = b;
160     this.prototype = prototype;
161 
162     exp = -1.0 / getOrder();
163 
164     // set the default values of the algorithm control parameters
165     setSafety(0.9);
166     setMinReduction(0.2);
167     setMaxGrowth(10.0);
168 
169   }
170 
171   /** Get the order of the method.
172    * @return order of the method
173    */
174   public abstract int getOrder();
175 
176   /** Get the safety factor for stepsize control.
177    * @return safety factor
178    */
179   public double getSafety() {
180     return safety;
181   }
182 
183   /** Set the safety factor for stepsize control.
184    * @param safety safety factor
185    */
186   public void setSafety(final double safety) {
187     this.safety = safety;
188   }
189 
190   /** {@inheritDoc} */
191   @Override
192   public void integrate(final ExpandableStatefulODE equations, final double t)
193       throws NumberIsTooSmallException, DimensionMismatchException,
194              MaxCountExceededException, NoBracketingException {
195 
196     sanityChecks(equations, t);
197     setEquations(equations);
198     final boolean forward = t > equations.getTime();
199 
200     // create some internal working arrays
201     final double[] y0  = equations.getCompleteState();
202     final double[] y = y0.clone();
203     final int stages = c.length + 1;
204     final double[][] yDotK = new double[stages][y.length];
205     final double[] yTmp    = y0.clone();
206     final double[] yDotTmp = new double[y.length];
207 
208     // set up an interpolator sharing the integrator arrays
209     final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
210     interpolator.reinitialize(this, yTmp, yDotK, forward,
211                               equations.getPrimaryMapper(), equations.getSecondaryMappers());
212     interpolator.storeTime(equations.getTime());
213 
214     // set up integration control objects
215     stepStart         = equations.getTime();
216     double  hNew      = 0;
217     boolean firstTime = true;
218     initIntegration(equations.getTime(), y0, t);
219 
220     // main integration loop
221     isLastStep = false;
222     do {
223 
224       interpolator.shift();
225 
226       // iterate over step size, ensuring local normalized error is smaller than 1
227       double error = 10;
228       while (error >= 1.0) {
229 
230         if (firstTime || !fsal) {
231           // first stage
232           computeDerivatives(stepStart, y, yDotK[0]);
233         }
234 
235         if (firstTime) {
236           final double[] scale = new double[mainSetDimension];
237           if (vecAbsoluteTolerance == null) {
238               for (int i = 0; i < scale.length; ++i) {
239                 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]);
240               }
241           } else {
242               for (int i = 0; i < scale.length; ++i) {
243                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
244               }
245           }
246           hNew = initializeStep(forward, getOrder(), scale,
247                                 stepStart, y, yDotK[0], yTmp, yDotK[1]);
248           firstTime = false;
249         }
250 
251         stepSize = hNew;
252         if (forward) {
253             if (stepStart + stepSize >= t) {
254                 stepSize = t - stepStart;
255             }
256         } else {
257             if (stepStart + stepSize <= t) {
258                 stepSize = t - stepStart;
259             }
260         }
261 
262         // next stages
263         for (int k = 1; k < stages; ++k) {
264 
265           for (int j = 0; j < y0.length; ++j) {
266             double sum = a[k-1][0] * yDotK[0][j];
267             for (int l = 1; l < k; ++l) {
268               sum += a[k-1][l] * yDotK[l][j];
269             }
270             yTmp[j] = y[j] + stepSize * sum;
271           }
272 
273           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
274 
275         }
276 
277         // estimate the state at the end of the step
278         for (int j = 0; j < y0.length; ++j) {
279           double sum    = b[0] * yDotK[0][j];
280           for (int l = 1; l < stages; ++l) {
281             sum    += b[l] * yDotK[l][j];
282           }
283           yTmp[j] = y[j] + stepSize * sum;
284         }
285 
286         // estimate the error at the end of the step
287         error = estimateError(yDotK, y, yTmp, stepSize);
288         if (error >= 1.0) {
289           // reject the step and attempt to reduce error by stepsize control
290           final double factor =
291               FastMath.min(maxGrowth,
292                            FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
293           hNew = filterStep(stepSize * factor, forward, false);
294         }
295 
296       }
297 
298       // local error is small enough: accept the step, trigger events and step handlers
299       interpolator.storeTime(stepStart + stepSize);
300       System.arraycopy(yTmp, 0, y, 0, y0.length);
301       System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
302       stepStart = acceptStep(interpolator, y, yDotTmp, t);
303       System.arraycopy(y, 0, yTmp, 0, y.length);
304 
305       if (!isLastStep) {
306 
307           // prepare next step
308           interpolator.storeTime(stepStart);
309 
310           if (fsal) {
311               // save the last evaluation for the next step
312               System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length);
313           }
314 
315           // stepsize control for next step
316           final double factor =
317               FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
318           final double  scaledH    = stepSize * factor;
319           final double  nextT      = stepStart + scaledH;
320           final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
321           hNew = filterStep(scaledH, forward, nextIsLast);
322 
323           final double  filteredNextT      = stepStart + hNew;
324           final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
325           if (filteredNextIsLast) {
326               hNew = t - stepStart;
327           }
328 
329       }
330 
331     } while (!isLastStep);
332 
333     // dispatch results
334     equations.setTime(stepStart);
335     equations.setCompleteState(y);
336 
337     resetInternalState();
338 
339   }
340 
341   /** Get the minimal reduction factor for stepsize control.
342    * @return minimal reduction factor
343    */
344   public double 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 double minReduction) {
352     this.minReduction = minReduction;
353   }
354 
355   /** Get the maximal growth factor for stepsize control.
356    * @return maximal growth factor
357    */
358   public double 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 double 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(double[][] yDotK,
377                                           double[] y0, double[] y1,
378                                           double h);
379 
380 }