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;
19
20 import java.util.Collection;
21
22 /**
23 * This class implements the common part of all fixed step Runge-Kutta
24 * integrators for Ordinary Differential Equations.
25 *
26 * <p>These methods are explicit Runge-Kutta methods, their Butcher
27 * arrays are as follows :
28 * <pre>
29 * 0 |
30 * c2 | a21
31 * c3 | a31 a32
32 * ... | ...
33 * cs | as1 as2 ... ass-1
34 * |--------------------------
35 * | b1 b2 ... bs-1 bs
36 * </pre>
37 * </p>
38 *
39 * @see EulerIntegrator
40 * @see ClassicalRungeKuttaIntegrator
41 * @see GillIntegrator
42 * @see MidpointIntegrator
43 * @version $Revision: 666292 $ $Date: 2008-06-10 21:32:52 +0200 (mar., 10 juin 2008) $
44 * @since 1.2
45 */
46
47 public abstract class RungeKuttaIntegrator
48 implements FirstOrderIntegrator {
49
50 /** Simple constructor.
51 * Build a Runge-Kutta integrator with the given
52 * step. The default step handler does nothing.
53 * @param c time steps from Butcher array (without the first zero)
54 * @param a internal weights from Butcher array (without the first empty row)
55 * @param b propagation weights for the high order method from Butcher array
56 * @param prototype prototype of the step interpolator to use
57 * @param step integration step
58 */
59 protected RungeKuttaIntegrator(double[] c, double[][] a, double[] b,
60 RungeKuttaStepInterpolator prototype,
61 double step) {
62 this.c = c;
63 this.a = a;
64 this.b = b;
65 this.prototype = prototype;
66 this.step = step;
67 handler = DummyStepHandler.getInstance();
68 switchesHandler = new SwitchingFunctionsHandler();
69 resetInternalState();
70 }
71
72 /** Get the name of the method.
73 * @return name of the method
74 */
75 public abstract String getName();
76
77 /** Set the step handler for this integrator.
78 * The handler will be called by the integrator for each accepted
79 * step.
80 * @param handler handler for the accepted steps
81 */
82 public void setStepHandler (StepHandler handler) {
83 this.handler = handler;
84 }
85
86 /** Get the step handler for this integrator.
87 * @return the step handler for this integrator
88 */
89 public StepHandler getStepHandler() {
90 return handler;
91 }
92
93 /** Add a switching function to the integrator.
94 * @param function switching function
95 * @param maxCheckInterval maximal time interval between switching
96 * function checks (this interval prevents missing sign changes in
97 * case the integration steps becomes very large)
98 * @param convergence convergence threshold in the event time search
99 * @param maxIterationCount upper limit of the iteration count in
100 * the event time search
101 * @see #getSwitchingFunctions()
102 * @see #clearSwitchingFunctions()
103 */
104 public void addSwitchingFunction(SwitchingFunction function,
105 double maxCheckInterval,
106 double convergence,
107 int maxIterationCount) {
108 switchesHandler.addSwitchingFunction(function, maxCheckInterval, convergence, maxIterationCount);
109 }
110
111 /** Get all the switching functions that have been added to the integrator.
112 * @return an unmodifiable collection of the added switching functions
113 * @see #addSwitchingFunction(SwitchingFunction, double, double, int)
114 * @see #clearSwitchingFunctions()
115 */
116 public Collection<SwitchState> getSwitchingFunctions() {
117 return switchesHandler.getSwitchingFunctions();
118 }
119
120 /** Remove all the switching functions that have been added to the integrator.
121 * @see #addSwitchingFunction(SwitchingFunction, double, double, int)
122 * @see #getSwitchingFunctions()
123 */
124 public void clearSwitchingFunctions() {
125 switchesHandler.clearSwitchingFunctions();
126 }
127
128 /** Perform some sanity checks on the integration parameters.
129 * @param equations differential equations set
130 * @param t0 start time
131 * @param y0 state vector at t0
132 * @param t target time for the integration
133 * @param y placeholder where to put the state vector
134 * @exception IntegratorException if some inconsistency is detected
135 */
136 private void sanityChecks(FirstOrderDifferentialEquations equations,
137 double t0, double[] y0, double t, double[] y)
138 throws IntegratorException {
139 if (equations.getDimension() != y0.length) {
140 throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0}," +
141 " initial state vector has dimension {1}",
142 new Object[] {
143 Integer.valueOf(equations.getDimension()),
144 Integer.valueOf(y0.length)
145 });
146 }
147 if (equations.getDimension() != y.length) {
148 throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0}," +
149 " final state vector has dimension {1}",
150 new Object[] {
151 Integer.valueOf(equations.getDimension()),
152 Integer.valueOf(y.length)
153 });
154 }
155 if (Math.abs(t - t0) <= 1.0e-12 * Math.max(Math.abs(t0), Math.abs(t))) {
156 throw new IntegratorException("too small integration interval: length = {0}",
157 new Object[] { Double.valueOf(Math.abs(t - t0)) });
158 }
159 }
160
161 /** Integrate the differential equations up to the given time.
162 * <p>This method solves an Initial Value Problem (IVP).</p>
163 * <p>Since this method stores some internal state variables made
164 * available in its public interface during integration ({@link
165 * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
166 * @param equations differential equations to integrate
167 * @param t0 initial time
168 * @param y0 initial value of the state vector at t0
169 * @param t target time for the integration
170 * (can be set to a value smaller than <code>t0</code> for backward integration)
171 * @param y placeholder where to put the state vector at each successful
172 * step (and hence at the end of integration), can be the same object as y0
173 * @throws IntegratorException if the integrator cannot perform integration
174 * @throws DerivativeException this exception is propagated to the caller if
175 * the underlying user function triggers one
176 */
177 public void integrate(FirstOrderDifferentialEquations equations,
178 double t0, double[] y0,
179 double t, double[] y)
180 throws DerivativeException, IntegratorException {
181
182 sanityChecks(equations, t0, y0, t, y);
183 boolean forward = (t > t0);
184
185 // create some internal working arrays
186 int stages = c.length + 1;
187 if (y != y0) {
188 System.arraycopy(y0, 0, y, 0, y0.length);
189 }
190 double[][] yDotK = new double[stages][];
191 for (int i = 0; i < stages; ++i) {
192 yDotK [i] = new double[y0.length];
193 }
194 double[] yTmp = new double[y0.length];
195
196 // set up an interpolator sharing the integrator arrays
197 AbstractStepInterpolator interpolator;
198 if (handler.requiresDenseOutput() || (! switchesHandler.isEmpty())) {
199 RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
200 rki.reinitialize(equations, yTmp, yDotK, forward);
201 interpolator = rki;
202 } else {
203 interpolator = new DummyStepInterpolator(yTmp, forward);
204 }
205 interpolator.storeTime(t0);
206
207 // recompute the step
208 long nbStep = Math.max(1l, Math.abs(Math.round((t - t0) / step)));
209 boolean lastStep = false;
210 stepStart = t0;
211 stepSize = (t - t0) / nbStep;
212 handler.reset();
213 for (long i = 0; ! lastStep; ++i) {
214
215 interpolator.shift();
216
217 boolean needUpdate = false;
218 for (boolean loop = true; loop;) {
219
220 // first stage
221 equations.computeDerivatives(stepStart, y, yDotK[0]);
222
223 // next stages
224 for (int k = 1; k < stages; ++k) {
225
226 for (int j = 0; j < y0.length; ++j) {
227 double sum = a[k-1][0] * yDotK[0][j];
228 for (int l = 1; l < k; ++l) {
229 sum += a[k-1][l] * yDotK[l][j];
230 }
231 yTmp[j] = y[j] + stepSize * sum;
232 }
233
234 equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
235
236 }
237
238 // estimate the state at the end of the step
239 for (int j = 0; j < y0.length; ++j) {
240 double sum = b[0] * yDotK[0][j];
241 for (int l = 1; l < stages; ++l) {
242 sum += b[l] * yDotK[l][j];
243 }
244 yTmp[j] = y[j] + stepSize * sum;
245 }
246
247 // Switching functions handling
248 interpolator.storeTime(stepStart + stepSize);
249 if (switchesHandler.evaluateStep(interpolator)) {
250 needUpdate = true;
251 stepSize = switchesHandler.getEventTime() - stepStart;
252 } else {
253 loop = false;
254 }
255
256 }
257
258 // the step has been accepted
259 double nextStep = stepStart + stepSize;
260 System.arraycopy(yTmp, 0, y, 0, y0.length);
261 switchesHandler.stepAccepted(nextStep, y);
262 if (switchesHandler.stop()) {
263 lastStep = true;
264 } else {
265 lastStep = (i == (nbStep - 1));
266 }
267
268 // provide the step data to the step handler
269 interpolator.storeTime(nextStep);
270 handler.handleStep(interpolator, lastStep);
271 stepStart = nextStep;
272
273 if (switchesHandler.reset(stepStart, y) && ! lastStep) {
274 // some switching function has triggered changes that
275 // invalidate the derivatives, we need to recompute them
276 equations.computeDerivatives(stepStart, y, yDotK[0]);
277 }
278
279 if (needUpdate) {
280 // a switching function has changed the step
281 // we need to recompute stepsize
282 nbStep = Math.max(1l, Math.abs(Math.round((t - stepStart) / step)));
283 stepSize = (t - stepStart) / nbStep;
284 i = -1;
285 }
286
287 }
288
289 resetInternalState();
290
291 }
292
293 /** Get the current value of the step start time t<sub>i</sub>.
294 * <p>This method can be called during integration (typically by
295 * the object implementing the {@link FirstOrderDifferentialEquations
296 * differential equations} problem) if the value of the current step that
297 * is attempted is needed.</p>
298 * <p>The result is undefined if the method is called outside of
299 * calls to {@link #integrate}</p>
300 * @return current value of the step start time t<sub>i</sub>
301 */
302 public double getCurrentStepStart() {
303 return stepStart;
304 }
305
306 /** Get the current signed value of the integration stepsize.
307 * <p>This method can be called during integration (typically by
308 * the object implementing the {@link FirstOrderDifferentialEquations
309 * differential equations} problem) if the signed value of the current stepsize
310 * that is tried is needed.</p>
311 * <p>The result is undefined if the method is called outside of
312 * calls to {@link #integrate}</p>
313 * @return current signed value of the stepsize
314 */
315 public double getCurrentSignedStepsize() {
316 return stepSize;
317 }
318
319 /** Reset internal state to dummy values. */
320 private void resetInternalState() {
321 stepStart = Double.NaN;
322 stepSize = Double.NaN;
323 }
324
325 /** Time steps from Butcher array (without the first zero). */
326 private double[] c;
327
328 /** Internal weights from Butcher array (without the first empty row). */
329 private double[][] a;
330
331 /** External weights for the high order method from Butcher array. */
332 private double[] b;
333
334 /** Prototype of the step interpolator. */
335 private RungeKuttaStepInterpolator prototype;
336
337 /** Integration step. */
338 private double step;
339
340 /** Step handler. */
341 private StepHandler handler;
342
343 /** Switching functions handler. */
344 protected SwitchingFunctionsHandler switchesHandler;
345
346 /** Current step start time. */
347 private double stepStart;
348
349 /** Current stepsize. */
350 private double stepSize;
351
352 }