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