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.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.ode.FieldEquationsMapper;
30 import org.hipparchus.ode.FieldExpandableODE;
31 import org.hipparchus.ode.FieldODEState;
32 import org.hipparchus.ode.FieldODEStateAndDerivative;
33 import org.hipparchus.util.FastMath;
34 import org.hipparchus.util.MathArrays;
35 import org.hipparchus.util.MathUtils;
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72 public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
73 extends AdaptiveStepsizeFieldIntegrator<T>
74 implements FieldButcherArrayProvider<T> {
75
76
77 private final int fsal;
78
79
80 private final T[] c;
81
82
83 private final T[][] a;
84
85
86 private final T[] b;
87
88
89 private final double exp;
90
91
92 private T safety;
93
94
95 private T minReduction;
96
97
98 private T maxGrowth;
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
115 final double minStep, final double maxStep,
116 final double scalAbsoluteTolerance,
117 final double scalRelativeTolerance) {
118
119 super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
120
121 this.fsal = fsal;
122 this.c = getC();
123 this.a = getA();
124 this.b = getB();
125
126 exp = -1.0 / getOrder();
127
128
129 setSafety(field.getZero().add(0.9));
130 setMinReduction(field.getZero().add(0.2));
131 setMaxGrowth(field.getZero().add(10.0));
132
133 }
134
135
136
137
138
139
140
141
142
143
144
145
146
147 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
148 final double minStep, final double maxStep,
149 final double[] vecAbsoluteTolerance,
150 final double[] vecRelativeTolerance) {
151
152 super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
153
154 this.fsal = fsal;
155 this.c = getC();
156 this.a = getA();
157 this.b = getB();
158
159 exp = -1.0 / getOrder();
160
161
162 setSafety(field.getZero().add(0.9));
163 setMinReduction(field.getZero().add(0.2));
164 setMaxGrowth(field.getZero().add(10.0));
165
166 }
167
168
169
170
171
172
173 protected T fraction(final int p, final int q) {
174 return getField().getOne().multiply(p).divide(q);
175 }
176
177
178
179
180
181
182 protected T fraction(final double p, final double q) {
183 return getField().getOne().multiply(p).divide(q);
184 }
185
186
187
188
189
190
191
192
193
194 protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
195 FieldODEStateAndDerivative<T> globalPreviousState,
196 FieldODEStateAndDerivative<T> globalCurrentState,
197 FieldEquationsMapper<T> mapper);
198
199
200
201 public abstract int getOrder();
202
203
204
205
206 public T getSafety() {
207 return safety;
208 }
209
210
211
212
213 public void setSafety(final T safety) {
214 this.safety = safety;
215 }
216
217
218 @Override
219 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
220 final FieldODEState<T> initialState, final T finalTime)
221 throws MathIllegalArgumentException, MathIllegalStateException {
222
223 sanityChecks(initialState, finalTime);
224 setStepStart(initIntegration(equations, initialState, finalTime));
225 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
226
227
228 final int stages = c.length + 1;
229 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
230 final T[] yTmp = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
231
232
233 T hNew = getField().getZero();
234 boolean firstTime = true;
235
236
237 setIsLastStep(false);
238 do {
239
240
241 double error = 10.0;
242 while (error >= 1.0) {
243
244
245 final T[] y = getStepStart().getCompleteState();
246 yDotK[0] = getStepStart().getCompleteDerivative();
247
248 if (firstTime) {
249 final StepsizeHelper helper = getStepSizeHelper();
250 final T[] scale = MathArrays.buildArray(getField(), helper.getMainSetDimension());
251 for (int i = 0; i < scale.length; ++i) {
252 scale[i] = helper.getTolerance(i, y[i].abs());
253 }
254 hNew = getField().getZero().add(initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper()));
255 firstTime = false;
256 }
257
258 setStepSize(hNew);
259 if (forward) {
260 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
261 setStepSize(finalTime.subtract(getStepStart().getTime()));
262 }
263 } else {
264 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
265 setStepSize(finalTime.subtract(getStepStart().getTime()));
266 }
267 }
268
269
270 for (int k = 1; k < stages; ++k) {
271
272 for (int j = 0; j < y.length; ++j) {
273 T sum = yDotK[0][j].multiply(a[k-1][0]);
274 for (int l = 1; l < k; ++l) {
275 sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
276 }
277 yTmp[j] = y[j].add(getStepSize().multiply(sum));
278 }
279
280 yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp);
281
282 }
283
284
285 for (int j = 0; j < y.length; ++j) {
286 T sum = yDotK[0][j].multiply(b[0]);
287 for (int l = 1; l < stages; ++l) {
288 sum = sum.add(yDotK[l][j].multiply(b[l]));
289 }
290 yTmp[j] = y[j].add(getStepSize().multiply(sum));
291 }
292
293
294 error = estimateError(yDotK, y, yTmp, getStepSize());
295 if (error >= 1.0) {
296
297 final T factor = MathUtils.min(maxGrowth,
298 MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
299 hNew = getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false);
300 }
301
302 }
303 final T stepEnd = getStepStart().getTime().add(getStepSize());
304 final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
305 final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
306
307
308 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
309 finalTime));
310
311 if (!isLastStep()) {
312
313
314 final T factor = MathUtils.min(maxGrowth,
315 MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
316 final T scaledH = getStepSize().multiply(factor);
317 final T nextT = getStepStart().getTime().add(scaledH);
318 final boolean nextIsLast = forward ?
319 nextT.subtract(finalTime).getReal() >= 0 :
320 nextT.subtract(finalTime).getReal() <= 0;
321 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
322
323 final T filteredNextT = getStepStart().getTime().add(hNew);
324 final boolean filteredNextIsLast = forward ?
325 filteredNextT.subtract(finalTime).getReal() >= 0 :
326 filteredNextT.subtract(finalTime).getReal() <= 0;
327 if (filteredNextIsLast) {
328 hNew = finalTime.subtract(getStepStart().getTime());
329 }
330
331 }
332
333 } while (!isLastStep());
334
335 final FieldODEStateAndDerivative<T> finalState = getStepStart();
336 resetInternalState();
337 return finalState;
338
339 }
340
341
342
343
344 public T getMinReduction() {
345 return minReduction;
346 }
347
348
349
350
351 public void setMinReduction(final T minReduction) {
352 this.minReduction = minReduction;
353 }
354
355
356
357
358 public T getMaxGrowth() {
359 return maxGrowth;
360 }
361
362
363
364
365 public void setMaxGrowth(final T maxGrowth) {
366 this.maxGrowth = maxGrowth;
367 }
368
369
370
371
372
373
374
375
376 protected abstract double estimateError(T[][] yDotK, T[] y0, T[] y1, T h);
377
378 }