1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.ode.nonstiff;
19
20 import java.io.IOException;
21 import java.io.ObjectInput;
22 import java.io.ObjectOutput;
23
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
26 import org.apache.commons.math4.legacy.ode.EquationsMapper;
27 import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
28
29
30
31
32
33
34
35
36
37
38 class DormandPrince853StepInterpolator
39 extends RungeKuttaStepInterpolator {
40
41
42 private static final long serialVersionUID = 20111120L;
43
44
45 private static final double B_01 = 104257.0 / 1920240.0;
46
47
48
49
50 private static final double B_06 = 3399327.0 / 763840.0;
51
52
53 private static final double B_07 = 66578432.0 / 35198415.0;
54
55
56 private static final double B_08 = -1674902723.0 / 288716400.0;
57
58
59 private static final double B_09 = 54980371265625.0 / 176692375811392.0;
60
61
62 private static final double B_10 = -734375.0 / 4826304.0;
63
64
65 private static final double B_11 = 171414593.0 / 851261400.0;
66
67
68 private static final double B_12 = 137909.0 / 3084480.0;
69
70
71 private static final double C14 = 1.0 / 10.0;
72
73
74 private static final double K14_01 = 13481885573.0 / 240030000000.0 - B_01;
75
76
77
78
79 private static final double K14_06 = 0.0 - B_06;
80
81
82 private static final double K14_07 = 139418837528.0 / 549975234375.0 - B_07;
83
84
85 private static final double K14_08 = -11108320068443.0 / 45111937500000.0 - B_08;
86
87
88 private static final double K14_09 = -1769651421925959.0 / 14249385146080000.0 - B_09;
89
90
91 private static final double K14_10 = 57799439.0 / 377055000.0 - B_10;
92
93
94 private static final double K14_11 = 793322643029.0 / 96734250000000.0 - B_11;
95
96
97 private static final double K14_12 = 1458939311.0 / 192780000000.0 - B_12;
98
99
100 private static final double K14_13 = -4149.0 / 500000.0;
101
102
103 private static final double C15 = 1.0 / 5.0;
104
105
106
107 private static final double K15_01 = 1595561272731.0 / 50120273500000.0 - B_01;
108
109
110
111
112 private static final double K15_06 = 975183916491.0 / 34457688031250.0 - B_06;
113
114
115 private static final double K15_07 = 38492013932672.0 / 718912673015625.0 - B_07;
116
117
118 private static final double K15_08 = -1114881286517557.0 / 20298710767500000.0 - B_08;
119
120
121 private static final double K15_09 = 0.0 - B_09;
122
123
124 private static final double K15_10 = 0.0 - B_10;
125
126
127 private static final double K15_11 = -2538710946863.0 / 23431227861250000.0 - B_11;
128
129
130 private static final double K15_12 = 8824659001.0 / 23066716781250.0 - B_12;
131
132
133 private static final double K15_13 = -11518334563.0 / 33831184612500.0;
134
135
136 private static final double K15_14 = 1912306948.0 / 13532473845.0;
137
138
139 private static final double C16 = 7.0 / 9.0;
140
141
142
143 private static final double K16_01 = -13613986967.0 / 31741908048.0 - B_01;
144
145
146
147
148 private static final double K16_06 = -4755612631.0 / 1012344804.0 - B_06;
149
150
151 private static final double K16_07 = 42939257944576.0 / 5588559685701.0 - B_07;
152
153
154 private static final double K16_08 = 77881972900277.0 / 19140370552944.0 - B_08;
155
156
157 private static final double K16_09 = 22719829234375.0 / 63689648654052.0 - B_09;
158
159
160 private static final double K16_10 = 0.0 - B_10;
161
162
163 private static final double K16_11 = 0.0 - B_11;
164
165
166 private static final double K16_12 = 0.0 - B_12;
167
168
169 private static final double K16_13 = -1199007803.0 / 857031517296.0;
170
171
172 private static final double K16_14 = 157882067000.0 / 53564469831.0;
173
174
175 private static final double K16_15 = -290468882375.0 / 31741908048.0;
176
177
178
179
180 private static final double[][] D = {
181
182 { -17751989329.0 / 2106076560.0, 4272954039.0 / 7539864640.0,
183 -118476319744.0 / 38604839385.0, 755123450731.0 / 316657731600.0,
184 3692384461234828125.0 / 1744130441634250432.0, -4612609375.0 / 5293382976.0,
185 2091772278379.0 / 933644586600.0, 2136624137.0 / 3382989120.0,
186 -126493.0 / 1421424.0, 98350000.0 / 5419179.0,
187 -18878125.0 / 2053168.0, -1944542619.0 / 438351368.0},
188
189 { 32941697297.0 / 3159114840.0, 456696183123.0 / 1884966160.0,
190 19132610714624.0 / 115814518155.0, -177904688592943.0 / 474986597400.0,
191 -4821139941836765625.0 / 218016305204281304.0, 30702015625.0 / 3970037232.0,
192 -85916079474274.0 / 2800933759800.0, -5919468007.0 / 634310460.0,
193 2479159.0 / 157936.0, -18750000.0 / 602131.0,
194 -19203125.0 / 2053168.0, 15700361463.0 / 438351368.0},
195
196 { 12627015655.0 / 631822968.0, -72955222965.0 / 188496616.0,
197 -13145744952320.0 / 69488710893.0, 30084216194513.0 / 56998391688.0,
198 -296858761006640625.0 / 25648977082856624.0, 569140625.0 / 82709109.0,
199 -18684190637.0 / 18672891732.0, 69644045.0 / 89549712.0,
200 -11847025.0 / 4264272.0, -978650000.0 / 16257537.0,
201 519371875.0 / 6159504.0, 5256837225.0 / 438351368.0},
202
203 { -450944925.0 / 17550638.0, -14532122925.0 / 94248308.0,
204 -595876966400.0 / 2573655959.0, 188748653015.0 / 527762886.0,
205 2545485458115234375.0 / 27252038150535163.0, -1376953125.0 / 36759604.0,
206 53995596795.0 / 518691437.0, 210311225.0 / 7047894.0,
207 -1718875.0 / 39484.0, 58000000.0 / 602131.0,
208 -1546875.0 / 39484.0, -1262172375.0 / 8429834.0}
209 };
210
211
212 private double[][] yDotKLast;
213
214
215 private double[][] v;
216
217
218 private boolean vectorsInitialized;
219
220
221
222
223
224
225
226
227
228
229
230
231 public DormandPrince853StepInterpolator() {
232 super();
233 yDotKLast = null;
234 v = null;
235 vectorsInitialized = false;
236 }
237
238
239
240
241
242
243
244 DormandPrince853StepInterpolator(final DormandPrince853StepInterpolator interpolator) {
245
246 super(interpolator);
247
248 if (interpolator.currentState == null) {
249
250 yDotKLast = null;
251 v = null;
252 vectorsInitialized = false;
253 } else {
254
255 final int dimension = interpolator.currentState.length;
256
257 yDotKLast = new double[3][];
258 for (int k = 0; k < yDotKLast.length; ++k) {
259 yDotKLast[k] = new double[dimension];
260 System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
261 dimension);
262 }
263
264 v = new double[7][];
265 for (int k = 0; k < v.length; ++k) {
266 v[k] = new double[dimension];
267 System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
268 }
269
270 vectorsInitialized = interpolator.vectorsInitialized;
271 }
272 }
273
274
275 @Override
276 protected StepInterpolator doCopy() {
277 return new DormandPrince853StepInterpolator(this);
278 }
279
280
281 @Override
282 public void reinitialize(final AbstractIntegrator integrator,
283 final double[] y, final double[][] yDotK, final boolean forward,
284 final EquationsMapper primaryMapper,
285 final EquationsMapper[] secondaryMappers) {
286
287 super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers);
288
289 final int dimension = currentState.length;
290
291 yDotKLast = new double[3][];
292 for (int k = 0; k < yDotKLast.length; ++k) {
293 yDotKLast[k] = new double[dimension];
294 }
295
296 v = new double[7][];
297 for (int k = 0; k < v.length; ++k) {
298 v[k] = new double[dimension];
299 }
300
301 vectorsInitialized = false;
302 }
303
304
305 @Override
306 public void storeTime(final double t) {
307 super.storeTime(t);
308 vectorsInitialized = false;
309 }
310
311
312 @Override
313 protected void computeInterpolatedStateAndDerivatives(final double theta,
314 final double oneMinusThetaH)
315 throws MaxCountExceededException {
316
317 if (! vectorsInitialized) {
318
319 if (v == null) {
320 v = new double[7][];
321 for (int k = 0; k < 7; ++k) {
322 v[k] = new double[interpolatedState.length];
323 }
324 }
325
326
327 finalizeStep();
328
329
330 for (int i = 0; i < interpolatedState.length; ++i) {
331 final double yDot1 = yDotK[0][i];
332 final double yDot6 = yDotK[5][i];
333 final double yDot7 = yDotK[6][i];
334 final double yDot8 = yDotK[7][i];
335 final double yDot9 = yDotK[8][i];
336 final double yDot10 = yDotK[9][i];
337 final double yDot11 = yDotK[10][i];
338 final double yDot12 = yDotK[11][i];
339 final double yDot13 = yDotK[12][i];
340 final double yDot14 = yDotKLast[0][i];
341 final double yDot15 = yDotKLast[1][i];
342 final double yDot16 = yDotKLast[2][i];
343 v[0][i] = B_01 * yDot1 + B_06 * yDot6 + B_07 * yDot7 +
344 B_08 * yDot8 + B_09 * yDot9 + B_10 * yDot10 +
345 B_11 * yDot11 + B_12 * yDot12;
346 v[1][i] = yDot1 - v[0][i];
347 v[2][i] = v[0][i] - v[1][i] - yDotK[12][i];
348 for (int k = 0; k < D.length; ++k) {
349 v[k+3][i] = D[k][0] * yDot1 + D[k][1] * yDot6 + D[k][2] * yDot7 +
350 D[k][3] * yDot8 + D[k][4] * yDot9 + D[k][5] * yDot10 +
351 D[k][6] * yDot11 + D[k][7] * yDot12 + D[k][8] * yDot13 +
352 D[k][9] * yDot14 + D[k][10] * yDot15 + D[k][11] * yDot16;
353 }
354 }
355
356 vectorsInitialized = true;
357 }
358
359 final double eta = 1 - theta;
360 final double twoTheta = 2 * theta;
361 final double theta2 = theta * theta;
362 final double dot1 = 1 - twoTheta;
363 final double dot2 = theta * (2 - 3 * theta);
364 final double dot3 = twoTheta * (1 + theta * (twoTheta -3));
365 final double dot4 = theta2 * (3 + theta * (5 * theta - 8));
366 final double dot5 = theta2 * (3 + theta * (-12 + theta * (15 - 6 * theta)));
367 final double dot6 = theta2 * theta * (4 + theta * (-15 + theta * (18 - 7 * theta)));
368
369 if (previousState != null && theta <= 0.5) {
370 for (int i = 0; i < interpolatedState.length; ++i) {
371 interpolatedState[i] = previousState[i] +
372 theta * h * (v[0][i] +
373 eta * (v[1][i] +
374 theta * (v[2][i] +
375 eta * (v[3][i] +
376 theta * (v[4][i] +
377 eta * (v[5][i] +
378 theta * (v[6][i])))))));
379 interpolatedDerivatives[i] = v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] +
380 dot3 * v[3][i] + dot4 * v[4][i] +
381 dot5 * v[5][i] + dot6 * v[6][i];
382 }
383 } else {
384 for (int i = 0; i < interpolatedState.length; ++i) {
385 interpolatedState[i] = currentState[i] -
386 oneMinusThetaH * (v[0][i] -
387 theta * (v[1][i] +
388 theta * (v[2][i] +
389 eta * (v[3][i] +
390 theta * (v[4][i] +
391 eta * (v[5][i] +
392 theta * (v[6][i])))))));
393 interpolatedDerivatives[i] = v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] +
394 dot3 * v[3][i] + dot4 * v[4][i] +
395 dot5 * v[5][i] + dot6 * v[6][i];
396 }
397 }
398 }
399
400
401 @Override
402 protected void doFinalize() throws MaxCountExceededException {
403
404 if (currentState == null) {
405
406 return;
407 }
408
409 double s;
410 final double[] yTmp = new double[currentState.length];
411 final double pT = getGlobalPreviousTime();
412
413
414 for (int j = 0; j < currentState.length; ++j) {
415 s = K14_01 * yDotK[0][j] + K14_06 * yDotK[5][j] + K14_07 * yDotK[6][j] +
416 K14_08 * yDotK[7][j] + K14_09 * yDotK[8][j] + K14_10 * yDotK[9][j] +
417 K14_11 * yDotK[10][j] + K14_12 * yDotK[11][j] + K14_13 * yDotK[12][j];
418 yTmp[j] = currentState[j] + h * s;
419 }
420 integrator.computeDerivatives(pT + C14 * h, yTmp, yDotKLast[0]);
421
422
423 for (int j = 0; j < currentState.length; ++j) {
424 s = K15_01 * yDotK[0][j] + K15_06 * yDotK[5][j] + K15_07 * yDotK[6][j] +
425 K15_08 * yDotK[7][j] + K15_09 * yDotK[8][j] + K15_10 * yDotK[9][j] +
426 K15_11 * yDotK[10][j] + K15_12 * yDotK[11][j] + K15_13 * yDotK[12][j] +
427 K15_14 * yDotKLast[0][j];
428 yTmp[j] = currentState[j] + h * s;
429 }
430 integrator.computeDerivatives(pT + C15 * h, yTmp, yDotKLast[1]);
431
432
433 for (int j = 0; j < currentState.length; ++j) {
434 s = K16_01 * yDotK[0][j] + K16_06 * yDotK[5][j] + K16_07 * yDotK[6][j] +
435 K16_08 * yDotK[7][j] + K16_09 * yDotK[8][j] + K16_10 * yDotK[9][j] +
436 K16_11 * yDotK[10][j] + K16_12 * yDotK[11][j] + K16_13 * yDotK[12][j] +
437 K16_14 * yDotKLast[0][j] + K16_15 * yDotKLast[1][j];
438 yTmp[j] = currentState[j] + h * s;
439 }
440 integrator.computeDerivatives(pT + C16 * h, yTmp, yDotKLast[2]);
441 }
442
443
444 @Override
445 public void writeExternal(final ObjectOutput out)
446 throws IOException {
447
448 try {
449
450 finalizeStep();
451 } catch (MaxCountExceededException mcee) {
452 final IOException ioe = new IOException(mcee.getLocalizedMessage());
453 ioe.initCause(mcee);
454 throw ioe;
455 }
456
457 final int dimension = (currentState == null) ? -1 : currentState.length;
458 out.writeInt(dimension);
459 for (int i = 0; i < dimension; ++i) {
460 out.writeDouble(yDotKLast[0][i]);
461 out.writeDouble(yDotKLast[1][i]);
462 out.writeDouble(yDotKLast[2][i]);
463 }
464
465
466 super.writeExternal(out);
467 }
468
469
470 @Override
471 public void readExternal(final ObjectInput in)
472 throws IOException, ClassNotFoundException {
473
474
475 yDotKLast = new double[3][];
476 final int dimension = in.readInt();
477 yDotKLast[0] = (dimension < 0) ? null : new double[dimension];
478 yDotKLast[1] = (dimension < 0) ? null : new double[dimension];
479 yDotKLast[2] = (dimension < 0) ? null : new double[dimension];
480
481 for (int i = 0; i < dimension; ++i) {
482 yDotKLast[0][i] = in.readDouble();
483 yDotKLast[1][i] = in.readDouble();
484 yDotKLast[2][i] = in.readDouble();
485 }
486
487
488 super.readExternal(in);
489 }
490 }