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.util.FastMath;
28  
29  /**
30   * This class implements the common part of all fixed step Runge-Kutta
31   * integrators for Ordinary Differential Equations.
32   *
33   * <p>These methods are explicit Runge-Kutta methods, their Butcher
34   * arrays are as follows :
35   * <pre>
36   *    0  |
37   *   c2  | a21
38   *   c3  | a31  a32
39   *   ... |        ...
40   *   cs  | as1  as2  ...  ass-1
41   *       |--------------------------
42   *       |  b1   b2  ...   bs-1  bs
43   * </pre>
44   * </p>
45   *
46   * @see EulerIntegrator
47   * @see ClassicalRungeKuttaIntegrator
48   * @see GillIntegrator
49   * @see MidpointIntegrator
50   * @version $Id: RungeKuttaIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
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     stepSize  = forward ? step : -step;
123     initIntegration(equations.getTime(), y0, t);
124 
125     // main integration loop
126     isLastStep = false;
127     do {
128 
129       interpolator.shift();
130 
131       // first stage
132       computeDerivatives(stepStart, y, yDotK[0]);
133 
134       // next stages
135       for (int k = 1; k < stages; ++k) {
136 
137           for (int j = 0; j < y0.length; ++j) {
138               double sum = a[k-1][0] * yDotK[0][j];
139               for (int l = 1; l < k; ++l) {
140                   sum += a[k-1][l] * yDotK[l][j];
141               }
142               yTmp[j] = y[j] + stepSize * sum;
143           }
144 
145           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
146 
147       }
148 
149       // estimate the state at the end of the step
150       for (int j = 0; j < y0.length; ++j) {
151           double sum    = b[0] * yDotK[0][j];
152           for (int l = 1; l < stages; ++l) {
153               sum    += b[l] * yDotK[l][j];
154           }
155           yTmp[j] = y[j] + stepSize * sum;
156       }
157 
158       // discrete events handling
159       interpolator.storeTime(stepStart + stepSize);
160       System.arraycopy(yTmp, 0, y, 0, y0.length);
161       System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
162       stepStart = acceptStep(interpolator, y, yDotTmp, t);
163 
164       if (!isLastStep) {
165 
166           // prepare next step
167           interpolator.storeTime(stepStart);
168 
169           // stepsize control for next step
170           final double  nextT      = stepStart + stepSize;
171           final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
172           if (nextIsLast) {
173               stepSize = t - stepStart;
174           }
175       }
176 
177     } while (!isLastStep);
178 
179     // dispatch results
180     equations.setTime(stepStart);
181     equations.setCompleteState(y);
182 
183     stepStart = Double.NaN;
184     stepSize  = Double.NaN;
185 
186   }
187 
188 }