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 import java.util.Arrays;
21
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 embedded Runge-Kutta
32 * integrators for Ordinary Differential Equations.
33 *
34 * <p>These methods are embedded explicit Runge-Kutta methods with two
35 * sets of coefficients allowing to estimate the error, their Butcher
36 * arrays are as follows :
37 * <pre>
38 * 0 |
39 * c2 | a21
40 * c3 | a31 a32
41 * ... | ...
42 * cs | as1 as2 ... ass-1
43 * |--------------------------
44 * | b1 b2 ... bs-1 bs
45 * | b'1 b'2 ... b's-1 b's
46 * </pre>
47 * </p>
48 *
49 * <p>In fact, we rather use the array defined by ej = bj - b'j to
50 * compute directly the error rather than computing two estimates and
51 * then comparing them.</p>
52 *
53 * <p>Some methods are qualified as <i>fsal</i> (first same as last)
54 * methods. This means the last evaluation of the derivatives in one
55 * step is the same as the first in the next step. Then, this
56 * evaluation can be reused from one step to the next one and the cost
57 * of such a method is really s-1 evaluations despite the method still
58 * has s stages. This behaviour is true only for successful steps, if
59 * the step is rejected after the error estimation phase, no
60 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
61 * asi = bi for all i.</p>
62 *
63 * @version $Revision: 811833 $ $Date: 2009-09-06 12:27:50 -0400 (Sun, 06 Sep 2009) $
64 * @since 1.2
65 */
66
67 public abstract class EmbeddedRungeKuttaIntegrator
68 extends AdaptiveStepsizeIntegrator {
69
70 /** Indicator for <i>fsal</i> methods. */
71 private final boolean fsal;
72
73 /** Time steps from Butcher array (without the first zero). */
74 private final double[] c;
75
76 /** Internal weights from Butcher array (without the first empty row). */
77 private final double[][] a;
78
79 /** External weights for the high order method from Butcher array. */
80 private final double[] b;
81
82 /** Prototype of the step interpolator. */
83 private final RungeKuttaStepInterpolator prototype;
84
85 /** Stepsize control exponent. */
86 private final double exp;
87
88 /** Safety factor for stepsize control. */
89 private double safety;
90
91 /** Minimal reduction factor for stepsize control. */
92 private double minReduction;
93
94 /** Maximal growth factor for stepsize control. */
95 private double maxGrowth;
96
97 /** Build a Runge-Kutta integrator with the given Butcher array.
98 * @param name name of the method
99 * @param fsal indicate that the method is an <i>fsal</i>
100 * @param c time steps from Butcher array (without the first zero)
101 * @param a internal weights from Butcher array (without the first empty row)
102 * @param b propagation weights for the high order method from Butcher array
103 * @param prototype prototype of the step interpolator to use
104 * @param minStep minimal step (must be positive even for backward
105 * integration), the last step can be smaller than this
106 * @param maxStep maximal step (must be positive even for backward
107 * integration)
108 * @param scalAbsoluteTolerance allowed absolute error
109 * @param scalRelativeTolerance allowed relative error
110 */
111 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
112 final double[] c, final double[][] a, final double[] b,
113 final RungeKuttaStepInterpolator prototype,
114 final double minStep, final double maxStep,
115 final double scalAbsoluteTolerance,
116 final double scalRelativeTolerance) {
117
118 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
119
120 this.fsal = fsal;
121 this.c = c;
122 this.a = a;
123 this.b = b;
124 this.prototype = prototype;
125
126 exp = -1.0 / getOrder();
127
128 // set the default values of the algorithm control parameters
129 setSafety(0.9);
130 setMinReduction(0.2);
131 setMaxGrowth(10.0);
132
133 }
134
135 /** Build a Runge-Kutta integrator with the given Butcher array.
136 * @param name name of the method
137 * @param fsal indicate that the method is an <i>fsal</i>
138 * @param c time steps from Butcher array (without the first zero)
139 * @param a internal weights from Butcher array (without the first empty row)
140 * @param b propagation weights for the high order method from Butcher array
141 * @param prototype prototype of the step interpolator to use
142 * @param minStep minimal step (must be positive even for backward
143 * integration), the last step can be smaller than this
144 * @param maxStep maximal step (must be positive even for backward
145 * integration)
146 * @param vecAbsoluteTolerance allowed absolute error
147 * @param vecRelativeTolerance allowed relative error
148 */
149 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
150 final double[] c, final double[][] a, final double[] b,
151 final RungeKuttaStepInterpolator prototype,
152 final double minStep, final double maxStep,
153 final double[] vecAbsoluteTolerance,
154 final double[] vecRelativeTolerance) {
155
156 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
157
158 this.fsal = fsal;
159 this.c = c;
160 this.a = a;
161 this.b = b;
162 this.prototype = prototype;
163
164 exp = -1.0 / getOrder();
165
166 // set the default values of the algorithm control parameters
167 setSafety(0.9);
168 setMinReduction(0.2);
169 setMaxGrowth(10.0);
170
171 }
172
173 /** Get the order of the method.
174 * @return order of the method
175 */
176 public abstract int getOrder();
177
178 /** Get the safety factor for stepsize control.
179 * @return safety factor
180 */
181 public double getSafety() {
182 return safety;
183 }
184
185 /** Set the safety factor for stepsize control.
186 * @param safety safety factor
187 */
188 public void setSafety(final double safety) {
189 this.safety = safety;
190 }
191
192 /** {@inheritDoc} */
193 @Override
194 public double integrate(final FirstOrderDifferentialEquations equations,
195 final double t0, final double[] y0,
196 final double t, final double[] y)
197 throws DerivativeException, IntegratorException {
198
199 sanityChecks(equations, t0, y0, t, y);
200 setEquations(equations);
201 resetEvaluations();
202 final boolean forward = t > t0;
203
204 // create some internal working arrays
205 final int stages = c.length + 1;
206 if (y != y0) {
207 System.arraycopy(y0, 0, y, 0, y0.length);
208 }
209 final double[][] yDotK = new double[stages][y0.length];
210 final double[] yTmp = new double[y0.length];
211
212 // set up an interpolator sharing the integrator arrays
213 AbstractStepInterpolator interpolator;
214 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
215 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
216 rki.reinitialize(this, yTmp, yDotK, forward);
217 interpolator = rki;
218 } else {
219 interpolator = new DummyStepInterpolator(yTmp, forward);
220 }
221 interpolator.storeTime(t0);
222
223 // set up integration control objects
224 stepStart = t0;
225 double hNew = 0;
226 boolean firstTime = true;
227 for (StepHandler handler : stepHandlers) {
228 handler.reset();
229 }
230 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
231 boolean lastStep = false;
232
233 // main integration loop
234 while (!lastStep) {
235
236 interpolator.shift();
237
238 double error = 0;
239 for (boolean loop = true; loop;) {
240
241 if (firstTime || !fsal) {
242 // first stage
243 computeDerivatives(stepStart, y, yDotK[0]);
244 }
245
246 if (firstTime) {
247 final double[] scale;
248 if (vecAbsoluteTolerance != null) {
249 scale = vecAbsoluteTolerance;
250 } else {
251 scale = new double[y0.length];
252 Arrays.fill(scale, scalAbsoluteTolerance);
253 }
254 hNew = initializeStep(equations, forward, getOrder(), scale,
255 stepStart, y, yDotK[0], yTmp, yDotK[1]);
256 firstTime = false;
257 }
258
259 stepSize = hNew;
260
261 // next stages
262 for (int k = 1; k < stages; ++k) {
263
264 for (int j = 0; j < y0.length; ++j) {
265 double sum = a[k-1][0] * yDotK[0][j];
266 for (int l = 1; l < k; ++l) {
267 sum += a[k-1][l] * yDotK[l][j];
268 }
269 yTmp[j] = y[j] + stepSize * sum;
270 }
271
272 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
273
274 }
275
276 // estimate the state at the end of the step
277 for (int j = 0; j < y0.length; ++j) {
278 double sum = b[0] * yDotK[0][j];
279 for (int l = 1; l < stages; ++l) {
280 sum += b[l] * yDotK[l][j];
281 }
282 yTmp[j] = y[j] + stepSize * sum;
283 }
284
285 // estimate the error at the end of the step
286 error = estimateError(yDotK, y, yTmp, stepSize);
287 if (error <= 1.0) {
288
289 // discrete events handling
290 interpolator.storeTime(stepStart + stepSize);
291 if (manager.evaluateStep(interpolator)) {
292 final double dt = manager.getEventTime() - stepStart;
293 if (Math.abs(dt) <= Math.ulp(stepStart)) {
294 // rejecting the step would lead to a too small next step, we accept it
295 loop = false;
296 } else {
297 // reject the step to match exactly the next switch time
298 hNew = dt;
299 }
300 } else {
301 // accept the step
302 loop = false;
303 }
304
305 } else {
306 // reject the step and attempt to reduce error by stepsize control
307 final double factor =
308 Math.min(maxGrowth,
309 Math.max(minReduction, safety * Math.pow(error, exp)));
310 hNew = filterStep(stepSize * factor, forward, false);
311 }
312
313 }
314
315 // the step has been accepted
316 final double nextStep = stepStart + stepSize;
317 System.arraycopy(yTmp, 0, y, 0, y0.length);
318 manager.stepAccepted(nextStep, y);
319 lastStep = manager.stop();
320
321 // provide the step data to the step handler
322 interpolator.storeTime(nextStep);
323 for (StepHandler handler : stepHandlers) {
324 handler.handleStep(interpolator, lastStep);
325 }
326 stepStart = nextStep;
327
328 if (fsal) {
329 // save the last evaluation for the next step
330 System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
331 }
332
333 if (manager.reset(stepStart, y) && ! lastStep) {
334 // some event handler has triggered changes that
335 // invalidate the derivatives, we need to recompute them
336 computeDerivatives(stepStart, y, yDotK[0]);
337 }
338
339 if (! lastStep) {
340 // in some rare cases we may get here with stepSize = 0, for example
341 // when an event occurs at integration start, reducing the first step
342 // to zero; we have to reset the step to some safe non zero value
343 stepSize = filterStep(stepSize, forward, true);
344
345 // stepsize control for next step
346 final double factor = Math.min(maxGrowth,
347 Math.max(minReduction,
348 safety * Math.pow(error, exp)));
349 final double scaledH = stepSize * factor;
350 final double nextT = stepStart + scaledH;
351 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
352 hNew = filterStep(scaledH, forward, nextIsLast);
353 }
354
355 }
356
357 final double stopTime = stepStart;
358 resetInternalState();
359 return stopTime;
360
361 }
362
363 /** Get the minimal reduction factor for stepsize control.
364 * @return minimal reduction factor
365 */
366 public double getMinReduction() {
367 return minReduction;
368 }
369
370 /** Set the minimal reduction factor for stepsize control.
371 * @param minReduction minimal reduction factor
372 */
373 public void setMinReduction(final double minReduction) {
374 this.minReduction = minReduction;
375 }
376
377 /** Get the maximal growth factor for stepsize control.
378 * @return maximal growth factor
379 */
380 public double getMaxGrowth() {
381 return maxGrowth;
382 }
383
384 /** Set the maximal growth factor for stepsize control.
385 * @param maxGrowth maximal growth factor
386 */
387 public void setMaxGrowth(final double maxGrowth) {
388 this.maxGrowth = maxGrowth;
389 }
390
391 /** Compute the error ratio.
392 * @param yDotK derivatives computed during the first stages
393 * @param y0 estimate of the step at the start of the step
394 * @param y1 estimate of the step at the end of the step
395 * @param h current step
396 * @return error ratio, greater than 1 if step should be rejected
397 */
398 protected abstract double estimateError(double[][] yDotK,
399 double[] y0, double[] y1,
400 double h);
401
402 }