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.exception.MaxCountExceededException;
23 import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
24 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
25 import org.apache.commons.math4.legacy.core.MathArrays;
26
27
28
29
30
31
32
33
34
35
36
37 class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
38 extends RungeKuttaFieldStepInterpolator<T> {
39
40
41
42
43 private final T[][] d;
44
45
46
47
48
49
50
51
52
53
54
55 DormandPrince853FieldStepInterpolator(final Field<T> field, final boolean forward,
56 final T[][] yDotK,
57 final FieldODEStateAndDerivative<T> globalPreviousState,
58 final FieldODEStateAndDerivative<T> globalCurrentState,
59 final FieldODEStateAndDerivative<T> softPreviousState,
60 final FieldODEStateAndDerivative<T> softCurrentState,
61 final FieldEquationsMapper<T> mapper) {
62 super(field, forward, yDotK,
63 globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
64 mapper);
65
66 d = MathArrays.buildArray(field, 7, 16);
67
68
69 d[0][ 0] = fraction(field, 104257, 1920240);
70 d[0][ 1] = field.getZero();
71 d[0][ 2] = field.getZero();
72 d[0][ 3] = field.getZero();
73 d[0][ 4] = field.getZero();
74 d[0][ 5] = fraction(field, 3399327.0, 763840.0);
75 d[0][ 6] = fraction(field, 66578432.0, 35198415.0);
76 d[0][ 7] = fraction(field, -1674902723.0, 288716400.0);
77 d[0][ 8] = fraction(field, 54980371265625.0, 176692375811392.0);
78 d[0][ 9] = fraction(field, -734375.0, 4826304.0);
79 d[0][10] = fraction(field, 171414593.0, 851261400.0);
80 d[0][11] = fraction(field, 137909.0, 3084480.0);
81 d[0][12] = field.getZero();
82 d[0][13] = field.getZero();
83 d[0][14] = field.getZero();
84 d[0][15] = field.getZero();
85
86 d[1][ 0] = d[0][ 0].negate().add(1);
87 d[1][ 1] = d[0][ 1].negate();
88 d[1][ 2] = d[0][ 2].negate();
89 d[1][ 3] = d[0][ 3].negate();
90 d[1][ 4] = d[0][ 4].negate();
91 d[1][ 5] = d[0][ 5].negate();
92 d[1][ 6] = d[0][ 6].negate();
93 d[1][ 7] = d[0][ 7].negate();
94 d[1][ 8] = d[0][ 8].negate();
95 d[1][ 9] = d[0][ 9].negate();
96 d[1][10] = d[0][10].negate();
97 d[1][11] = d[0][11].negate();
98 d[1][12] = d[0][12].negate();
99 d[1][13] = d[0][13].negate();
100 d[1][14] = d[0][14].negate();
101 d[1][15] = d[0][15].negate();
102
103 d[2][ 0] = d[0][ 0].multiply(2).subtract(1);
104 d[2][ 1] = d[0][ 1].multiply(2);
105 d[2][ 2] = d[0][ 2].multiply(2);
106 d[2][ 3] = d[0][ 3].multiply(2);
107 d[2][ 4] = d[0][ 4].multiply(2);
108 d[2][ 5] = d[0][ 5].multiply(2);
109 d[2][ 6] = d[0][ 6].multiply(2);
110 d[2][ 7] = d[0][ 7].multiply(2);
111 d[2][ 8] = d[0][ 8].multiply(2);
112 d[2][ 9] = d[0][ 9].multiply(2);
113 d[2][10] = d[0][10].multiply(2);
114 d[2][11] = d[0][11].multiply(2);
115 d[2][12] = d[0][12].multiply(2).subtract(1);
116 d[2][13] = d[0][13].multiply(2);
117 d[2][14] = d[0][14].multiply(2);
118 d[2][15] = d[0][15].multiply(2);
119
120 d[3][ 0] = fraction(field, -17751989329.0, 2106076560.0);
121 d[3][ 1] = field.getZero();
122 d[3][ 2] = field.getZero();
123 d[3][ 3] = field.getZero();
124 d[3][ 4] = field.getZero();
125 d[3][ 5] = fraction(field, 4272954039.0, 7539864640.0);
126 d[3][ 6] = fraction(field, -118476319744.0, 38604839385.0);
127 d[3][ 7] = fraction(field, 755123450731.0, 316657731600.0);
128 d[3][ 8] = fraction(field, 3692384461234828125.0, 1744130441634250432.0);
129 d[3][ 9] = fraction(field, -4612609375.0, 5293382976.0);
130 d[3][10] = fraction(field, 2091772278379.0, 933644586600.0);
131 d[3][11] = fraction(field, 2136624137.0, 3382989120.0);
132 d[3][12] = fraction(field, -126493.0, 1421424.0);
133 d[3][13] = fraction(field, 98350000.0, 5419179.0);
134 d[3][14] = fraction(field, -18878125.0, 2053168.0);
135 d[3][15] = fraction(field, -1944542619.0, 438351368.0);
136
137 d[4][ 0] = fraction(field, 32941697297.0, 3159114840.0);
138 d[4][ 1] = field.getZero();
139 d[4][ 2] = field.getZero();
140 d[4][ 3] = field.getZero();
141 d[4][ 4] = field.getZero();
142 d[4][ 5] = fraction(field, 456696183123.0, 1884966160.0);
143 d[4][ 6] = fraction(field, 19132610714624.0, 115814518155.0);
144 d[4][ 7] = fraction(field, -177904688592943.0, 474986597400.0);
145 d[4][ 8] = fraction(field, -4821139941836765625.0, 218016305204281304.0);
146 d[4][ 9] = fraction(field, 30702015625.0, 3970037232.0);
147 d[4][10] = fraction(field, -85916079474274.0, 2800933759800.0);
148 d[4][11] = fraction(field, -5919468007.0, 634310460.0);
149 d[4][12] = fraction(field, 2479159.0, 157936.0);
150 d[4][13] = fraction(field, -18750000.0, 602131.0);
151 d[4][14] = fraction(field, -19203125.0, 2053168.0);
152 d[4][15] = fraction(field, 15700361463.0, 438351368.0);
153
154 d[5][ 0] = fraction(field, 12627015655.0, 631822968.0);
155 d[5][ 1] = field.getZero();
156 d[5][ 2] = field.getZero();
157 d[5][ 3] = field.getZero();
158 d[5][ 4] = field.getZero();
159 d[5][ 5] = fraction(field, -72955222965.0, 188496616.0);
160 d[5][ 6] = fraction(field, -13145744952320.0, 69488710893.0);
161 d[5][ 7] = fraction(field, 30084216194513.0, 56998391688.0);
162 d[5][ 8] = fraction(field, -296858761006640625.0, 25648977082856624.0);
163 d[5][ 9] = fraction(field, 569140625.0, 82709109.0);
164 d[5][10] = fraction(field, -18684190637.0, 18672891732.0);
165 d[5][11] = fraction(field, 69644045.0, 89549712.0);
166 d[5][12] = fraction(field, -11847025.0, 4264272.0);
167 d[5][13] = fraction(field, -978650000.0, 16257537.0);
168 d[5][14] = fraction(field, 519371875.0, 6159504.0);
169 d[5][15] = fraction(field, 5256837225.0, 438351368.0);
170
171 d[6][ 0] = fraction(field, -450944925.0, 17550638.0);
172 d[6][ 1] = field.getZero();
173 d[6][ 2] = field.getZero();
174 d[6][ 3] = field.getZero();
175 d[6][ 4] = field.getZero();
176 d[6][ 5] = fraction(field, -14532122925.0, 94248308.0);
177 d[6][ 6] = fraction(field, -595876966400.0, 2573655959.0);
178 d[6][ 7] = fraction(field, 188748653015.0, 527762886.0);
179 d[6][ 8] = fraction(field, 2545485458115234375.0, 27252038150535163.0);
180 d[6][ 9] = fraction(field, -1376953125.0, 36759604.0);
181 d[6][10] = fraction(field, 53995596795.0, 518691437.0);
182 d[6][11] = fraction(field, 210311225.0, 7047894.0);
183 d[6][12] = fraction(field, -1718875.0, 39484.0);
184 d[6][13] = fraction(field, 58000000.0, 602131.0);
185 d[6][14] = fraction(field, -1546875.0, 39484.0);
186 d[6][15] = fraction(field, -1262172375.0, 8429834.0);
187 }
188
189
190 @Override
191 protected DormandPrince853FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
192 final FieldODEStateAndDerivative<T> newGlobalPreviousState,
193 final FieldODEStateAndDerivative<T> newGlobalCurrentState,
194 final FieldODEStateAndDerivative<T> newSoftPreviousState,
195 final FieldODEStateAndDerivative<T> newSoftCurrentState,
196 final FieldEquationsMapper<T> newMapper) {
197 return new DormandPrince853FieldStepInterpolator<>(newField, newForward, newYDotK,
198 newGlobalPreviousState, newGlobalCurrentState,
199 newSoftPreviousState, newSoftCurrentState,
200 newMapper);
201 }
202
203
204
205
206
207
208
209 private T fraction(final Field<T> field, final double p, final double q) {
210 return field.getZero().add(p).divide(q);
211 }
212
213
214 @SuppressWarnings("unchecked")
215 @Override
216 protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
217 final T time, final T theta,
218 final T thetaH, final T oneMinusThetaH)
219 throws MaxCountExceededException {
220
221 final T one = time.getField().getOne();
222 final T eta = one.subtract(theta);
223 final T twoTheta = theta.multiply(2);
224 final T theta2 = theta.multiply(theta);
225 final T dot1 = one.subtract(twoTheta);
226 final T dot2 = theta.multiply(theta.multiply(-3).add(2));
227 final T dot3 = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
228 final T dot4 = theta2.multiply(theta.multiply(theta.multiply(5).subtract(8)).add(3));
229 final T dot5 = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(-6).add(15)).subtract(12)).add(3));
230 final T dot6 = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-7).add(18)).subtract(15)).add(4)));
231 final T[] interpolatedState;
232 final T[] interpolatedDerivatives;
233
234
235 if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
236 final T f0 = thetaH;
237 final T f1 = f0.multiply(eta);
238 final T f2 = f1.multiply(theta);
239 final T f3 = f2.multiply(eta);
240 final T f4 = f3.multiply(theta);
241 final T f5 = f4.multiply(eta);
242 final T f6 = f5.multiply(theta);
243 final T[] p = MathArrays.buildArray(time.getField(), 16);
244 final T[] q = MathArrays.buildArray(time.getField(), 16);
245 for (int i = 0; i < p.length; ++i) {
246 p[i] = f0.multiply(d[0][i]).
247 add(f1.multiply(d[1][i])).
248 add(f2.multiply(d[2][i])).
249 add(f3.multiply(d[3][i])).
250 add(f4.multiply(d[4][i])).
251 add(f5.multiply(d[5][i])).
252 add(f6.multiply(d[6][i]));
253 q[i] = d[0][i].
254 add(dot1.multiply(d[1][i])).
255 add(dot2.multiply(d[2][i])).
256 add(dot3.multiply(d[3][i])).
257 add(dot4.multiply(d[4][i])).
258 add(dot5.multiply(d[5][i])).
259 add(dot6.multiply(d[6][i]));
260 }
261 interpolatedState = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
262 p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
263 interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
264 q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
265 } else {
266 final T f0 = oneMinusThetaH.negate();
267 final T f1 = f0.multiply(theta).negate();
268 final T f2 = f1.multiply(theta);
269 final T f3 = f2.multiply(eta);
270 final T f4 = f3.multiply(theta);
271 final T f5 = f4.multiply(eta);
272 final T f6 = f5.multiply(theta);
273 final T[] p = MathArrays.buildArray(time.getField(), 16);
274 final T[] q = MathArrays.buildArray(time.getField(), 16);
275 for (int i = 0; i < p.length; ++i) {
276 p[i] = f0.multiply(d[0][i]).
277 add(f1.multiply(d[1][i])).
278 add(f2.multiply(d[2][i])).
279 add(f3.multiply(d[3][i])).
280 add(f4.multiply(d[4][i])).
281 add(f5.multiply(d[5][i])).
282 add(f6.multiply(d[6][i]));
283 q[i] = d[0][i].
284 add(dot1.multiply(d[1][i])).
285 add(dot2.multiply(d[2][i])).
286 add(dot3.multiply(d[3][i])).
287 add(dot4.multiply(d[4][i])).
288 add(dot5.multiply(d[5][i])).
289 add(dot6.multiply(d[6][i]));
290 }
291 interpolatedState = currentStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
292 p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
293 interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
294 q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
295 }
296
297 return new FieldODEStateAndDerivative<>(time, interpolatedState, interpolatedDerivatives);
298 }
299 }