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