001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.ode.nonstiff;
019
020 import org.apache.commons.math.ode.DerivativeException;
021 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
022 import org.apache.commons.math.ode.IntegratorException;
023 import org.apache.commons.math.ode.events.CombinedEventsManager;
024 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
025 import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
026 import org.apache.commons.math.ode.sampling.StepHandler;
027
028 /**
029 * This class implements the common part of all embedded Runge-Kutta
030 * integrators for Ordinary Differential Equations.
031 *
032 * <p>These methods are embedded explicit Runge-Kutta methods with two
033 * sets of coefficients allowing to estimate the error, their Butcher
034 * arrays are as follows :
035 * <pre>
036 * 0 |
037 * c2 | a21
038 * c3 | a31 a32
039 * ... | ...
040 * cs | as1 as2 ... ass-1
041 * |--------------------------
042 * | b1 b2 ... bs-1 bs
043 * | b'1 b'2 ... b's-1 b's
044 * </pre>
045 * </p>
046 *
047 * <p>In fact, we rather use the array defined by ej = bj - b'j to
048 * compute directly the error rather than computing two estimates and
049 * then comparing them.</p>
050 *
051 * <p>Some methods are qualified as <i>fsal</i> (first same as last)
052 * methods. This means the last evaluation of the derivatives in one
053 * step is the same as the first in the next step. Then, this
054 * evaluation can be reused from one step to the next one and the cost
055 * of such a method is really s-1 evaluations despite the method still
056 * has s stages. This behaviour is true only for successful steps, if
057 * the step is rejected after the error estimation phase, no
058 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
059 * asi = bi for all i.</p>
060 *
061 * @version $Revision: 927202 $ $Date: 2010-03-24 18:11:51 -0400 (Wed, 24 Mar 2010) $
062 * @since 1.2
063 */
064
065 public abstract class EmbeddedRungeKuttaIntegrator
066 extends AdaptiveStepsizeIntegrator {
067
068 /** Indicator for <i>fsal</i> methods. */
069 private final boolean fsal;
070
071 /** Time steps from Butcher array (without the first zero). */
072 private final double[] c;
073
074 /** Internal weights from Butcher array (without the first empty row). */
075 private final double[][] a;
076
077 /** External weights for the high order method from Butcher array. */
078 private final double[] b;
079
080 /** Prototype of the step interpolator. */
081 private final RungeKuttaStepInterpolator prototype;
082
083 /** Stepsize control exponent. */
084 private final double exp;
085
086 /** Safety factor for stepsize control. */
087 private double safety;
088
089 /** Minimal reduction factor for stepsize control. */
090 private double minReduction;
091
092 /** Maximal growth factor for stepsize control. */
093 private double maxGrowth;
094
095 /** Build a Runge-Kutta integrator with the given Butcher array.
096 * @param name name of the method
097 * @param fsal indicate that the method is an <i>fsal</i>
098 * @param c time steps from Butcher array (without the first zero)
099 * @param a internal weights from Butcher array (without the first empty row)
100 * @param b propagation weights for the high order method from Butcher array
101 * @param prototype prototype of the step interpolator to use
102 * @param minStep minimal step (must be positive even for backward
103 * integration), the last step can be smaller than this
104 * @param maxStep maximal step (must be positive even for backward
105 * integration)
106 * @param scalAbsoluteTolerance allowed absolute error
107 * @param scalRelativeTolerance allowed relative error
108 */
109 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
110 final double[] c, final double[][] a, final double[] b,
111 final RungeKuttaStepInterpolator prototype,
112 final double minStep, final double maxStep,
113 final double scalAbsoluteTolerance,
114 final double scalRelativeTolerance) {
115
116 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
117
118 this.fsal = fsal;
119 this.c = c;
120 this.a = a;
121 this.b = b;
122 this.prototype = prototype;
123
124 exp = -1.0 / getOrder();
125
126 // set the default values of the algorithm control parameters
127 setSafety(0.9);
128 setMinReduction(0.2);
129 setMaxGrowth(10.0);
130
131 }
132
133 /** Build a Runge-Kutta integrator with the given Butcher array.
134 * @param name name of the method
135 * @param fsal indicate that the method is an <i>fsal</i>
136 * @param c time steps from Butcher array (without the first zero)
137 * @param a internal weights from Butcher array (without the first empty row)
138 * @param b propagation weights for the high order method from Butcher array
139 * @param prototype prototype of the step interpolator to use
140 * @param minStep minimal step (must be positive even for backward
141 * integration), the last step can be smaller than this
142 * @param maxStep maximal step (must be positive even for backward
143 * integration)
144 * @param vecAbsoluteTolerance allowed absolute error
145 * @param vecRelativeTolerance allowed relative error
146 */
147 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
148 final double[] c, final double[][] a, final double[] b,
149 final RungeKuttaStepInterpolator prototype,
150 final double minStep, final double maxStep,
151 final double[] vecAbsoluteTolerance,
152 final double[] vecRelativeTolerance) {
153
154 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
155
156 this.fsal = fsal;
157 this.c = c;
158 this.a = a;
159 this.b = b;
160 this.prototype = prototype;
161
162 exp = -1.0 / getOrder();
163
164 // set the default values of the algorithm control parameters
165 setSafety(0.9);
166 setMinReduction(0.2);
167 setMaxGrowth(10.0);
168
169 }
170
171 /** Get the order of the method.
172 * @return order of the method
173 */
174 public abstract int getOrder();
175
176 /** Get the safety factor for stepsize control.
177 * @return safety factor
178 */
179 public double getSafety() {
180 return safety;
181 }
182
183 /** Set the safety factor for stepsize control.
184 * @param safety safety factor
185 */
186 public void setSafety(final double safety) {
187 this.safety = safety;
188 }
189
190 /** {@inheritDoc} */
191 @Override
192 public double integrate(final FirstOrderDifferentialEquations equations,
193 final double t0, final double[] y0,
194 final double t, final double[] y)
195 throws DerivativeException, IntegratorException {
196
197 sanityChecks(equations, t0, y0, t, y);
198 setEquations(equations);
199 resetEvaluations();
200 final boolean forward = t > t0;
201
202 // create some internal working arrays
203 final int stages = c.length + 1;
204 if (y != y0) {
205 System.arraycopy(y0, 0, y, 0, y0.length);
206 }
207 final double[][] yDotK = new double[stages][y0.length];
208 final double[] yTmp = new double[y0.length];
209
210 // set up an interpolator sharing the integrator arrays
211 AbstractStepInterpolator interpolator;
212 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
213 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
214 rki.reinitialize(this, yTmp, yDotK, forward);
215 interpolator = rki;
216 } else {
217 interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
218 }
219 interpolator.storeTime(t0);
220
221 // set up integration control objects
222 stepStart = t0;
223 double hNew = 0;
224 boolean firstTime = true;
225 for (StepHandler handler : stepHandlers) {
226 handler.reset();
227 }
228 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
229 boolean lastStep = false;
230
231 // main integration loop
232 while (!lastStep) {
233
234 interpolator.shift();
235
236 double error = 0;
237 for (boolean loop = true; loop;) {
238
239 if (firstTime || !fsal) {
240 // first stage
241 computeDerivatives(stepStart, y, yDotK[0]);
242 }
243
244 if (firstTime) {
245 final double[] scale = new double[y0.length];
246 if (vecAbsoluteTolerance == null) {
247 for (int i = 0; i < scale.length; ++i) {
248 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * Math.abs(y[i]);
249 }
250 } else {
251 for (int i = 0; i < scale.length; ++i) {
252 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * Math.abs(y[i]);
253 }
254 }
255 hNew = initializeStep(equations, forward, getOrder(), scale,
256 stepStart, y, yDotK[0], yTmp, yDotK[1]);
257 firstTime = false;
258 }
259
260 stepSize = hNew;
261
262 // next stages
263 for (int k = 1; k < stages; ++k) {
264
265 for (int j = 0; j < y0.length; ++j) {
266 double sum = a[k-1][0] * yDotK[0][j];
267 for (int l = 1; l < k; ++l) {
268 sum += a[k-1][l] * yDotK[l][j];
269 }
270 yTmp[j] = y[j] + stepSize * sum;
271 }
272
273 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
274
275 }
276
277 // estimate the state at the end of the step
278 for (int j = 0; j < y0.length; ++j) {
279 double sum = b[0] * yDotK[0][j];
280 for (int l = 1; l < stages; ++l) {
281 sum += b[l] * yDotK[l][j];
282 }
283 yTmp[j] = y[j] + stepSize * sum;
284 }
285
286 // estimate the error at the end of the step
287 error = estimateError(yDotK, y, yTmp, stepSize);
288 if (error <= 1.0) {
289
290 // discrete events handling
291 interpolator.storeTime(stepStart + stepSize);
292 if (manager.evaluateStep(interpolator)) {
293 final double dt = manager.getEventTime() - stepStart;
294 if (Math.abs(dt) <= Math.ulp(stepStart)) {
295 // we cannot simply truncate the step, reject the current computation
296 // and let the loop compute another state with the truncated step.
297 // it is so small (much probably exactly 0 due to limited accuracy)
298 // that the code above would fail handling it.
299 // So we set up an artificial 0 size step by copying states
300 interpolator.storeTime(stepStart);
301 System.arraycopy(y, 0, yTmp, 0, y0.length);
302 hNew = 0;
303 stepSize = 0;
304 loop = false;
305 } else {
306 // reject the step to match exactly the next switch time
307 hNew = dt;
308 }
309 } else {
310 // accept the step
311 loop = false;
312 }
313
314 } else {
315 // reject the step and attempt to reduce error by stepsize control
316 final double factor =
317 Math.min(maxGrowth,
318 Math.max(minReduction, safety * Math.pow(error, exp)));
319 hNew = filterStep(stepSize * factor, forward, false);
320 }
321
322 }
323
324 // the step has been accepted
325 final double nextStep = stepStart + stepSize;
326 System.arraycopy(yTmp, 0, y, 0, y0.length);
327 manager.stepAccepted(nextStep, y);
328 lastStep = manager.stop();
329
330 // provide the step data to the step handler
331 interpolator.storeTime(nextStep);
332 for (StepHandler handler : stepHandlers) {
333 handler.handleStep(interpolator, lastStep);
334 }
335 stepStart = nextStep;
336
337 if (fsal) {
338 // save the last evaluation for the next step
339 System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
340 }
341
342 if (manager.reset(stepStart, y) && ! lastStep) {
343 // some event handler has triggered changes that
344 // invalidate the derivatives, we need to recompute them
345 computeDerivatives(stepStart, y, yDotK[0]);
346 }
347
348 if (! lastStep) {
349 // in some rare cases we may get here with stepSize = 0, for example
350 // when an event occurs at integration start, reducing the first step
351 // to zero; we have to reset the step to some safe non zero value
352 stepSize = filterStep(stepSize, forward, true);
353
354 // stepsize control for next step
355 final double factor = Math.min(maxGrowth,
356 Math.max(minReduction,
357 safety * Math.pow(error, exp)));
358 final double scaledH = stepSize * factor;
359 final double nextT = stepStart + scaledH;
360 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
361 hNew = filterStep(scaledH, forward, nextIsLast);
362 }
363
364 }
365
366 final double stopTime = stepStart;
367 resetInternalState();
368 return stopTime;
369
370 }
371
372 /** Get the minimal reduction factor for stepsize control.
373 * @return minimal reduction factor
374 */
375 public double getMinReduction() {
376 return minReduction;
377 }
378
379 /** Set the minimal reduction factor for stepsize control.
380 * @param minReduction minimal reduction factor
381 */
382 public void setMinReduction(final double minReduction) {
383 this.minReduction = minReduction;
384 }
385
386 /** Get the maximal growth factor for stepsize control.
387 * @return maximal growth factor
388 */
389 public double getMaxGrowth() {
390 return maxGrowth;
391 }
392
393 /** Set the maximal growth factor for stepsize control.
394 * @param maxGrowth maximal growth factor
395 */
396 public void setMaxGrowth(final double maxGrowth) {
397 this.maxGrowth = maxGrowth;
398 }
399
400 /** Compute the error ratio.
401 * @param yDotK derivatives computed during the first stages
402 * @param y0 estimate of the step at the start of the step
403 * @param y1 estimate of the step at the end of the step
404 * @param h current step
405 * @return error ratio, greater than 1 if step should be rejected
406 */
407 protected abstract double estimateError(double[][] yDotK,
408 double[] y0, double[] y1,
409 double h);
410
411 }