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