1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.hipparchus.ode.nonstiff;
19
20 import org.hipparchus.ode.EquationsMapper;
21 import org.hipparchus.ode.ODEStateAndDerivative;
22
23
24
25
26
27
28
29
30
31 class DormandPrince853StateInterpolator
32 extends RungeKuttaStateInterpolator {
33
34
35 private static final long serialVersionUID = 20160328L;
36
37
38 private static final double[][] D = {
39 {
40
41 104257.0 / 1920240.0,
42 0.0,
43 0.0,
44 0.0,
45 0.0,
46 3399327.0 / 763840.0,
47 66578432.0 / 35198415.0,
48 -1674902723.0 / 288716400.0,
49 54980371265625.0 / 176692375811392.0,
50 -734375.0 / 4826304.0,
51 171414593.0 / 851261400.0,
52 137909.0 / 3084480.0,
53 0.0,
54 0.0,
55 0.0,
56 0.0,
57 }, {
58 1815983.0 / 1920240.0,
59 0.0,
60 0.0,
61 0.0,
62 0.0,
63 -3399327.0 / 763840.0,
64 -66578432.0 / 35198415.0,
65 1674902723.0 / 288716400.0,
66 -54980371265625.0 / 176692375811392.0,
67 734375.0 / 4826304.0,
68 -171414593.0 / 851261400.0,
69 -137909.0 / 3084480.0,
70 0.0,
71 0.0,
72 0.0,
73 0.0,
74 }, {
75 -855863.0 / 960120.0,
76 0.0,
77 0.0,
78 0.0,
79 0.0,
80 3399327.0 / 381920.0,
81 133156864.0 / 35198415.0,
82 -1674902723.0 / 144358200.0,
83 54980371265625.0 / 88346187905696.0,
84 -734375.0 / 2413152.0,
85 171414593.0 / 425630700.0,
86 137909.0 / 1542240.0,
87 -1.0,
88 0.0,
89 0.0,
90 0.0
91 }, {
92 -17751989329.0 / 2106076560.0,
93 0.0,
94 0.0,
95 0.0,
96 0.0,
97 4272954039.0 / 7539864640.0,
98 -118476319744.0 / 38604839385.0,
99 755123450731.0 / 316657731600.0,
100 3692384461234828125.0 / 1744130441634250432.0,
101 -4612609375.0 / 5293382976.0,
102 2091772278379.0 / 933644586600.0,
103 2136624137.0 / 3382989120.0,
104 -126493.0 / 1421424.0,
105 98350000.0 / 5419179.0,
106 -18878125.0 / 2053168.0,
107 -1944542619.0 / 438351368.0
108 }, {
109 32941697297.0 / 3159114840.0,
110 0.0,
111 0.0,
112 0.0,
113 0.0,
114 456696183123.0 / 1884966160.0,
115 19132610714624.0 / 115814518155.0,
116 -177904688592943.0 / 474986597400.0,
117 -4821139941836765625.0 / 218016305204281304.0,
118 30702015625.0 / 3970037232.0,
119 -85916079474274.0 / 2800933759800.0,
120 -5919468007.0 / 634310460.0,
121 2479159.0 / 157936.0,
122 -18750000.0 / 602131.0,
123 -19203125.0 / 2053168.0,
124 15700361463.0 / 438351368.0
125 }, {
126 12627015655.0 / 631822968.0,
127 0.0,
128 0.0,
129 0.0,
130 0.0,
131 -72955222965.0 / 188496616.0,
132 -13145744952320.0 / 69488710893.0,
133 30084216194513.0 / 56998391688.0,
134 -296858761006640625.0 / 25648977082856624.0,
135 569140625.0 / 82709109.0,
136 -18684190637.0 / 18672891732.0,
137 69644045.0 / 89549712.0,
138 -11847025.0 / 4264272.0,
139 -978650000.0 / 16257537.0,
140 519371875.0 / 6159504.0,
141 5256837225.0 / 438351368.0
142 }, {
143 -450944925.0 / 17550638.0,
144 0.0,
145 0.0,
146 0.0,
147 0.0,
148 -14532122925.0 / 94248308.0,
149 -595876966400.0 / 2573655959.0,
150 188748653015.0 / 527762886.0,
151 2545485458115234375.0 / 27252038150535163.0,
152 -1376953125.0 / 36759604.0,
153 53995596795.0 / 518691437.0,
154 210311225.0 / 7047894.0,
155 -1718875.0 / 39484.0,
156 58000000.0 / 602131.0,
157 -1546875.0 / 39484.0,
158 -1262172375.0 / 8429834.0
159 }
160 };
161
162
163
164
165
166
167
168
169
170
171 DormandPrince853StateInterpolator(final boolean forward,
172 final double[][] yDotK,
173 final ODEStateAndDerivative globalPreviousState,
174 final ODEStateAndDerivative globalCurrentState,
175 final ODEStateAndDerivative softPreviousState,
176 final ODEStateAndDerivative softCurrentState,
177 final EquationsMapper mapper) {
178 super(forward, yDotK,
179 globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
180 mapper);
181 }
182
183
184 @Override
185 protected DormandPrince853StateInterpolator create(final boolean newForward, final double[][] newYDotK,
186 final ODEStateAndDerivative newGlobalPreviousState,
187 final ODEStateAndDerivative newGlobalCurrentState,
188 final ODEStateAndDerivative newSoftPreviousState,
189 final ODEStateAndDerivative newSoftCurrentState,
190 final EquationsMapper newMapper) {
191 return new DormandPrince853StateInterpolator(newForward, newYDotK,
192 newGlobalPreviousState, newGlobalCurrentState,
193 newSoftPreviousState, newSoftCurrentState,
194 newMapper);
195 }
196
197
198 @Override
199 protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
200 final double time, final double theta,
201 final double thetaH, final double oneMinusThetaH) {
202
203 final double eta = 1.0 - theta;
204 final double twoTheta = 2 * theta;
205 final double theta2 = theta * theta;
206 final double dot1 = 1.0 - twoTheta;
207 final double dot2 = theta * (2 - 3 * theta);
208 final double dot3 = twoTheta * (theta * (twoTheta - 3) + 1);
209 final double dot4 = theta2 * (theta * (5 * theta - 8) + 3);
210 final double dot5 = theta2 * (theta * (theta * (15 - 6 * theta) - 12) + 3);
211 final double dot6 = theta2 * (theta * (theta * (theta * (18 - 7 * theta) - 15) + 4));
212 final double[] interpolatedState;
213 final double[] interpolatedDerivatives;
214
215
216 if (getGlobalPreviousState() != null && theta <= 0.5) {
217 final double f0 = thetaH;
218 final double f1 = f0 * eta;
219 final double f2 = f1 * theta;
220 final double f3 = f2 * eta;
221 final double f4 = f3 * theta;
222 final double f5 = f4 * eta;
223 final double f6 = f5 * theta;
224 final double[] p = new double[16];
225 final double[] q = new double[16];
226 for (int i = 0; i < p.length; ++i) {
227 p[i] = f0 * D[0][i] + f1 * D[1][i] + f2 * D[2][i] + f3 * D[3][i] +
228 f4 * D[4][i] + f5 * D[5][i] + f6 * D[6][i];
229 q[i] = D[0][i] + dot1 * D[1][i] + dot2 * D[2][i] + dot3 * D[3][i] +
230 dot4 * D[4][i] + dot5 * D[5][i] + dot6 * D[6][i];
231 }
232 interpolatedState = previousStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
233 p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
234 interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
235 q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
236 } else {
237 final double f0 = -oneMinusThetaH;
238 final double f1 = -f0 * theta;
239 final double f2 = f1 * theta;
240 final double f3 = f2 * eta;
241 final double f4 = f3 * theta;
242 final double f5 = f4 * eta;
243 final double f6 = f5 * theta;
244 final double[] p = new double[16];
245 final double[] q = new double[16];
246 for (int i = 0; i < p.length; ++i) {
247 p[i] = f0 * D[0][i] + f1 * D[1][i] + f2 * D[2][i] + f3 * D[3][i] +
248 f4 * D[4][i] + f5 * D[5][i] + f6 * D[6][i];
249 q[i] = D[0][i] + dot1 * D[1][i] + dot2 * D[2][i] + dot3 * D[3][i] +
250 dot4 * D[4][i] + dot5 * D[5][i] + dot6 * D[6][i];
251 }
252 interpolatedState = currentStateLinearCombination(p[0], p[1], p[ 2], p[ 3], p[ 4], p[ 5], p[ 6], p[ 7],
253 p[8], p[9], p[10], p[11], p[12], p[13], p[14], p[15]);
254 interpolatedDerivatives = derivativeLinearCombination(q[0], q[1], q[ 2], q[ 3], q[ 4], q[ 5], q[ 6], q[ 7],
255 q[8], q[9], q[10], q[11], q[12], q[13], q[14], q[15]);
256 }
257
258 return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);
259
260 }
261
262 }