1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.ode.nonstiff;
24
25 import org.hipparchus.CalculusFieldElement;
26 import org.hipparchus.Field;
27 import org.hipparchus.ode.FieldEquationsMapper;
28 import org.hipparchus.ode.FieldODEStateAndDerivative;
29
30
31
32
33
34
35
36
37
38
39 class DormandPrince54FieldStateInterpolator<T extends CalculusFieldElement<T>>
40 extends RungeKuttaFieldStateInterpolator<T> {
41
42
43 private final T a70;
44
45
46
47
48 private final T a72;
49
50
51 private final T a73;
52
53
54 private final T a74;
55
56
57 private final T a75;
58
59
60 private final T d0;
61
62
63
64
65 private final T d2;
66
67
68 private final T d3;
69
70
71 private final T d4;
72
73
74 private final T d5;
75
76
77 private final T d6;
78
79
80
81
82
83
84
85
86
87
88
89 DormandPrince54FieldStateInterpolator(final Field<T> field, final boolean forward,
90 final T[][] yDotK,
91 final FieldODEStateAndDerivative<T> globalPreviousState,
92 final FieldODEStateAndDerivative<T> globalCurrentState,
93 final FieldODEStateAndDerivative<T> softPreviousState,
94 final FieldODEStateAndDerivative<T> softCurrentState,
95 final FieldEquationsMapper<T> mapper) {
96 super(field, forward, yDotK,
97 globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
98 mapper);
99 final T one = field.getOne();
100 a70 = one.newInstance( 35.0 / 384.0);
101 a72 = one.newInstance( 500.0 / 1113.0);
102 a73 = one.newInstance( 125.0 / 192.0);
103 a74 = one.newInstance(-2187.0 / 6784.0);
104 a75 = one.newInstance( 11.0 / 84.0);
105 d0 = one.newInstance(-12715105075.0 / 11282082432.0);
106 d2 = one.newInstance( 87487479700.0 / 32700410799.0);
107 d3 = one.newInstance(-10690763975.0 / 1880347072.0);
108 d4 = one.newInstance(701980252875.0 / 199316789632.0);
109 d5 = one.newInstance( -1453857185.0 / 822651844.0);
110 d6 = one.newInstance( 69997945.0 / 29380423.0);
111 }
112
113
114 @Override
115 protected DormandPrince54FieldStateInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
116 final FieldODEStateAndDerivative<T> newGlobalPreviousState,
117 final FieldODEStateAndDerivative<T> newGlobalCurrentState,
118 final FieldODEStateAndDerivative<T> newSoftPreviousState,
119 final FieldODEStateAndDerivative<T> newSoftCurrentState,
120 final FieldEquationsMapper<T> newMapper) {
121 return new DormandPrince54FieldStateInterpolator<T>(newField, newForward, newYDotK,
122 newGlobalPreviousState, newGlobalCurrentState,
123 newSoftPreviousState, newSoftCurrentState,
124 newMapper);
125 }
126
127 @SuppressWarnings("unchecked")
128 @Override
129 protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
130 final T time, final T theta,
131 final T thetaH, final T oneMinusThetaH) {
132
133
134 final T one = time.getField().getOne();
135 final T eta = one.subtract(theta);
136 final T twoTheta = theta.multiply(2);
137 final T dot2 = one.subtract(twoTheta);
138 final T dot3 = theta.multiply(theta.multiply(-3).add(2));
139 final T dot4 = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
140 final T[] interpolatedState;
141 final T[] interpolatedDerivatives;
142 if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
143 final T f1 = thetaH;
144 final T f2 = f1.multiply(eta);
145 final T f3 = f2.multiply(theta);
146 final T f4 = f3.multiply(eta);
147 final T coeff0 = f1.multiply(a70).
148 subtract(f2.multiply(a70.subtract(1))).
149 add(f3.multiply(a70.multiply(2).subtract(1))).
150 add(f4.multiply(d0));
151 final T coeff1 = time.getField().getZero();
152 final T coeff2 = f1.multiply(a72).
153 subtract(f2.multiply(a72)).
154 add(f3.multiply(a72.multiply(2))).
155 add(f4.multiply(d2));
156 final T coeff3 = f1.multiply(a73).
157 subtract(f2.multiply(a73)).
158 add(f3.multiply(a73.multiply(2))).
159 add(f4.multiply(d3));
160 final T coeff4 = f1.multiply(a74).
161 subtract(f2.multiply(a74)).
162 add(f3.multiply(a74.multiply(2))).
163 add(f4.multiply(d4));
164 final T coeff5 = f1.multiply(a75).
165 subtract(f2.multiply(a75)).
166 add(f3.multiply(a75.multiply(2))).
167 add(f4.multiply(d5));
168 final T coeff6 = f4.multiply(d6).subtract(f3);
169 final T coeffDot0 = a70.
170 subtract(dot2.multiply(a70.subtract(1))).
171 add(dot3.multiply(a70.multiply(2).subtract(1))).
172 add(dot4.multiply(d0));
173 final T coeffDot1 = time.getField().getZero();
174 final T coeffDot2 = a72.
175 subtract(dot2.multiply(a72)).
176 add(dot3.multiply(a72.multiply(2))).
177 add(dot4.multiply(d2));
178 final T coeffDot3 = a73.
179 subtract(dot2.multiply(a73)).
180 add(dot3.multiply(a73.multiply(2))).
181 add(dot4.multiply(d3));
182 final T coeffDot4 = a74.
183 subtract(dot2.multiply(a74)).
184 add(dot3.multiply(a74.multiply(2))).
185 add(dot4.multiply(d4));
186 final T coeffDot5 = a75.
187 subtract(dot2.multiply(a75)).
188 add(dot3.multiply(a75.multiply(2))).
189 add(dot4.multiply(d5));
190 final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
191 interpolatedState = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
192 coeff4, coeff5, coeff6);
193 interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
194 coeffDot4, coeffDot5, coeffDot6);
195 } else {
196 final T f1 = oneMinusThetaH.negate();
197 final T f2 = oneMinusThetaH.multiply(theta);
198 final T f3 = f2.multiply(theta);
199 final T f4 = f3.multiply(eta);
200 final T coeff0 = f1.multiply(a70).
201 subtract(f2.multiply(a70.subtract(1))).
202 add(f3.multiply(a70.multiply(2).subtract(1))).
203 add(f4.multiply(d0));
204 final T coeff1 = time.getField().getZero();
205 final T coeff2 = f1.multiply(a72).
206 subtract(f2.multiply(a72)).
207 add(f3.multiply(a72.multiply(2))).
208 add(f4.multiply(d2));
209 final T coeff3 = f1.multiply(a73).
210 subtract(f2.multiply(a73)).
211 add(f3.multiply(a73.multiply(2))).
212 add(f4.multiply(d3));
213 final T coeff4 = f1.multiply(a74).
214 subtract(f2.multiply(a74)).
215 add(f3.multiply(a74.multiply(2))).
216 add(f4.multiply(d4));
217 final T coeff5 = f1.multiply(a75).
218 subtract(f2.multiply(a75)).
219 add(f3.multiply(a75.multiply(2))).
220 add(f4.multiply(d5));
221 final T coeff6 = f4.multiply(d6).subtract(f3);
222 final T coeffDot0 = a70.
223 subtract(dot2.multiply(a70.subtract(1))).
224 add(dot3.multiply(a70.multiply(2).subtract(1))).
225 add(dot4.multiply(d0));
226 final T coeffDot1 = time.getField().getZero();
227 final T coeffDot2 = a72.
228 subtract(dot2.multiply(a72)).
229 add(dot3.multiply(a72.multiply(2))).
230 add(dot4.multiply(d2));
231 final T coeffDot3 = a73.
232 subtract(dot2.multiply(a73)).
233 add(dot3.multiply(a73.multiply(2))).
234 add(dot4.multiply(d3));
235 final T coeffDot4 = a74.
236 subtract(dot2.multiply(a74)).
237 add(dot3.multiply(a74.multiply(2))).
238 add(dot4.multiply(d4));
239 final T coeffDot5 = a75.
240 subtract(dot2.multiply(a75)).
241 add(dot3.multiply(a75.multiply(2))).
242 add(dot4.multiply(d5));
243 final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
244 interpolatedState = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
245 coeff4, coeff5, coeff6);
246 interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
247 coeffDot4, coeffDot5, coeffDot6);
248 }
249 return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);
250
251 }
252
253 }