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.math.ode.nonstiff;
19  
20  import java.util.Arrays;
21  
22  import org.apache.commons.math.ode.DerivativeException;
23  import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
24  import org.apache.commons.math.ode.IntegratorException;
25  import org.apache.commons.math.ode.events.CombinedEventsManager;
26  import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
27  import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
28  import org.apache.commons.math.ode.sampling.StepHandler;
29  
30  /**
31   * This class implements the common part of all embedded Runge-Kutta
32   * integrators for Ordinary Differential Equations.
33   *
34   * <p>These methods are embedded explicit Runge-Kutta methods with two
35   * sets of coefficients allowing to estimate the error, their Butcher
36   * arrays are as follows :
37   * <pre>
38   *    0  |
39   *   c2  | a21
40   *   c3  | a31  a32
41   *   ... |        ...
42   *   cs  | as1  as2  ...  ass-1
43   *       |--------------------------
44   *       |  b1   b2  ...   bs-1  bs
45   *       |  b'1  b'2 ...   b's-1 b's
46   * </pre>
47   * </p>
48   *
49   * <p>In fact, we rather use the array defined by ej = bj - b'j to
50   * compute directly the error rather than computing two estimates and
51   * then comparing them.</p>
52   *
53   * <p>Some methods are qualified as <i>fsal</i> (first same as last)
54   * methods. This means the last evaluation of the derivatives in one
55   * step is the same as the first in the next step. Then, this
56   * evaluation can be reused from one step to the next one and the cost
57   * of such a method is really s-1 evaluations despite the method still
58   * has s stages. This behaviour is true only for successful steps, if
59   * the step is rejected after the error estimation phase, no
60   * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
61   * asi = bi for all i.</p>
62   *
63   * @version $Revision: 786874 $ $Date: 2009-06-20 14:09:16 -0400 (Sat, 20 Jun 2009) $
64   * @since 1.2
65   */
66  
67  public abstract class EmbeddedRungeKuttaIntegrator
68    extends AdaptiveStepsizeIntegrator {
69  
70    /** Build a Runge-Kutta integrator with the given Butcher array.
71     * @param name name of the method
72     * @param fsal indicate that the method is an <i>fsal</i>
73     * @param c time steps from Butcher array (without the first zero)
74     * @param a internal weights from Butcher array (without the first empty row)
75     * @param b propagation weights for the high order method from Butcher array
76     * @param prototype prototype of the step interpolator to use
77     * @param minStep minimal step (must be positive even for backward
78     * integration), the last step can be smaller than this
79     * @param maxStep maximal step (must be positive even for backward
80     * integration)
81     * @param scalAbsoluteTolerance allowed absolute error
82     * @param scalRelativeTolerance allowed relative error
83     */
84    protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
85                                           final double[] c, final double[][] a, final double[] b,
86                                           final RungeKuttaStepInterpolator prototype,
87                                           final double minStep, final double maxStep,
88                                           final double scalAbsoluteTolerance,
89                                           final double scalRelativeTolerance) {
90  
91      super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
92  
93      this.fsal      = fsal;
94      this.c         = c;
95      this.a         = a;
96      this.b         = b;
97      this.prototype = prototype;
98  
99      exp = -1.0 / getOrder();
100 
101     // set the default values of the algorithm control parameters
102     setSafety(0.9);
103     setMinReduction(0.2);
104     setMaxGrowth(10.0);
105 
106   }
107 
108   /** Build a Runge-Kutta integrator with the given Butcher array.
109    * @param name name of the method
110    * @param fsal indicate that the method is an <i>fsal</i>
111    * @param c time steps from Butcher array (without the first zero)
112    * @param a internal weights from Butcher array (without the first empty row)
113    * @param b propagation weights for the high order method from Butcher array
114    * @param prototype prototype of the step interpolator to use
115    * @param minStep minimal step (must be positive even for backward
116    * integration), the last step can be smaller than this
117    * @param maxStep maximal step (must be positive even for backward
118    * integration)
119    * @param vecAbsoluteTolerance allowed absolute error
120    * @param vecRelativeTolerance allowed relative error
121    */
122   protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
123                                          final double[] c, final double[][] a, final double[] b,
124                                          final RungeKuttaStepInterpolator prototype,
125                                          final double   minStep, final double maxStep,
126                                          final double[] vecAbsoluteTolerance,
127                                          final double[] vecRelativeTolerance) {
128 
129     super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
130 
131     this.fsal      = fsal;
132     this.c         = c;
133     this.a         = a;
134     this.b         = b;
135     this.prototype = prototype;
136 
137     exp = -1.0 / getOrder();
138 
139     // set the default values of the algorithm control parameters
140     setSafety(0.9);
141     setMinReduction(0.2);
142     setMaxGrowth(10.0);
143 
144   }
145 
146   /** Get the order of the method.
147    * @return order of the method
148    */
149   public abstract int getOrder();
150 
151   /** Get the safety factor for stepsize control.
152    * @return safety factor
153    */
154   public double getSafety() {
155     return safety;
156   }
157 
158   /** Set the safety factor for stepsize control.
159    * @param safety safety factor
160    */
161   public void setSafety(final double safety) {
162     this.safety = safety;
163   }
164 
165   /** {@inheritDoc} */
166   @Override
167   public double integrate(final FirstOrderDifferentialEquations equations,
168                           final double t0, final double[] y0,
169                           final double t, final double[] y)
170   throws DerivativeException, IntegratorException {
171 
172     sanityChecks(equations, t0, y0, t, y);
173     setEquations(equations);
174     resetEvaluations();
175     final boolean forward = (t > t0);
176 
177     // create some internal working arrays
178     final int stages = c.length + 1;
179     if (y != y0) {
180       System.arraycopy(y0, 0, y, 0, y0.length);
181     }
182     final double[][] yDotK = new double[stages][y0.length];
183     final double[] yTmp = new double[y0.length];
184 
185     // set up an interpolator sharing the integrator arrays
186     AbstractStepInterpolator interpolator;
187     if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
188       final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
189       rki.reinitialize(this, yTmp, yDotK, forward);
190       interpolator = rki;
191     } else {
192       interpolator = new DummyStepInterpolator(yTmp, forward);
193     }
194     interpolator.storeTime(t0);
195 
196     // set up integration control objects
197     stepStart         = t0;
198     double  hNew      = 0;
199     boolean firstTime = true;
200     for (StepHandler handler : stepHandlers) {
201         handler.reset();
202     }
203     CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
204     boolean lastStep = false;
205 
206     // main integration loop
207     while (!lastStep) {
208 
209       interpolator.shift();
210 
211       double error = 0;
212       for (boolean loop = true; loop;) {
213 
214         if (firstTime || !fsal) {
215           // first stage
216           computeDerivatives(stepStart, y, yDotK[0]);
217         }
218 
219         if (firstTime) {
220           final double[] scale;
221           if (vecAbsoluteTolerance != null) {
222             scale = vecAbsoluteTolerance;
223           } else {
224             scale = new double[y0.length];
225             Arrays.fill(scale, scalAbsoluteTolerance);
226           }
227           hNew = initializeStep(equations, forward, getOrder(), scale,
228                                 stepStart, y, yDotK[0], yTmp, yDotK[1]);
229           firstTime = false;
230         }
231 
232         stepSize = hNew;
233 
234         // next stages
235         for (int k = 1; k < stages; ++k) {
236 
237           for (int j = 0; j < y0.length; ++j) {
238             double sum = a[k-1][0] * yDotK[0][j];
239             for (int l = 1; l < k; ++l) {
240               sum += a[k-1][l] * yDotK[l][j];
241             }
242             yTmp[j] = y[j] + stepSize * sum;
243           }
244 
245           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
246 
247         }
248 
249         // estimate the state at the end of the step
250         for (int j = 0; j < y0.length; ++j) {
251           double sum    = b[0] * yDotK[0][j];
252           for (int l = 1; l < stages; ++l) {
253             sum    += b[l] * yDotK[l][j];
254           }
255           yTmp[j] = y[j] + stepSize * sum;
256         }
257 
258         // estimate the error at the end of the step
259         error = estimateError(yDotK, y, yTmp, stepSize);
260         if (error <= 1.0) {
261 
262           // discrete events handling
263           interpolator.storeTime(stepStart + stepSize);
264           if (manager.evaluateStep(interpolator)) {
265               final double dt = manager.getEventTime() - stepStart;
266               if (Math.abs(dt) <= Math.ulp(stepStart)) {
267                   // rejecting the step would lead to a too small next step, we accept it
268                   loop = false;
269               } else {
270                   // reject the step to match exactly the next switch time
271                   hNew = dt;
272               }
273           } else {
274             // accept the step
275             loop = false;
276           }
277 
278         } else {
279           // reject the step and attempt to reduce error by stepsize control
280           final double factor =
281               Math.min(maxGrowth,
282                        Math.max(minReduction, safety * Math.pow(error, exp)));
283           hNew = filterStep(stepSize * factor, forward, false);
284         }
285 
286       }
287 
288       // the step has been accepted
289       final double nextStep = stepStart + stepSize;
290       System.arraycopy(yTmp, 0, y, 0, y0.length);
291       manager.stepAccepted(nextStep, y);
292       lastStep = manager.stop();
293 
294       // provide the step data to the step handler
295       interpolator.storeTime(nextStep);
296       for (StepHandler handler : stepHandlers) {
297           handler.handleStep(interpolator, lastStep);
298       }
299       stepStart = nextStep;
300 
301       if (fsal) {
302         // save the last evaluation for the next step
303         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
304       }
305 
306       if (manager.reset(stepStart, y) && ! lastStep) {
307         // some event handler has triggered changes that
308         // invalidate the derivatives, we need to recompute them
309         computeDerivatives(stepStart, y, yDotK[0]);
310       }
311 
312       if (! lastStep) {
313         // in some rare cases we may get here with stepSize = 0, for example
314         // when an event occurs at integration start, reducing the first step
315         // to zero; we have to reset the step to some safe non zero value
316           stepSize = filterStep(stepSize, forward, true);
317 
318         // stepsize control for next step
319         final double factor = Math.min(maxGrowth,
320                                        Math.max(minReduction,
321                                                 safety * Math.pow(error, exp)));
322         final double  scaledH    = stepSize * factor;
323         final double  nextT      = stepStart + scaledH;
324         final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
325         hNew = filterStep(scaledH, forward, nextIsLast);
326       }
327 
328     }
329 
330     final double stopTime = stepStart;
331     resetInternalState();
332     return stopTime;
333 
334   }
335 
336   /** Get the minimal reduction factor for stepsize control.
337    * @return minimal reduction factor
338    */
339   public double getMinReduction() {
340     return minReduction;
341   }
342 
343   /** Set the minimal reduction factor for stepsize control.
344    * @param minReduction minimal reduction factor
345    */
346   public void setMinReduction(final double minReduction) {
347     this.minReduction = minReduction;
348   }
349 
350   /** Get the maximal growth factor for stepsize control.
351    * @return maximal growth factor
352    */
353   public double getMaxGrowth() {
354     return maxGrowth;
355   }
356 
357   /** Set the maximal growth factor for stepsize control.
358    * @param maxGrowth maximal growth factor
359    */
360   public void setMaxGrowth(final double maxGrowth) {
361     this.maxGrowth = maxGrowth;
362   }
363 
364   /** Compute the error ratio.
365    * @param yDotK derivatives computed during the first stages
366    * @param y0 estimate of the step at the start of the step
367    * @param y1 estimate of the step at the end of the step
368    * @param h  current step
369    * @return error ratio, greater than 1 if step should be rejected
370    */
371   protected abstract double estimateError(double[][] yDotK,
372                                           double[] y0, double[] y1,
373                                           double h);
374 
375   /** Indicator for <i>fsal</i> methods. */
376   private boolean fsal;
377 
378   /** Time steps from Butcher array (without the first zero). */
379   private double[] c;
380 
381   /** Internal weights from Butcher array (without the first empty row). */
382   private double[][] a;
383 
384   /** External weights for the high order method from Butcher array. */
385   private double[] b;
386 
387   /** Prototype of the step interpolator. */
388   private RungeKuttaStepInterpolator prototype;
389                                          
390   /** Stepsize control exponent. */
391   private double exp;
392 
393   /** Safety factor for stepsize control. */
394   private double safety;
395 
396   /** Minimal reduction factor for stepsize control. */
397   private double minReduction;
398 
399   /** Maximal growth factor for stepsize control. */
400   private double maxGrowth;
401 
402 }