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 java.util.Arrays;
021    
022    import org.apache.commons.math.ode.DerivativeException;
023    import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
024    import org.apache.commons.math.ode.IntegratorException;
025    import org.apache.commons.math.ode.events.CombinedEventsManager;
026    import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
027    import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
028    import org.apache.commons.math.ode.sampling.StepHandler;
029    
030    /**
031     * This class implements the common part of all embedded Runge-Kutta
032     * integrators for Ordinary Differential Equations.
033     *
034     * <p>These methods are embedded explicit Runge-Kutta methods with two
035     * sets of coefficients allowing to estimate the error, their Butcher
036     * arrays are as follows :
037     * <pre>
038     *    0  |
039     *   c2  | a21
040     *   c3  | a31  a32
041     *   ... |        ...
042     *   cs  | as1  as2  ...  ass-1
043     *       |--------------------------
044     *       |  b1   b2  ...   bs-1  bs
045     *       |  b'1  b'2 ...   b's-1 b's
046     * </pre>
047     * </p>
048     *
049     * <p>In fact, we rather use the array defined by ej = bj - b'j to
050     * compute directly the error rather than computing two estimates and
051     * then comparing them.</p>
052     *
053     * <p>Some methods are qualified as <i>fsal</i> (first same as last)
054     * methods. This means the last evaluation of the derivatives in one
055     * step is the same as the first in the next step. Then, this
056     * evaluation can be reused from one step to the next one and the cost
057     * of such a method is really s-1 evaluations despite the method still
058     * has s stages. This behaviour is true only for successful steps, if
059     * the step is rejected after the error estimation phase, no
060     * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
061     * asi = bi for all i.</p>
062     *
063     * @version $Revision: 786874 $ $Date: 2009-06-20 14:09:16 -0400 (Sat, 20 Jun 2009) $
064     * @since 1.2
065     */
066    
067    public abstract class EmbeddedRungeKuttaIntegrator
068      extends AdaptiveStepsizeIntegrator {
069    
070      /** Build a Runge-Kutta integrator with the given Butcher array.
071       * @param name name of the method
072       * @param fsal indicate that the method is an <i>fsal</i>
073       * @param c time steps from Butcher array (without the first zero)
074       * @param a internal weights from Butcher array (without the first empty row)
075       * @param b propagation weights for the high order method from Butcher array
076       * @param prototype prototype of the step interpolator to use
077       * @param minStep minimal step (must be positive even for backward
078       * integration), the last step can be smaller than this
079       * @param maxStep maximal step (must be positive even for backward
080       * integration)
081       * @param scalAbsoluteTolerance allowed absolute error
082       * @param scalRelativeTolerance allowed relative error
083       */
084      protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
085                                             final double[] c, final double[][] a, final double[] b,
086                                             final RungeKuttaStepInterpolator prototype,
087                                             final double minStep, final double maxStep,
088                                             final double scalAbsoluteTolerance,
089                                             final double scalRelativeTolerance) {
090    
091        super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
092    
093        this.fsal      = fsal;
094        this.c         = c;
095        this.a         = a;
096        this.b         = b;
097        this.prototype = prototype;
098    
099        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    }