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.core.jdkmath.JdkMath;
21
22
23
24
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 public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
57
58
59 private static final String METHOD_NAME = "Dormand-Prince 8 (5, 3)";
60
61
62 private static final double[] STATIC_C = {
63 (12.0 - 2.0 * JdkMath.sqrt(6.0)) / 135.0, (6.0 - JdkMath.sqrt(6.0)) / 45.0, (6.0 - JdkMath.sqrt(6.0)) / 30.0,
64 (6.0 + JdkMath.sqrt(6.0)) / 30.0, 1.0/3.0, 1.0/4.0, 4.0/13.0, 127.0/195.0, 3.0/5.0,
65 6.0/7.0, 1.0, 1.0
66 };
67
68
69 private static final double[][] STATIC_A = {
70
71
72 {(12.0 - 2.0 * JdkMath.sqrt(6.0)) / 135.0},
73
74
75 {(6.0 - JdkMath.sqrt(6.0)) / 180.0, (6.0 - JdkMath.sqrt(6.0)) / 60.0},
76
77
78 {(6.0 - JdkMath.sqrt(6.0)) / 120.0, 0.0, (6.0 - JdkMath.sqrt(6.0)) / 40.0},
79
80
81 {(462.0 + 107.0 * JdkMath.sqrt(6.0)) / 3000.0, 0.0,
82 (-402.0 - 197.0 * JdkMath.sqrt(6.0)) / 1000.0, (168.0 + 73.0 * JdkMath.sqrt(6.0)) / 375.0},
83
84
85 {1.0 / 27.0, 0.0, 0.0, (16.0 + JdkMath.sqrt(6.0)) / 108.0, (16.0 - JdkMath.sqrt(6.0)) / 108.0},
86
87
88 {19.0 / 512.0, 0.0, 0.0, (118.0 + 23.0 * JdkMath.sqrt(6.0)) / 1024.0,
89 (118.0 - 23.0 * JdkMath.sqrt(6.0)) / 1024.0, -9.0 / 512.0},
90
91
92 {13772.0 / 371293.0, 0.0, 0.0, (51544.0 + 4784.0 * JdkMath.sqrt(6.0)) / 371293.0,
93 (51544.0 - 4784.0 * JdkMath.sqrt(6.0)) / 371293.0, -5688.0 / 371293.0, 3072.0 / 371293.0},
94
95
96 {58656157643.0 / 93983540625.0, 0.0, 0.0,
97 (-1324889724104.0 - 318801444819.0 * JdkMath.sqrt(6.0)) / 626556937500.0,
98 (-1324889724104.0 + 318801444819.0 * JdkMath.sqrt(6.0)) / 626556937500.0,
99 96044563816.0 / 3480871875.0, 5682451879168.0 / 281950621875.0,
100 -165125654.0 / 3796875.0},
101
102
103 {8909899.0 / 18653125.0, 0.0, 0.0,
104 (-4521408.0 - 1137963.0 * JdkMath.sqrt(6.0)) / 2937500.0,
105 (-4521408.0 + 1137963.0 * JdkMath.sqrt(6.0)) / 2937500.0,
106 96663078.0 / 4553125.0, 2107245056.0 / 137915625.0,
107 -4913652016.0 / 147609375.0, -78894270.0 / 3880452869.0},
108
109
110 {-20401265806.0 / 21769653311.0, 0.0, 0.0,
111 (354216.0 + 94326.0 * JdkMath.sqrt(6.0)) / 112847.0,
112 (354216.0 - 94326.0 * JdkMath.sqrt(6.0)) / 112847.0,
113 -43306765128.0 / 5313852383.0, -20866708358144.0 / 1126708119789.0,
114 14886003438020.0 / 654632330667.0, 35290686222309375.0 / 14152473387134411.0,
115 -1477884375.0 / 485066827.0},
116
117
118 {39815761.0 / 17514443.0, 0.0, 0.0,
119 (-3457480.0 - 960905.0 * JdkMath.sqrt(6.0)) / 551636.0,
120 (-3457480.0 + 960905.0 * JdkMath.sqrt(6.0)) / 551636.0,
121 -844554132.0 / 47026969.0, 8444996352.0 / 302158619.0,
122 -2509602342.0 / 877790785.0, -28388795297996250.0 / 3199510091356783.0,
123 226716250.0 / 18341897.0, 1371316744.0 / 2131383595.0},
124
125
126
127
128 {104257.0/1920240.0, 0.0, 0.0, 0.0, 0.0, 3399327.0/763840.0,
129 66578432.0/35198415.0, -1674902723.0/288716400.0,
130 54980371265625.0/176692375811392.0, -734375.0/4826304.0,
131 171414593.0/851261400.0, 137909.0/3084480.0}
132 };
133
134
135 private static final double[] STATIC_B = {
136 104257.0/1920240.0,
137 0.0,
138 0.0,
139 0.0,
140 0.0,
141 3399327.0/763840.0,
142 66578432.0/35198415.0,
143 -1674902723.0/288716400.0,
144 54980371265625.0/176692375811392.0,
145 -734375.0/4826304.0,
146 171414593.0/851261400.0,
147 137909.0/3084480.0,
148 0.0
149 };
150
151
152 private static final double E1_01 = 116092271.0 / 8848465920.0;
153
154
155
156
157 private static final double E1_06 = -1871647.0 / 1527680.0;
158
159
160 private static final double E1_07 = -69799717.0 / 140793660.0;
161
162
163 private static final double E1_08 = 1230164450203.0 / 739113984000.0;
164
165
166 private static final double E1_09 = -1980813971228885.0 / 5654156025964544.0;
167
168
169 private static final double E1_10 = 464500805.0 / 1389975552.0;
170
171
172 private static final double E1_11 = 1606764981773.0 / 19613062656000.0;
173
174
175 private static final double E1_12 = -137909.0 / 6168960.0;
176
177
178
179 private static final double E2_01 = -364463.0 / 1920240.0;
180
181
182
183
184 private static final double E2_06 = 3399327.0 / 763840.0;
185
186
187 private static final double E2_07 = 66578432.0 / 35198415.0;
188
189
190 private static final double E2_08 = -1674902723.0 / 288716400.0;
191
192
193 private static final double E2_09 = -74684743568175.0 / 176692375811392.0;
194
195
196 private static final double E2_10 = -734375.0 / 4826304.0;
197
198
199 private static final double E2_11 = 171414593.0 / 851261400.0;
200
201
202 private static final double E2_12 = 69869.0 / 3084480.0;
203
204
205
206
207
208
209
210
211
212
213
214
215 public DormandPrince853Integrator(final double minStep, final double maxStep,
216 final double scalAbsoluteTolerance,
217 final double scalRelativeTolerance) {
218 super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
219 new DormandPrince853StepInterpolator(),
220 minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
221 }
222
223
224
225
226
227
228
229
230
231
232
233
234 public DormandPrince853Integrator(final double minStep, final double maxStep,
235 final double[] vecAbsoluteTolerance,
236 final double[] vecRelativeTolerance) {
237 super(METHOD_NAME, true, STATIC_C, STATIC_A, STATIC_B,
238 new DormandPrince853StepInterpolator(),
239 minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
240 }
241
242
243 @Override
244 public int getOrder() {
245 return 8;
246 }
247
248
249 @Override
250 protected double estimateError(final double[][] yDotK,
251 final double[] y0, final double[] y1,
252 final double h) {
253 double error1 = 0;
254 double error2 = 0;
255
256 for (int j = 0; j < mainSetDimension; ++j) {
257 final double errSum1 = E1_01 * yDotK[0][j] + E1_06 * yDotK[5][j] +
258 E1_07 * yDotK[6][j] + E1_08 * yDotK[7][j] +
259 E1_09 * yDotK[8][j] + E1_10 * yDotK[9][j] +
260 E1_11 * yDotK[10][j] + E1_12 * yDotK[11][j];
261 final double errSum2 = E2_01 * yDotK[0][j] + E2_06 * yDotK[5][j] +
262 E2_07 * yDotK[6][j] + E2_08 * yDotK[7][j] +
263 E2_09 * yDotK[8][j] + E2_10 * yDotK[9][j] +
264 E2_11 * yDotK[10][j] + E2_12 * yDotK[11][j];
265
266 final double yScale = JdkMath.max(JdkMath.abs(y0[j]), JdkMath.abs(y1[j]));
267 final double tol = (vecAbsoluteTolerance == null) ?
268 (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
269 (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
270 final double ratio1 = errSum1 / tol;
271 error1 += ratio1 * ratio1;
272 final double ratio2 = errSum2 / tol;
273 error2 += ratio2 * ratio2;
274 }
275
276 double den = error1 + 0.01 * error2;
277 if (den <= 0.0) {
278 den = 1.0;
279 }
280
281 return JdkMath.abs(h) * error1 / JdkMath.sqrt(mainSetDimension * den);
282 }
283 }