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   * @version $Id: EmbeddedRungeKuttaIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
61   * @since 1.2
62   */
63  
64  public abstract class EmbeddedRungeKuttaIntegrator
65    extends AdaptiveStepsizeIntegrator {
66  
67      /** Indicator for <i>fsal</i> methods. */
68      private final boolean fsal;
69  
70      /** Time steps from Butcher array (without the first zero). */
71      private final double[] c;
72  
73      /** Internal weights from Butcher array (without the first empty row). */
74      private final double[][] a;
75  
76      /** External weights for the high order method from Butcher array. */
77      private final double[] b;
78  
79      /** Prototype of the step interpolator. */
80      private final RungeKuttaStepInterpolator prototype;
81  
82      /** Stepsize control exponent. */
83      private final double exp;
84  
85      /** Safety factor for stepsize control. */
86      private double safety;
87  
88      /** Minimal reduction factor for stepsize control. */
89      private double minReduction;
90  
91      /** Maximal growth factor for stepsize control. */
92      private double maxGrowth;
93  
94    /** Build a Runge-Kutta integrator with the given Butcher array.
95     * @param name name of the method
96     * @param fsal indicate that the method is an <i>fsal</i>
97     * @param c time steps from Butcher array (without the first zero)
98     * @param a internal weights from Butcher array (without the first empty row)
99     * @param b propagation weights for the high order method from Butcher array
100    * @param prototype prototype of the step interpolator to use
101    * @param minStep minimal step (sign is irrelevant, regardless of
102    * integration direction, forward or backward), the last step can
103    * be smaller than this
104    * @param maxStep maximal step (sign is irrelevant, regardless of
105    * integration direction, forward or backward), the last step can
106    * be smaller than this
107    * @param scalAbsoluteTolerance allowed absolute error
108    * @param scalRelativeTolerance allowed relative error
109    */
110   protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
111                                          final double[] c, final double[][] a, final double[] b,
112                                          final RungeKuttaStepInterpolator prototype,
113                                          final double minStep, final double maxStep,
114                                          final double scalAbsoluteTolerance,
115                                          final double scalRelativeTolerance) {
116 
117     super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
118 
119     this.fsal      = fsal;
120     this.c         = c;
121     this.a         = a;
122     this.b         = b;
123     this.prototype = prototype;
124 
125     exp = -1.0 / getOrder();
126 
127     // set the default values of the algorithm control parameters
128     setSafety(0.9);
129     setMinReduction(0.2);
130     setMaxGrowth(10.0);
131 
132   }
133 
134   /** Build a Runge-Kutta integrator with the given Butcher array.
135    * @param name name of the method
136    * @param fsal indicate that the method is an <i>fsal</i>
137    * @param c time steps from Butcher array (without the first zero)
138    * @param a internal weights from Butcher array (without the first empty row)
139    * @param b propagation weights for the high order method from Butcher array
140    * @param prototype prototype of the step interpolator to use
141    * @param minStep minimal step (must be positive even for backward
142    * integration), the last step can be smaller than this
143    * @param maxStep maximal step (must be positive even for backward
144    * integration)
145    * @param vecAbsoluteTolerance allowed absolute error
146    * @param vecRelativeTolerance allowed relative error
147    */
148   protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
149                                          final double[] c, final double[][] a, final double[] b,
150                                          final RungeKuttaStepInterpolator prototype,
151                                          final double   minStep, final double maxStep,
152                                          final double[] vecAbsoluteTolerance,
153                                          final double[] vecRelativeTolerance) {
154 
155     super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
156 
157     this.fsal      = fsal;
158     this.c         = c;
159     this.a         = a;
160     this.b         = b;
161     this.prototype = prototype;
162 
163     exp = -1.0 / getOrder();
164 
165     // set the default values of the algorithm control parameters
166     setSafety(0.9);
167     setMinReduction(0.2);
168     setMaxGrowth(10.0);
169 
170   }
171 
172   /** Get the order of the method.
173    * @return order of the method
174    */
175   public abstract int getOrder();
176 
177   /** Get the safety factor for stepsize control.
178    * @return safety factor
179    */
180   public double getSafety() {
181     return safety;
182   }
183 
184   /** Set the safety factor for stepsize control.
185    * @param safety safety factor
186    */
187   public void setSafety(final double safety) {
188     this.safety = safety;
189   }
190 
191   /** {@inheritDoc} */
192   @Override
193   public void integrate(final ExpandableStatefulODE equations, final double t)
194       throws NumberIsTooSmallException, DimensionMismatchException,
195              MaxCountExceededException, NoBracketingException {
196 
197     sanityChecks(equations, t);
198     setEquations(equations);
199     final boolean forward = t > equations.getTime();
200 
201     // create some internal working arrays
202     final double[] y0  = equations.getCompleteState();
203     final double[] y = y0.clone();
204     final int stages = c.length + 1;
205     final double[][] yDotK = new double[stages][y.length];
206     final double[] yTmp    = y0.clone();
207     final double[] yDotTmp = new double[y.length];
208 
209     // set up an interpolator sharing the integrator arrays
210     final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
211     interpolator.reinitialize(this, yTmp, yDotK, forward,
212                               equations.getPrimaryMapper(), equations.getSecondaryMappers());
213     interpolator.storeTime(equations.getTime());
214 
215     // set up integration control objects
216     stepStart         = equations.getTime();
217     double  hNew      = 0;
218     boolean firstTime = true;
219     initIntegration(equations.getTime(), y0, t);
220 
221     // main integration loop
222     isLastStep = false;
223     do {
224 
225       interpolator.shift();
226 
227       // iterate over step size, ensuring local normalized error is smaller than 1
228       double error = 10;
229       while (error >= 1.0) {
230 
231         if (firstTime || !fsal) {
232           // first stage
233           computeDerivatives(stepStart, y, yDotK[0]);
234         }
235 
236         if (firstTime) {
237           final double[] scale = new double[mainSetDimension];
238           if (vecAbsoluteTolerance == null) {
239               for (int i = 0; i < scale.length; ++i) {
240                 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]);
241               }
242           } else {
243               for (int i = 0; i < scale.length; ++i) {
244                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
245               }
246           }
247           hNew = initializeStep(forward, getOrder(), scale,
248                                 stepStart, y, yDotK[0], yTmp, yDotK[1]);
249           firstTime = false;
250         }
251 
252         stepSize = hNew;
253         if (forward) {
254             if (stepStart + stepSize >= t) {
255                 stepSize = t - stepStart;
256             }
257         } else {
258             if (stepStart + stepSize <= t) {
259                 stepSize = t - stepStart;
260             }
261         }
262 
263         // next stages
264         for (int k = 1; k < stages; ++k) {
265 
266           for (int j = 0; j < y0.length; ++j) {
267             double sum = a[k-1][0] * yDotK[0][j];
268             for (int l = 1; l < k; ++l) {
269               sum += a[k-1][l] * yDotK[l][j];
270             }
271             yTmp[j] = y[j] + stepSize * sum;
272           }
273 
274           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
275 
276         }
277 
278         // estimate the state at the end of the step
279         for (int j = 0; j < y0.length; ++j) {
280           double sum    = b[0] * yDotK[0][j];
281           for (int l = 1; l < stages; ++l) {
282             sum    += b[l] * yDotK[l][j];
283           }
284           yTmp[j] = y[j] + stepSize * sum;
285         }
286 
287         // estimate the error at the end of the step
288         error = estimateError(yDotK, y, yTmp, stepSize);
289         if (error >= 1.0) {
290           // reject the step and attempt to reduce error by stepsize control
291           final double factor =
292               FastMath.min(maxGrowth,
293                            FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
294           hNew = filterStep(stepSize * factor, forward, false);
295         }
296 
297       }
298 
299       // local error is small enough: accept the step, trigger events and step handlers
300       interpolator.storeTime(stepStart + stepSize);
301       System.arraycopy(yTmp, 0, y, 0, y0.length);
302       System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
303       stepStart = acceptStep(interpolator, y, yDotTmp, t);
304       System.arraycopy(y, 0, yTmp, 0, y.length);
305 
306       if (!isLastStep) {
307 
308           // prepare next step
309           interpolator.storeTime(stepStart);
310 
311           if (fsal) {
312               // save the last evaluation for the next step
313               System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length);
314           }
315 
316           // stepsize control for next step
317           final double factor =
318               FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
319           final double  scaledH    = stepSize * factor;
320           final double  nextT      = stepStart + scaledH;
321           final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
322           hNew = filterStep(scaledH, forward, nextIsLast);
323 
324           final double  filteredNextT      = stepStart + hNew;
325           final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
326           if (filteredNextIsLast) {
327               hNew = t - stepStart;
328           }
329 
330       }
331 
332     } while (!isLastStep);
333 
334     // dispatch results
335     equations.setTime(stepStart);
336     equations.setCompleteState(y);
337 
338     resetInternalState();
339 
340   }
341 
342   /** Get the minimal reduction factor for stepsize control.
343    * @return minimal reduction factor
344    */
345   public double getMinReduction() {
346     return minReduction;
347   }
348 
349   /** Set the minimal reduction factor for stepsize control.
350    * @param minReduction minimal reduction factor
351    */
352   public void setMinReduction(final double minReduction) {
353     this.minReduction = minReduction;
354   }
355 
356   /** Get the maximal growth factor for stepsize control.
357    * @return maximal growth factor
358    */
359   public double getMaxGrowth() {
360     return maxGrowth;
361   }
362 
363   /** Set the maximal growth factor for stepsize control.
364    * @param maxGrowth maximal growth factor
365    */
366   public void setMaxGrowth(final double maxGrowth) {
367     this.maxGrowth = maxGrowth;
368   }
369 
370   /** Compute the error ratio.
371    * @param yDotK derivatives computed during the first stages
372    * @param y0 estimate of the step at the start of the step
373    * @param y1 estimate of the step at the end of the step
374    * @param h  current step
375    * @return error ratio, greater than 1 if step should be rejected
376    */
377   protected abstract double estimateError(double[][] yDotK,
378                                           double[] y0, double[] y1,
379                                           double h);
380 
381 }