001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.ode.nonstiff;
019    
020    import org.apache.commons.math.ode.DerivativeException;
021    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
022    import org.apache.commons.math.ode.IntegratorException;
023    import org.apache.commons.math.ode.events.CombinedEventsManager;
024    import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
025    import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
026    import org.apache.commons.math.ode.sampling.StepHandler;
027    
028    /**
029     * This class implements the common part of all embedded Runge-Kutta
030     * integrators for Ordinary Differential Equations.
031     *
032     * <p>These methods are embedded explicit Runge-Kutta methods with two
033     * sets of coefficients allowing to estimate the error, their Butcher
034     * arrays are as follows :
035     * <pre>
036     *    0  |
037     *   c2  | a21
038     *   c3  | a31  a32
039     *   ... |        ...
040     *   cs  | as1  as2  ...  ass-1
041     *       |--------------------------
042     *       |  b1   b2  ...   bs-1  bs
043     *       |  b'1  b'2 ...   b's-1 b's
044     * </pre>
045     * </p>
046     *
047     * <p>In fact, we rather use the array defined by ej = bj - b'j to
048     * compute directly the error rather than computing two estimates and
049     * then comparing them.</p>
050     *
051     * <p>Some methods are qualified as <i>fsal</i> (first same as last)
052     * methods. This means the last evaluation of the derivatives in one
053     * step is the same as the first in the next step. Then, this
054     * evaluation can be reused from one step to the next one and the cost
055     * of such a method is really s-1 evaluations despite the method still
056     * has s stages. This behaviour is true only for successful steps, if
057     * the step is rejected after the error estimation phase, no
058     * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
059     * asi = bi for all i.</p>
060     *
061     * @version $Revision: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $
062     * @since 1.2
063     */
064    
065    public abstract class EmbeddedRungeKuttaIntegrator
066      extends AdaptiveStepsizeIntegrator {
067    
068        /** Indicator for <i>fsal</i> methods. */
069        private final boolean fsal;
070    
071        /** Time steps from Butcher array (without the first zero). */
072        private final double[] c;
073    
074        /** Internal weights from Butcher array (without the first empty row). */
075        private final double[][] a;
076    
077        /** External weights for the high order method from Butcher array. */
078        private final double[] b;
079    
080        /** Prototype of the step interpolator. */
081        private final RungeKuttaStepInterpolator prototype;
082    
083        /** Stepsize control exponent. */
084        private final double exp;
085    
086        /** Safety factor for stepsize control. */
087        private double safety;
088    
089        /** Minimal reduction factor for stepsize control. */
090        private double minReduction;
091    
092        /** Maximal growth factor for stepsize control. */
093        private double maxGrowth;
094    
095      /** Build a Runge-Kutta integrator with the given Butcher array.
096       * @param name name of the method
097       * @param fsal indicate that the method is an <i>fsal</i>
098       * @param c time steps from Butcher array (without the first zero)
099       * @param a internal weights from Butcher array (without the first empty row)
100       * @param b propagation weights for the high order method from Butcher array
101       * @param prototype prototype of the step interpolator to use
102       * @param minStep minimal step (must be positive even for backward
103       * integration), the last step can be smaller than this
104       * @param maxStep maximal step (must be positive even for backward
105       * integration)
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 double integrate(final FirstOrderDifferentialEquations equations,
193                              final double t0, final double[] y0,
194                              final double t, final double[] y)
195      throws DerivativeException, IntegratorException {
196    
197        sanityChecks(equations, t0, y0, t, y);
198        setEquations(equations);
199        resetEvaluations();
200        final boolean forward = t > t0;
201    
202        // create some internal working arrays
203        final int stages = c.length + 1;
204        if (y != y0) {
205          System.arraycopy(y0, 0, y, 0, y0.length);
206        }
207        final double[][] yDotK = new double[stages][y0.length];
208        final double[] yTmp = new double[y0.length];
209    
210        // set up an interpolator sharing the integrator arrays
211        AbstractStepInterpolator interpolator;
212        if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
213          final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
214          rki.reinitialize(this, yTmp, yDotK, forward);
215          interpolator = rki;
216        } else {
217          interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
218        }
219        interpolator.storeTime(t0);
220    
221        // set up integration control objects
222        stepStart         = t0;
223        double  hNew      = 0;
224        boolean firstTime = true;
225        for (StepHandler handler : stepHandlers) {
226            handler.reset();
227        }
228        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
229        boolean lastStep = false;
230    
231        // main integration loop
232        while (!lastStep) {
233    
234          interpolator.shift();
235    
236          double error = 0;
237          for (boolean loop = true; loop;) {
238    
239            if (firstTime || !fsal) {
240              // first stage
241              computeDerivatives(stepStart, y, yDotK[0]);
242            }
243    
244            if (firstTime) {
245              final double[] scale = new double[y0.length];
246              if (vecAbsoluteTolerance == null) {
247                  for (int i = 0; i < scale.length; ++i) {
248                    scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * Math.abs(y[i]);
249                  }
250                } else {
251                  for (int i = 0; i < scale.length; ++i) {
252                    scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * Math.abs(y[i]);
253                  }
254                }
255              hNew = initializeStep(equations, forward, getOrder(), scale,
256                                    stepStart, y, yDotK[0], yTmp, yDotK[1]);
257              firstTime = false;
258            }
259    
260            stepSize = hNew;
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    
290              // discrete events handling
291              interpolator.storeTime(stepStart + stepSize);
292              if (manager.evaluateStep(interpolator)) {
293                  final double dt = manager.getEventTime() - stepStart;
294                  if (Math.abs(dt) <= Math.ulp(stepStart)) {
295                      // we cannot simply truncate the step, reject the current computation
296                      // and let the loop compute another state with the truncated step.
297                      // it is so small (much probably exactly 0 due to limited accuracy)
298                      // that the code above would fail handling it.
299                      // So we set up an artificial 0 size step by copying states
300                      interpolator.storeTime(stepStart);
301                      System.arraycopy(y, 0, yTmp, 0, y0.length);
302                      hNew     = 0;
303                      stepSize = 0;
304                      loop     = false;
305                  } else {
306                      // reject the step to match exactly the next switch time
307                      hNew = dt;
308                  }
309              } else {
310                // accept the step
311                loop = false;
312              }
313    
314            } else {
315              // reject the step and attempt to reduce error by stepsize control
316              final double factor =
317                  Math.min(maxGrowth,
318                           Math.max(minReduction, safety * Math.pow(error, exp)));
319              hNew = filterStep(stepSize * factor, forward, false);
320            }
321    
322          }
323    
324          // the step has been accepted
325          final double nextStep = stepStart + stepSize;
326          System.arraycopy(yTmp, 0, y, 0, y0.length);
327          manager.stepAccepted(nextStep, y);
328          lastStep = manager.stop();
329    
330          // provide the step data to the step handler
331          interpolator.storeTime(nextStep);
332          for (StepHandler handler : stepHandlers) {
333              handler.handleStep(interpolator, lastStep);
334          }
335          stepStart = nextStep;
336    
337          if (fsal) {
338            // save the last evaluation for the next step
339            System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
340          }
341    
342          if (manager.reset(stepStart, y) && ! lastStep) {
343            // some event handler has triggered changes that
344            // invalidate the derivatives, we need to recompute them
345            computeDerivatives(stepStart, y, yDotK[0]);
346          }
347    
348          if (! lastStep) {
349            // in some rare cases we may get here with stepSize = 0, for example
350            // when an event occurs at integration start, reducing the first step
351            // to zero; we have to reset the step to some safe non zero value
352              stepSize = filterStep(stepSize, forward, true);
353    
354            // stepsize control for next step
355            final double factor = Math.min(maxGrowth,
356                                           Math.max(minReduction,
357                                                    safety * Math.pow(error, exp)));
358            final double  scaledH    = stepSize * factor;
359            final double  nextT      = stepStart + scaledH;
360            final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
361            hNew = filterStep(scaledH, forward, nextIsLast);
362          }
363    
364        }
365    
366        final double stopTime = stepStart;
367        resetInternalState();
368        return stopTime;
369    
370      }
371    
372      /** Get the minimal reduction factor for stepsize control.
373       * @return minimal reduction factor
374       */
375      public double getMinReduction() {
376        return minReduction;
377      }
378    
379      /** Set the minimal reduction factor for stepsize control.
380       * @param minReduction minimal reduction factor
381       */
382      public void setMinReduction(final double minReduction) {
383        this.minReduction = minReduction;
384      }
385    
386      /** Get the maximal growth factor for stepsize control.
387       * @return maximal growth factor
388       */
389      public double getMaxGrowth() {
390        return maxGrowth;
391      }
392    
393      /** Set the maximal growth factor for stepsize control.
394       * @param maxGrowth maximal growth factor
395       */
396      public void setMaxGrowth(final double maxGrowth) {
397        this.maxGrowth = maxGrowth;
398      }
399    
400      /** Compute the error ratio.
401       * @param yDotK derivatives computed during the first stages
402       * @param y0 estimate of the step at the start of the step
403       * @param y1 estimate of the step at the end of the step
404       * @param h  current step
405       * @return error ratio, greater than 1 if step should be rejected
406       */
407      protected abstract double estimateError(double[][] yDotK,
408                                              double[] y0, double[] y1,
409                                              double h);
410    
411    }