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 org.apache.commons.math4.legacy.core.Field;
21 import org.apache.commons.math4.legacy.core.RealFieldElement;
22 import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
23 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
24 import org.apache.commons.math4.legacy.core.MathArrays;
25
26
27
28
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 public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
62 extends EmbeddedRungeKuttaFieldIntegrator<T> {
63
64
65 private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
66
67
68 private final T e1_01;
69
70
71
72
73 private final T e1_06;
74
75
76 private final T e1_07;
77
78
79 private final T e1_08;
80
81
82 private final T e1_09;
83
84
85 private final T e1_10;
86
87
88 private final T e1_11;
89
90
91 private final T e1_12;
92
93
94
95 private final T e2_01;
96
97
98
99
100 private final T e2_06;
101
102
103 private final T e2_07;
104
105
106 private final T e2_08;
107
108
109 private final T e2_09;
110
111
112 private final T e2_10;
113
114
115 private final T e2_11;
116
117
118 private final T e2_12;
119
120
121
122
123
124
125
126
127
128
129
130
131
132 public DormandPrince853FieldIntegrator(final Field<T> field,
133 final double minStep, final double maxStep,
134 final double scalAbsoluteTolerance,
135 final double scalRelativeTolerance) {
136 super(field, METHOD_NAME, 12,
137 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
138 e1_01 = fraction( 116092271.0, 8848465920.0);
139 e1_06 = fraction( -1871647.0, 1527680.0);
140 e1_07 = fraction( -69799717.0, 140793660.0);
141 e1_08 = fraction( 1230164450203.0, 739113984000.0);
142 e1_09 = fraction(-1980813971228885.0, 5654156025964544.0);
143 e1_10 = fraction( 464500805.0, 1389975552.0);
144 e1_11 = fraction( 1606764981773.0, 19613062656000.0);
145 e1_12 = fraction( -137909.0, 6168960.0);
146 e2_01 = fraction( -364463.0, 1920240.0);
147 e2_06 = fraction( 3399327.0, 763840.0);
148 e2_07 = fraction( 66578432.0, 35198415.0);
149 e2_08 = fraction( -1674902723.0, 288716400.0);
150 e2_09 = fraction( -74684743568175.0, 176692375811392.0);
151 e2_10 = fraction( -734375.0, 4826304.0);
152 e2_11 = fraction( 171414593.0, 851261400.0);
153 e2_12 = fraction( 69869.0, 3084480.0);
154 }
155
156
157
158
159
160
161
162
163
164
165
166
167
168 public DormandPrince853FieldIntegrator(final Field<T> field,
169 final double minStep, final double maxStep,
170 final double[] vecAbsoluteTolerance,
171 final double[] vecRelativeTolerance) {
172 super(field, METHOD_NAME, 12,
173 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
174 e1_01 = fraction( 116092271.0, 8848465920.0);
175 e1_06 = fraction( -1871647.0, 1527680.0);
176 e1_07 = fraction( -69799717.0, 140793660.0);
177 e1_08 = fraction( 1230164450203.0, 739113984000.0);
178 e1_09 = fraction(-1980813971228885.0, 5654156025964544.0);
179 e1_10 = fraction( 464500805.0, 1389975552.0);
180 e1_11 = fraction( 1606764981773.0, 19613062656000.0);
181 e1_12 = fraction( -137909.0, 6168960.0);
182 e2_01 = fraction( -364463.0, 1920240.0);
183 e2_06 = fraction( 3399327.0, 763840.0);
184 e2_07 = fraction( 66578432.0, 35198415.0);
185 e2_08 = fraction( -1674902723.0, 288716400.0);
186 e2_09 = fraction( -74684743568175.0, 176692375811392.0);
187 e2_10 = fraction( -734375.0, 4826304.0);
188 e2_11 = fraction( 171414593.0, 851261400.0);
189 e2_12 = fraction( 69869.0, 3084480.0);
190 }
191
192
193 @Override
194 public T[] getC() {
195
196 final T sqrt6 = getField().getOne().multiply(6).sqrt();
197
198 final T[] c = MathArrays.buildArray(getField(), 15);
199 c[ 0] = sqrt6.add(-6).divide(-67.5);
200 c[ 1] = sqrt6.add(-6).divide(-45.0);
201 c[ 2] = sqrt6.add(-6).divide(-30.0);
202 c[ 3] = sqrt6.add( 6).divide( 30.0);
203 c[ 4] = fraction(1, 3);
204 c[ 5] = fraction(1, 4);
205 c[ 6] = fraction(4, 13);
206 c[ 7] = fraction(127, 195);
207 c[ 8] = fraction(3, 5);
208 c[ 9] = fraction(6, 7);
209 c[10] = getField().getOne();
210 c[11] = getField().getOne();
211 c[12] = fraction(1.0, 10.0);
212 c[13] = fraction(1.0, 5.0);
213 c[14] = fraction(7.0, 9.0);
214
215 return c;
216 }
217
218
219 @Override
220 public T[][] getA() {
221
222 final T sqrt6 = getField().getOne().multiply(6).sqrt();
223
224 final T[][] a = MathArrays.buildArray(getField(), 15, -1);
225 for (int i = 0; i < a.length; ++i) {
226 a[i] = MathArrays.buildArray(getField(), i + 1);
227 }
228
229 a[ 0][ 0] = sqrt6.add(-6).divide(-67.5);
230
231 a[ 1][ 0] = sqrt6.add(-6).divide(-180);
232 a[ 1][ 1] = sqrt6.add(-6).divide( -60);
233
234 a[ 2][ 0] = sqrt6.add(-6).divide(-120);
235 a[ 2][ 1] = getField().getZero();
236 a[ 2][ 2] = sqrt6.add(-6).divide( -40);
237
238 a[ 3][ 0] = sqrt6.multiply(107).add(462).divide( 3000);
239 a[ 3][ 1] = getField().getZero();
240 a[ 3][ 2] = sqrt6.multiply(197).add(402).divide(-1000);
241 a[ 3][ 3] = sqrt6.multiply( 73).add(168).divide( 375);
242
243 a[ 4][ 0] = fraction(1, 27);
244 a[ 4][ 1] = getField().getZero();
245 a[ 4][ 2] = getField().getZero();
246 a[ 4][ 3] = sqrt6.add( 16).divide( 108);
247 a[ 4][ 4] = sqrt6.add(-16).divide(-108);
248
249 a[ 5][ 0] = fraction(19, 512);
250 a[ 5][ 1] = getField().getZero();
251 a[ 5][ 2] = getField().getZero();
252 a[ 5][ 3] = sqrt6.multiply( 23).add(118).divide(1024);
253 a[ 5][ 4] = sqrt6.multiply(-23).add(118).divide(1024);
254 a[ 5][ 5] = fraction(-9, 512);
255
256 a[ 6][ 0] = fraction(13772, 371293);
257 a[ 6][ 1] = getField().getZero();
258 a[ 6][ 2] = getField().getZero();
259 a[ 6][ 3] = sqrt6.multiply( 4784).add(51544).divide(371293);
260 a[ 6][ 4] = sqrt6.multiply(-4784).add(51544).divide(371293);
261 a[ 6][ 5] = fraction(-5688, 371293);
262 a[ 6][ 6] = fraction( 3072, 371293);
263
264 a[ 7][ 0] = fraction(58656157643.0, 93983540625.0);
265 a[ 7][ 1] = getField().getZero();
266 a[ 7][ 2] = getField().getZero();
267 a[ 7][ 3] = sqrt6.multiply(-318801444819.0).add(-1324889724104.0).divide(626556937500.0);
268 a[ 7][ 4] = sqrt6.multiply( 318801444819.0).add(-1324889724104.0).divide(626556937500.0);
269 a[ 7][ 5] = fraction(96044563816.0, 3480871875.0);
270 a[ 7][ 6] = fraction(5682451879168.0, 281950621875.0);
271 a[ 7][ 7] = fraction(-165125654.0, 3796875.0);
272
273 a[ 8][ 0] = fraction(8909899.0, 18653125.0);
274 a[ 8][ 1] = getField().getZero();
275 a[ 8][ 2] = getField().getZero();
276 a[ 8][ 3] = sqrt6.multiply(-1137963.0).add(-4521408.0).divide(2937500.0);
277 a[ 8][ 4] = sqrt6.multiply( 1137963.0).add(-4521408.0).divide(2937500.0);
278 a[ 8][ 5] = fraction(96663078.0, 4553125.0);
279 a[ 8][ 6] = fraction(2107245056.0, 137915625.0);
280 a[ 8][ 7] = fraction(-4913652016.0, 147609375.0);
281 a[ 8][ 8] = fraction(-78894270.0, 3880452869.0);
282
283 a[ 9][ 0] = fraction(-20401265806.0, 21769653311.0);
284 a[ 9][ 1] = getField().getZero();
285 a[ 9][ 2] = getField().getZero();
286 a[ 9][ 3] = sqrt6.multiply( 94326.0).add(354216.0).divide(112847.0);
287 a[ 9][ 4] = sqrt6.multiply(-94326.0).add(354216.0).divide(112847.0);
288 a[ 9][ 5] = fraction(-43306765128.0, 5313852383.0);
289 a[ 9][ 6] = fraction(-20866708358144.0, 1126708119789.0);
290 a[ 9][ 7] = fraction(14886003438020.0, 654632330667.0);
291 a[ 9][ 8] = fraction(35290686222309375.0, 14152473387134411.0);
292 a[ 9][ 9] = fraction(-1477884375.0, 485066827.0);
293
294 a[10][ 0] = fraction(39815761.0, 17514443.0);
295 a[10][ 1] = getField().getZero();
296 a[10][ 2] = getField().getZero();
297 a[10][ 3] = sqrt6.multiply(-960905.0).add(-3457480.0).divide(551636.0);
298 a[10][ 4] = sqrt6.multiply( 960905.0).add(-3457480.0).divide(551636.0);
299 a[10][ 5] = fraction(-844554132.0, 47026969.0);
300 a[10][ 6] = fraction(8444996352.0, 302158619.0);
301 a[10][ 7] = fraction(-2509602342.0, 877790785.0);
302 a[10][ 8] = fraction(-28388795297996250.0, 3199510091356783.0);
303 a[10][ 9] = fraction(226716250.0, 18341897.0);
304 a[10][10] = fraction(1371316744.0, 2131383595.0);
305
306
307
308 a[11][ 0] = fraction(104257.0, 1920240.0);
309 a[11][ 1] = getField().getZero();
310 a[11][ 2] = getField().getZero();
311 a[11][ 3] = getField().getZero();
312 a[11][ 4] = getField().getZero();
313 a[11][ 5] = fraction(3399327.0, 763840.0);
314 a[11][ 6] = fraction(66578432.0, 35198415.0);
315 a[11][ 7] = fraction(-1674902723.0, 288716400.0);
316 a[11][ 8] = fraction(54980371265625.0, 176692375811392.0);
317 a[11][ 9] = fraction(-734375.0, 4826304.0);
318 a[11][10] = fraction(171414593.0, 851261400.0);
319 a[11][11] = fraction(137909.0, 3084480.0);
320
321
322 a[12][ 0] = fraction( 13481885573.0, 240030000000.0);
323 a[12][ 1] = getField().getZero();
324 a[12][ 2] = getField().getZero();
325 a[12][ 3] = getField().getZero();
326 a[12][ 4] = getField().getZero();
327 a[12][ 5] = getField().getZero();
328 a[12][ 6] = fraction( 139418837528.0, 549975234375.0);
329 a[12][ 7] = fraction( -11108320068443.0, 45111937500000.0);
330 a[12][ 8] = fraction(-1769651421925959.0, 14249385146080000.0);
331 a[12][ 9] = fraction( 57799439.0, 377055000.0);
332 a[12][10] = fraction( 793322643029.0, 96734250000000.0);
333 a[12][11] = fraction( 1458939311.0, 192780000000.0);
334 a[12][12] = fraction( -4149.0, 500000.0);
335
336 a[13][ 0] = fraction( 1595561272731.0, 50120273500000.0);
337 a[13][ 1] = getField().getZero();
338 a[13][ 2] = getField().getZero();
339 a[13][ 3] = getField().getZero();
340 a[13][ 4] = getField().getZero();
341 a[13][ 5] = fraction( 975183916491.0, 34457688031250.0);
342 a[13][ 6] = fraction( 38492013932672.0, 718912673015625.0);
343 a[13][ 7] = fraction(-1114881286517557.0, 20298710767500000.0);
344 a[13][ 8] = getField().getZero();
345 a[13][ 9] = getField().getZero();
346 a[13][10] = fraction( -2538710946863.0, 23431227861250000.0);
347 a[13][11] = fraction( 8824659001.0, 23066716781250.0);
348 a[13][12] = fraction( -11518334563.0, 33831184612500.0);
349 a[13][13] = fraction( 1912306948.0, 13532473845.0);
350
351 a[14][ 0] = fraction( -13613986967.0, 31741908048.0);
352 a[14][ 1] = getField().getZero();
353 a[14][ 2] = getField().getZero();
354 a[14][ 3] = getField().getZero();
355 a[14][ 4] = getField().getZero();
356 a[14][ 5] = fraction( -4755612631.0, 1012344804.0);
357 a[14][ 6] = fraction( 42939257944576.0, 5588559685701.0);
358 a[14][ 7] = fraction( 77881972900277.0, 19140370552944.0);
359 a[14][ 8] = fraction( 22719829234375.0, 63689648654052.0);
360 a[14][ 9] = getField().getZero();
361 a[14][10] = getField().getZero();
362 a[14][11] = getField().getZero();
363 a[14][12] = fraction( -1199007803.0, 857031517296.0);
364 a[14][13] = fraction( 157882067000.0, 53564469831.0);
365 a[14][14] = fraction( -290468882375.0, 31741908048.0);
366
367 return a;
368 }
369
370
371 @Override
372 public T[] getB() {
373 final T[] b = MathArrays.buildArray(getField(), 16);
374 b[ 0] = fraction(104257, 1920240);
375 b[ 1] = getField().getZero();
376 b[ 2] = getField().getZero();
377 b[ 3] = getField().getZero();
378 b[ 4] = getField().getZero();
379 b[ 5] = fraction( 3399327.0, 763840.0);
380 b[ 6] = fraction( 66578432.0, 35198415.0);
381 b[ 7] = fraction( -1674902723.0, 288716400.0);
382 b[ 8] = fraction( 54980371265625.0, 176692375811392.0);
383 b[ 9] = fraction( -734375.0, 4826304.0);
384 b[10] = fraction( 171414593.0, 851261400.0);
385 b[11] = fraction( 137909.0, 3084480.0);
386 b[12] = getField().getZero();
387 b[13] = getField().getZero();
388 b[14] = getField().getZero();
389 b[15] = getField().getZero();
390 return b;
391 }
392
393
394 @Override
395 protected DormandPrince853FieldStepInterpolator<T>
396 createInterpolator(final boolean forward, T[][] yDotK,
397 final FieldODEStateAndDerivative<T> globalPreviousState,
398 final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
399 return new DormandPrince853FieldStepInterpolator<>(getField(), forward, yDotK,
400 globalPreviousState, globalCurrentState,
401 globalPreviousState, globalCurrentState,
402 mapper);
403 }
404
405
406 @Override
407 public int getOrder() {
408 return 8;
409 }
410
411
412 @Override
413 protected T estimateError(final T[][] yDotK, final T[] y0, final T[] y1, final T h) {
414 T error1 = h.getField().getZero();
415 T error2 = h.getField().getZero();
416
417 for (int j = 0; j < mainSetDimension; ++j) {
418 final T errSum1 = yDotK[ 0][j].multiply(e1_01).
419 add(yDotK[ 5][j].multiply(e1_06)).
420 add(yDotK[ 6][j].multiply(e1_07)).
421 add(yDotK[ 7][j].multiply(e1_08)).
422 add(yDotK[ 8][j].multiply(e1_09)).
423 add(yDotK[ 9][j].multiply(e1_10)).
424 add(yDotK[10][j].multiply(e1_11)).
425 add(yDotK[11][j].multiply(e1_12));
426 final T errSum2 = yDotK[ 0][j].multiply(e2_01).
427 add(yDotK[ 5][j].multiply(e2_06)).
428 add(yDotK[ 6][j].multiply(e2_07)).
429 add(yDotK[ 7][j].multiply(e2_08)).
430 add(yDotK[ 8][j].multiply(e2_09)).
431 add(yDotK[ 9][j].multiply(e2_10)).
432 add(yDotK[10][j].multiply(e2_11)).
433 add(yDotK[11][j].multiply(e2_12));
434
435 final T yScale = RealFieldElement.max(y0[j].abs(), y1[j].abs());
436 final T tol = vecAbsoluteTolerance == null ?
437 yScale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) :
438 yScale.multiply(vecRelativeTolerance[j]).add(vecAbsoluteTolerance[j]);
439 final T ratio1 = errSum1.divide(tol);
440 error1 = error1.add(ratio1.multiply(ratio1));
441 final T ratio2 = errSum2.divide(tol);
442 error2 = error2.add(ratio2.multiply(ratio2));
443 }
444
445 T den = error1.add(error2.multiply(0.01));
446 if (den.getReal() <= 0.0) {
447 den = h.getField().getOne();
448 }
449
450 return h.abs().multiply(error1).divide(den.multiply(mainSetDimension).sqrt());
451 }
452 }