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  
21  import org.apache.commons.math.ode.AbstractIntegrator;
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 fixed step Runge-Kutta
32   * integrators for Ordinary Differential Equations.
33   *
34   * <p>These methods are explicit Runge-Kutta methods, their Butcher
35   * arrays are as follows :
36   * <pre>
37   *    0  |
38   *   c2  | a21
39   *   c3  | a31  a32
40   *   ... |        ...
41   *   cs  | as1  as2  ...  ass-1
42   *       |--------------------------
43   *       |  b1   b2  ...   bs-1  bs
44   * </pre>
45   * </p>
46   *
47   * @see EulerIntegrator
48   * @see ClassicalRungeKuttaIntegrator
49   * @see GillIntegrator
50   * @see MidpointIntegrator
51   * @version $Revision: 785473 $ $Date: 2009-06-17 00:02:35 -0400 (Wed, 17 Jun 2009) $
52   * @since 1.2
53   */
54  
55  public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
56  
57    /** Simple constructor.
58     * Build a Runge-Kutta integrator with the given
59     * step. The default step handler does nothing.
60     * @param name name of the method
61     * @param c time steps from Butcher array (without the first zero)
62     * @param a internal weights from Butcher array (without the first empty row)
63     * @param b propagation weights for the high order method from Butcher array
64     * @param prototype prototype of the step interpolator to use
65     * @param step integration step
66     */
67    protected RungeKuttaIntegrator(final String name,
68                                   final double[] c, final double[][] a, final double[] b,
69                                   final RungeKuttaStepInterpolator prototype,
70                                   final double step) {
71      super(name);
72      this.c          = c;
73      this.a          = a;
74      this.b          = b;
75      this.prototype  = prototype;
76      this.step       = Math.abs(step);
77    }
78  
79    /** {@inheritDoc} */
80    public double integrate(final FirstOrderDifferentialEquations equations,
81                            final double t0, final double[] y0,
82                            final double t, final double[] y)
83    throws DerivativeException, IntegratorException {
84  
85      sanityChecks(equations, t0, y0, t, y);
86      setEquations(equations);
87      resetEvaluations();
88      final boolean forward = (t > t0);
89  
90      // create some internal working arrays
91      final int stages = c.length + 1;
92      if (y != y0) {
93        System.arraycopy(y0, 0, y, 0, y0.length);
94      }
95      final double[][] yDotK = new double[stages][];
96      for (int i = 0; i < stages; ++i) {
97        yDotK [i] = new double[y0.length];
98      }
99      final double[] yTmp = new double[y0.length];
100 
101     // set up an interpolator sharing the integrator arrays
102     AbstractStepInterpolator interpolator;
103     if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
104       final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
105       rki.reinitialize(this, yTmp, yDotK, forward);
106       interpolator = rki;
107     } else {
108       interpolator = new DummyStepInterpolator(yTmp, forward);
109     }
110     interpolator.storeTime(t0);
111 
112     // set up integration control objects
113     stepStart = t0;
114     stepSize  = forward ? step : -step;
115     for (StepHandler handler : stepHandlers) {
116         handler.reset();
117     }
118     CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
119     boolean lastStep = false;
120 
121     // main integration loop
122     while (!lastStep) {
123 
124       interpolator.shift();
125 
126       for (boolean loop = true; loop;) {
127 
128         // first stage
129         computeDerivatives(stepStart, y, yDotK[0]);
130 
131         // next stages
132         for (int k = 1; k < stages; ++k) {
133 
134           for (int j = 0; j < y0.length; ++j) {
135             double sum = a[k-1][0] * yDotK[0][j];
136             for (int l = 1; l < k; ++l) {
137               sum += a[k-1][l] * yDotK[l][j];
138             }
139             yTmp[j] = y[j] + stepSize * sum;
140           }
141 
142           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
143 
144         }
145 
146         // estimate the state at the end of the step
147         for (int j = 0; j < y0.length; ++j) {
148           double sum    = b[0] * yDotK[0][j];
149           for (int l = 1; l < stages; ++l) {
150             sum    += b[l] * yDotK[l][j];
151           }
152           yTmp[j] = y[j] + stepSize * sum;
153         }
154 
155         // discrete events handling
156         interpolator.storeTime(stepStart + stepSize);
157         if (manager.evaluateStep(interpolator)) {
158             final double dt = manager.getEventTime() - stepStart;
159             if (Math.abs(dt) <= Math.ulp(stepStart)) {
160                 // rejecting the step would lead to a too small next step, we accept it
161                 loop = false;
162             } else {
163                 // reject the step to match exactly the next switch time
164                 stepSize = dt;
165             }
166         } else {
167           loop = false;
168         }
169 
170       }
171 
172       // the step has been accepted
173       final double nextStep = stepStart + stepSize;
174       System.arraycopy(yTmp, 0, y, 0, y0.length);
175       manager.stepAccepted(nextStep, y);
176       lastStep = manager.stop();
177 
178       // provide the step data to the step handler
179       interpolator.storeTime(nextStep);
180       for (StepHandler handler : stepHandlers) {
181           handler.handleStep(interpolator, lastStep);
182       }
183       stepStart = nextStep;
184 
185       if (manager.reset(stepStart, y) && ! lastStep) {
186         // some events handler has triggered changes that
187         // invalidate the derivatives, we need to recompute them
188         computeDerivatives(stepStart, y, yDotK[0]);
189       }
190 
191       // make sure step size is set to default before next step
192       stepSize = forward ? step : -step;
193 
194     }
195 
196     final double stopTime = stepStart;
197     stepStart = Double.NaN;
198     stepSize  = Double.NaN;
199     return stopTime;
200 
201   }
202 
203   /** Time steps from Butcher array (without the first zero). */
204   private double[] c;
205 
206   /** Internal weights from Butcher array (without the first empty row). */
207   private double[][] a;
208 
209   /** External weights for the high order method from Butcher array. */
210   private double[] b;
211 
212   /** Prototype of the step interpolator. */
213   private RungeKuttaStepInterpolator prototype;
214                                          
215   /** Integration step. */
216   private double step;
217 
218 }