1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67 public abstract class EmbeddedRungeKuttaIntegrator
68 extends AdaptiveStepsizeIntegrator {
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
85 final double[] c, final double[][] a, final double[] b,
86 final RungeKuttaStepInterpolator prototype,
87 final double minStep, final double maxStep,
88 final double scalAbsoluteTolerance,
89 final double scalRelativeTolerance) {
90
91 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
92
93 this.fsal = fsal;
94 this.c = c;
95 this.a = a;
96 this.b = b;
97 this.prototype = prototype;
98
99 exp = -1.0 / getOrder();
100
101
102 setSafety(0.9);
103 setMinReduction(0.2);
104 setMaxGrowth(10.0);
105
106 }
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122 protected EmbeddedRungeKuttaIntegrator(final String name, final boolean fsal,
123 final double[] c, final double[][] a, final double[] b,
124 final RungeKuttaStepInterpolator prototype,
125 final double minStep, final double maxStep,
126 final double[] vecAbsoluteTolerance,
127 final double[] vecRelativeTolerance) {
128
129 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
130
131 this.fsal = fsal;
132 this.c = c;
133 this.a = a;
134 this.b = b;
135 this.prototype = prototype;
136
137 exp = -1.0 / getOrder();
138
139
140 setSafety(0.9);
141 setMinReduction(0.2);
142 setMaxGrowth(10.0);
143
144 }
145
146
147
148
149 public abstract int getOrder();
150
151
152
153
154 public double getSafety() {
155 return safety;
156 }
157
158
159
160
161 public void setSafety(final double safety) {
162 this.safety = safety;
163 }
164
165
166 @Override
167 public double integrate(final FirstOrderDifferentialEquations equations,
168 final double t0, final double[] y0,
169 final double t, final double[] y)
170 throws DerivativeException, IntegratorException {
171
172 sanityChecks(equations, t0, y0, t, y);
173 setEquations(equations);
174 resetEvaluations();
175 final boolean forward = (t > t0);
176
177
178 final int stages = c.length + 1;
179 if (y != y0) {
180 System.arraycopy(y0, 0, y, 0, y0.length);
181 }
182 final double[][] yDotK = new double[stages][y0.length];
183 final double[] yTmp = new double[y0.length];
184
185
186 AbstractStepInterpolator interpolator;
187 if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
188 final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
189 rki.reinitialize(this, yTmp, yDotK, forward);
190 interpolator = rki;
191 } else {
192 interpolator = new DummyStepInterpolator(yTmp, forward);
193 }
194 interpolator.storeTime(t0);
195
196
197 stepStart = t0;
198 double hNew = 0;
199 boolean firstTime = true;
200 for (StepHandler handler : stepHandlers) {
201 handler.reset();
202 }
203 CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
204 boolean lastStep = false;
205
206
207 while (!lastStep) {
208
209 interpolator.shift();
210
211 double error = 0;
212 for (boolean loop = true; loop;) {
213
214 if (firstTime || !fsal) {
215
216 computeDerivatives(stepStart, y, yDotK[0]);
217 }
218
219 if (firstTime) {
220 final double[] scale;
221 if (vecAbsoluteTolerance != null) {
222 scale = vecAbsoluteTolerance;
223 } else {
224 scale = new double[y0.length];
225 Arrays.fill(scale, scalAbsoluteTolerance);
226 }
227 hNew = initializeStep(equations, forward, getOrder(), scale,
228 stepStart, y, yDotK[0], yTmp, yDotK[1]);
229 firstTime = false;
230 }
231
232 stepSize = hNew;
233
234
235 for (int k = 1; k < stages; ++k) {
236
237 for (int j = 0; j < y0.length; ++j) {
238 double sum = a[k-1][0] * yDotK[0][j];
239 for (int l = 1; l < k; ++l) {
240 sum += a[k-1][l] * yDotK[l][j];
241 }
242 yTmp[j] = y[j] + stepSize * sum;
243 }
244
245 computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
246
247 }
248
249
250 for (int j = 0; j < y0.length; ++j) {
251 double sum = b[0] * yDotK[0][j];
252 for (int l = 1; l < stages; ++l) {
253 sum += b[l] * yDotK[l][j];
254 }
255 yTmp[j] = y[j] + stepSize * sum;
256 }
257
258
259 error = estimateError(yDotK, y, yTmp, stepSize);
260 if (error <= 1.0) {
261
262
263 interpolator.storeTime(stepStart + stepSize);
264 if (manager.evaluateStep(interpolator)) {
265 final double dt = manager.getEventTime() - stepStart;
266 if (Math.abs(dt) <= Math.ulp(stepStart)) {
267
268 loop = false;
269 } else {
270
271 hNew = dt;
272 }
273 } else {
274
275 loop = false;
276 }
277
278 } else {
279
280 final double factor =
281 Math.min(maxGrowth,
282 Math.max(minReduction, safety * Math.pow(error, exp)));
283 hNew = filterStep(stepSize * factor, forward, false);
284 }
285
286 }
287
288
289 final double nextStep = stepStart + stepSize;
290 System.arraycopy(yTmp, 0, y, 0, y0.length);
291 manager.stepAccepted(nextStep, y);
292 lastStep = manager.stop();
293
294
295 interpolator.storeTime(nextStep);
296 for (StepHandler handler : stepHandlers) {
297 handler.handleStep(interpolator, lastStep);
298 }
299 stepStart = nextStep;
300
301 if (fsal) {
302
303 System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
304 }
305
306 if (manager.reset(stepStart, y) && ! lastStep) {
307
308
309 computeDerivatives(stepStart, y, yDotK[0]);
310 }
311
312 if (! lastStep) {
313
314
315
316 stepSize = filterStep(stepSize, forward, true);
317
318
319 final double factor = Math.min(maxGrowth,
320 Math.max(minReduction,
321 safety * Math.pow(error, exp)));
322 final double scaledH = stepSize * factor;
323 final double nextT = stepStart + scaledH;
324 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
325 hNew = filterStep(scaledH, forward, nextIsLast);
326 }
327
328 }
329
330 final double stopTime = stepStart;
331 resetInternalState();
332 return stopTime;
333
334 }
335
336
337
338
339 public double getMinReduction() {
340 return minReduction;
341 }
342
343
344
345
346 public void setMinReduction(final double minReduction) {
347 this.minReduction = minReduction;
348 }
349
350
351
352
353 public double getMaxGrowth() {
354 return maxGrowth;
355 }
356
357
358
359
360 public void setMaxGrowth(final double maxGrowth) {
361 this.maxGrowth = maxGrowth;
362 }
363
364
365
366
367
368
369
370
371 protected abstract double estimateError(double[][] yDotK,
372 double[] y0, double[] y1,
373 double h);
374
375
376 private boolean fsal;
377
378
379 private double[] c;
380
381
382 private double[][] a;
383
384
385 private double[] b;
386
387
388 private RungeKuttaStepInterpolator prototype;
389
390
391 private double exp;
392
393
394 private double safety;
395
396
397 private double minReduction;
398
399
400 private double maxGrowth;
401
402 }