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  
21  import org.apache.commons.math3.exception.DimensionMismatchException;
22  import org.apache.commons.math3.exception.MaxCountExceededException;
23  import org.apache.commons.math3.exception.NoBracketingException;
24  import org.apache.commons.math3.exception.NumberIsTooSmallException;
25  import org.apache.commons.math3.ode.AbstractIntegrator;
26  import org.apache.commons.math3.ode.ExpandableStatefulODE;
27  import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
28  import org.apache.commons.math3.util.FastMath;
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   * @since 1.2
52   */
53  
54  public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
55  
56      /** Time steps from Butcher array (without the first zero). */
57      private final double[] c;
58  
59      /** Internal weights from Butcher array (without the first empty row). */
60      private final double[][] a;
61  
62      /** External weights for the high order method from Butcher array. */
63      private final double[] b;
64  
65      /** Prototype of the step interpolator. */
66      private final RungeKuttaStepInterpolator prototype;
67  
68      /** Integration step. */
69      private final double step;
70  
71    /** Simple constructor.
72     * Build a Runge-Kutta integrator with the given
73     * step. The default step handler does nothing.
74     * @param name name of the method
75     * @param c time steps from Butcher array (without the first zero)
76     * @param a internal weights from Butcher array (without the first empty row)
77     * @param b propagation weights for the high order method from Butcher array
78     * @param prototype prototype of the step interpolator to use
79     * @param step integration step
80     */
81    protected RungeKuttaIntegrator(final String name,
82                                   final double[] c, final double[][] a, final double[] b,
83                                   final RungeKuttaStepInterpolator prototype,
84                                   final double step) {
85      super(name);
86      this.c          = c;
87      this.a          = a;
88      this.b          = b;
89      this.prototype  = prototype;
90      this.step       = FastMath.abs(step);
91    }
92  
93    /** {@inheritDoc} */
94    @Override
95    public void integrate(final ExpandableStatefulODE equations, final double t)
96        throws NumberIsTooSmallException, DimensionMismatchException,
97               MaxCountExceededException, NoBracketingException {
98  
99      sanityChecks(equations, t);
100     setEquations(equations);
101     final boolean forward = t > equations.getTime();
102 
103     // create some internal working arrays
104     final double[] y0      = equations.getCompleteState();
105     final double[] y       = y0.clone();
106     final int stages       = c.length + 1;
107     final double[][] yDotK = new double[stages][];
108     for (int i = 0; i < stages; ++i) {
109       yDotK [i] = new double[y0.length];
110     }
111     final double[] yTmp    = y0.clone();
112     final double[] yDotTmp = new double[y0.length];
113 
114     // set up an interpolator sharing the integrator arrays
115     final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
116     interpolator.reinitialize(this, yTmp, yDotK, forward,
117                               equations.getPrimaryMapper(), equations.getSecondaryMappers());
118     interpolator.storeTime(equations.getTime());
119 
120     // set up integration control objects
121     stepStart = equations.getTime();
122     if (forward) {
123         if (stepStart + step >= t) {
124             stepSize = t - stepStart;
125         } else {
126             stepSize = step;
127         }
128     } else {
129         if (stepStart - step <= t) {
130             stepSize = t - stepStart;
131         } else {
132             stepSize = -step;
133         }
134     }
135     initIntegration(equations.getTime(), y0, t);
136 
137     // main integration loop
138     isLastStep = false;
139     do {
140 
141       interpolator.shift();
142 
143       // first stage
144       computeDerivatives(stepStart, y, yDotK[0]);
145 
146       // next stages
147       for (int k = 1; k < stages; ++k) {
148 
149           for (int j = 0; j < y0.length; ++j) {
150               double sum = a[k-1][0] * yDotK[0][j];
151               for (int l = 1; l < k; ++l) {
152                   sum += a[k-1][l] * yDotK[l][j];
153               }
154               yTmp[j] = y[j] + stepSize * sum;
155           }
156 
157           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
158 
159       }
160 
161       // estimate the state at the end of the step
162       for (int j = 0; j < y0.length; ++j) {
163           double sum    = b[0] * yDotK[0][j];
164           for (int l = 1; l < stages; ++l) {
165               sum    += b[l] * yDotK[l][j];
166           }
167           yTmp[j] = y[j] + stepSize * sum;
168       }
169 
170       // discrete events handling
171       interpolator.storeTime(stepStart + stepSize);
172       System.arraycopy(yTmp, 0, y, 0, y0.length);
173       System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
174       stepStart = acceptStep(interpolator, y, yDotTmp, t);
175 
176       if (!isLastStep) {
177 
178           // prepare next step
179           interpolator.storeTime(stepStart);
180 
181           // stepsize control for next step
182           final double  nextT      = stepStart + stepSize;
183           final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
184           if (nextIsLast) {
185               stepSize = t - stepStart;
186           }
187       }
188 
189     } while (!isLastStep);
190 
191     // dispatch results
192     equations.setTime(stepStart);
193     equations.setCompleteState(y);
194 
195     stepStart = Double.NaN;
196     stepSize  = Double.NaN;
197 
198   }
199 
200   /** Fast computation of a single step of ODE integration.
201    * <p>This method is intended for the limited use case of
202    * very fast computation of only one step without using any of the
203    * rich features of general integrators that may take some time
204    * to set up (i.e. no step handlers, no events handlers, no additional
205    * states, no interpolators, no error control, no evaluations count,
206    * no sanity checks ...). It handles the strict minimum of computation,
207    * so it can be embedded in outer loops.</p>
208    * <p>
209    * This method is <em>not</em> used at all by the {@link #integrate(ExpandableStatefulODE, double)}
210    * method. It also completely ignores the step set at construction time, and
211    * uses only a single step to go from {@code t0} to {@code t}.
212    * </p>
213    * <p>
214    * As this method does not use any of the state-dependent features of the integrator,
215    * it should be reasonably thread-safe <em>if and only if</em> the provided differential
216    * equations are themselves thread-safe.
217    * </p>
218    * @param equations differential equations to integrate
219    * @param t0 initial time
220    * @param y0 initial value of the state vector at t0
221    * @param t target time for the integration
222    * (can be set to a value smaller than {@code t0} for backward integration)
223    * @return state vector at {@code t}
224    */
225   public double[] singleStep(final FirstOrderDifferentialEquations equations,
226                              final double t0, final double[] y0, final double t) {
227 
228       // create some internal working arrays
229       final double[] y       = y0.clone();
230       final int stages       = c.length + 1;
231       final double[][] yDotK = new double[stages][];
232       for (int i = 0; i < stages; ++i) {
233           yDotK [i] = new double[y0.length];
234       }
235       final double[] yTmp    = y0.clone();
236 
237       // first stage
238       final double h = t - t0;
239       equations.computeDerivatives(t0, y, yDotK[0]);
240 
241       // next stages
242       for (int k = 1; k < stages; ++k) {
243 
244           for (int j = 0; j < y0.length; ++j) {
245               double sum = a[k-1][0] * yDotK[0][j];
246               for (int l = 1; l < k; ++l) {
247                   sum += a[k-1][l] * yDotK[l][j];
248               }
249               yTmp[j] = y[j] + h * sum;
250           }
251 
252           equations.computeDerivatives(t0 + c[k-1] * h, yTmp, yDotK[k]);
253 
254       }
255 
256       // estimate the state at the end of the step
257       for (int j = 0; j < y0.length; ++j) {
258           double sum = b[0] * yDotK[0][j];
259           for (int l = 1; l < stages; ++l) {
260               sum += b[l] * yDotK[l][j];
261           }
262           y[j] += h * sum;
263       }
264 
265       return y;
266 
267   }
268 
269 }