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 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
21 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
22 import org.apache.commons.math4.legacy.exception.NoBracketingException;
23 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
24 import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE;
25 import org.apache.commons.math4.core.jdkmath.JdkMath;
26
27 /**
28 * This class implements the common part of all embedded Runge-Kutta
29 * integrators for Ordinary Differential Equations.
30 *
31 * <p>These methods are embedded explicit Runge-Kutta methods with two
32 * sets of coefficients allowing to estimate the error, their Butcher
33 * arrays are as follows :
34 * <pre>
35 * 0 |
36 * c2 | a21
37 * c3 | a31 a32
38 * ... | ...
39 * cs | as1 as2 ... ass-1
40 * |--------------------------
41 * | b1 b2 ... bs-1 bs
42 * | b'1 b'2 ... b's-1 b's
43 * </pre>
44 *
45 * <p>In fact, we rather use the array defined by ej = bj - b'j to
46 * compute directly the error rather than computing two estimates and
47 * then comparing them.</p>
48 *
49 * <p>Some methods are qualified as <i>fsal</i> (first same as last)
50 * methods. This means the last evaluation of the derivatives in one
51 * step is the same as the first in the next step. Then, this
52 * evaluation can be reused from one step to the next one and the cost
53 * of such a method is really s-1 evaluations despite the method still
54 * has s stages. This behaviour is true only for successful steps, if
55 * the step is rejected after the error estimation phase, no
56 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
57 * asi = bi for all i.</p>
58 *
59 * @since 1.2
60 */
61
62 public abstract class EmbeddedRungeKuttaIntegrator
63 extends AdaptiveStepsizeIntegrator {
64
65 /** Indicator for <i>fsal</i> methods. */
66 private final boolean fsal;
67
68 /** Time steps from Butcher array (without the first zero). */
69 private final double[] c;
70
71 /** Internal weights from Butcher array (without the first empty row). */
72 private final double[][] a;
73
74 /** External weights for the high order method from Butcher array. */
75 private final double[] b;
76
77 /** Prototype of the step interpolator. */
78 private final RungeKuttaStepInterpolator prototype;
79
80 /** Stepsize control exponent. */
81 private final double exp;
82
83 /** Safety factor for stepsize control. */
84 private double safety;
85
86 /** Minimal reduction factor for stepsize control. */
87 private double minReduction;
88
89 /** Maximal growth factor for stepsize control. */
90 private double maxGrowth;
91
92 /** Build a Runge-Kutta integrator with the given Butcher array.
93 * @param name name of the method
94 * @param fsal indicate that the method is an <i>fsal</i>
95 * @param c time steps from Butcher array (without the first zero)
96 * @param a internal weights from Butcher array (without the first empty row)
97 * @param b propagation weights for the high order method from Butcher array
98 * @param prototype prototype of the step interpolator to use
99 * @param minStep minimal step (sign is irrelevant, regardless of
100 * integration direction, forward or backward), the last step can
101 * be smaller than this
102 * @param maxStep maximal step (sign is irrelevant, regardless of
103 * integration direction, forward or backward), the last step can
104 * be smaller than this
105 * @param scalAbsoluteTolerance allowed absolute error
106 * @param scalRelativeTolerance allowed relative error
107 */
108 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
109 final double[] c, final double[][] a, final double[] b,
110 final RungeKuttaStepInterpolator prototype,
111 final double minStep, final double maxStep,
112 final double scalAbsoluteTolerance,
113 final double scalRelativeTolerance) {
114
115 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
116
117 this.fsal = fsal;
118 this.c = c;
119 this.a = a;
120 this.b = b;
121 this.prototype = prototype;
122
123 exp = -1.0 / getOrder();
124
125 // set the default values of the algorithm control parameters
126 setSafety(0.9);
127 setMinReduction(0.2);
128 setMaxGrowth(10.0);
129 }
130
131 /** Build a Runge-Kutta integrator with the given Butcher array.
132 * @param name name of the method
133 * @param fsal indicate that the method is an <i>fsal</i>
134 * @param c time steps from Butcher array (without the first zero)
135 * @param a internal weights from Butcher array (without the first empty row)
136 * @param b propagation weights for the high order method from Butcher array
137 * @param prototype prototype of the step interpolator to use
138 * @param minStep minimal step (must be positive even for backward
139 * integration), the last step can be smaller than this
140 * @param maxStep maximal step (must be positive even for backward
141 * integration)
142 * @param vecAbsoluteTolerance allowed absolute error
143 * @param vecRelativeTolerance allowed relative error
144 */
145 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
146 final double[] c, final double[][] a, final double[] b,
147 final RungeKuttaStepInterpolator prototype,
148 final double minStep, final double maxStep,
149 final double[] vecAbsoluteTolerance,
150 final double[] vecRelativeTolerance) {
151
152 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
153
154 this.fsal = fsal;
155 this.c = c;
156 this.a = a;
157 this.b = b;
158 this.prototype = prototype;
159
160 exp = -1.0 / getOrder();
161
162 // set the default values of the algorithm control parameters
163 setSafety(0.9);
164 setMinReduction(0.2);
165 setMaxGrowth(10.0);
166 }
167
168 /** Get the order of the method.
169 * @return order of the method
170 */
171 public abstract int getOrder();
172
173 /** Get the safety factor for stepsize control.
174 * @return safety factor
175 */
176 public double getSafety() {
177 return safety;
178 }
179
180 /** Set the safety factor for stepsize control.
181 * @param safety safety factor
182 */
183 public void setSafety(final double safety) {
184 this.safety = safety;
185 }
186
187 /** {@inheritDoc} */
188 @Override
189 public void integrate(final ExpandableStatefulODE equations, final double t)
190 throws NumberIsTooSmallException, DimensionMismatchException,
191 MaxCountExceededException, NoBracketingException {
192
193 sanityChecks(equations, t);
194 setEquations(equations);
195 final boolean forward = t > equations.getTime();
196
197 // create some internal working arrays
198 final double[] y0 = equations.getCompleteState();
199 final double[] y = y0.clone();
200 final int stages = c.length + 1;
201 final double[][] yDotK = new double[stages][y.length];
202 final double[] yTmp = y0.clone();
203 final double[] yDotTmp = new double[y.length];
204
205 // set up an interpolator sharing the integrator arrays
206 final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
207 interpolator.reinitialize(this, yTmp, yDotK, forward,
208 equations.getPrimaryMapper(), equations.getSecondaryMappers());
209 interpolator.storeTime(equations.getTime());
210
211 // set up integration control objects
212 stepStart = equations.getTime();
213 double hNew = 0;
214 boolean firstTime = true;
215 initIntegration(equations.getTime(), y0, t);
216
217 // main integration loop
218 isLastStep = false;
219 do {
220
221 interpolator.shift();
222
223 // iterate over step size, ensuring local normalized error is smaller than 1
224 double error = 10;
225 while (error >= 1.0) {
226
227 if (firstTime || !fsal) {
228 // first stage
229 computeDerivatives(stepStart, y, yDotK[0]);
230 }
231
232 if (firstTime) {
233 final double[] scale = new double[mainSetDimension];
234 if (vecAbsoluteTolerance == null) {
235 for (int i = 0; i < scale.length; ++i) {
236 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * JdkMath.abs(y[i]);
237 }
238 } else {
239 for (int i = 0; i < scale.length; ++i) {
240 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * JdkMath.abs(y[i]);
241 }
242 }
243 hNew = initializeStep(forward, getOrder(), scale,
244 stepStart, y, yDotK[0], yTmp, yDotK[1]);
245 firstTime = false;
246 }
247
248 stepSize = hNew;
249 if (forward) {
250 if (stepStart + stepSize >= t) {
251 stepSize = t - stepStart;
252 }
253 } else {
254 if (stepStart + stepSize <= t) {
255 stepSize = t - stepStart;
256 }
257 }
258
259 // next stages
260 for (int k = 1; k < stages; ++k) {
261
262 for (int j = 0; j < y0.length; ++j) {
263 double sum = a[k-1][0] * yDotK[0][j];
264 for (int l = 1; l < k; ++l) {
265 sum += a[k-1][l] * yDotK[l][j];
266 }
267 yTmp[j] = y[j] + stepSize * sum;
268 }
269
270 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
271 }
272
273 // estimate the state at the end of the step
274 for (int j = 0; j < y0.length; ++j) {
275 double sum = b[0] * yDotK[0][j];
276 for (int l = 1; l < stages; ++l) {
277 sum += b[l] * yDotK[l][j];
278 }
279 yTmp[j] = y[j] + stepSize * sum;
280 }
281
282 // estimate the error at the end of the step
283 error = estimateError(yDotK, y, yTmp, stepSize);
284 if (error >= 1.0) {
285 // reject the step and attempt to reduce error by stepsize control
286 final double factor =
287 JdkMath.min(maxGrowth,
288 JdkMath.max(minReduction, safety * JdkMath.pow(error, exp)));
289 hNew = filterStep(stepSize * factor, forward, false);
290 }
291 }
292
293 // local error is small enough: accept the step, trigger events and step handlers
294 interpolator.storeTime(stepStart + stepSize);
295 System.arraycopy(yTmp, 0, y, 0, y0.length);
296 System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
297 stepStart = acceptStep(interpolator, y, yDotTmp, t);
298 System.arraycopy(y, 0, yTmp, 0, y.length);
299
300 if (!isLastStep) {
301
302 // prepare next step
303 interpolator.storeTime(stepStart);
304
305 if (fsal) {
306 // save the last evaluation for the next step
307 System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length);
308 }
309
310 // stepsize control for next step
311 final double factor =
312 JdkMath.min(maxGrowth, JdkMath.max(minReduction, safety * JdkMath.pow(error, exp)));
313 final double scaledH = stepSize * factor;
314 final double nextT = stepStart + scaledH;
315 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
316 hNew = filterStep(scaledH, forward, nextIsLast);
317
318 final double filteredNextT = stepStart + hNew;
319 final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
320 if (filteredNextIsLast) {
321 hNew = t - stepStart;
322 }
323 }
324 } while (!isLastStep);
325
326 // dispatch results
327 equations.setTime(stepStart);
328 equations.setCompleteState(y);
329
330 resetInternalState();
331 }
332
333 /** Get the minimal reduction factor for stepsize control.
334 * @return minimal reduction factor
335 */
336 public double getMinReduction() {
337 return minReduction;
338 }
339
340 /** Set the minimal reduction factor for stepsize control.
341 * @param minReduction minimal reduction factor
342 */
343 public void setMinReduction(final double minReduction) {
344 this.minReduction = minReduction;
345 }
346
347 /** Get the maximal growth factor for stepsize control.
348 * @return maximal growth factor
349 */
350 public double getMaxGrowth() {
351 return maxGrowth;
352 }
353
354 /** Set the maximal growth factor for stepsize control.
355 * @param maxGrowth maximal growth factor
356 */
357 public void setMaxGrowth(final double maxGrowth) {
358 this.maxGrowth = maxGrowth;
359 }
360
361 /** Compute the error ratio.
362 * @param yDotK derivatives computed during the first stages
363 * @param y0 estimate of the step at the start of the step
364 * @param y1 estimate of the step at the end of the step
365 * @param h current step
366 * @return error ratio, greater than 1 if step should be rejected
367 */
368 protected abstract double estimateError(double[][] yDotK,
369 double[] y0, double[] y1,
370 double h);
371 }