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 FieldExplicitRungeKuttaIntegrator<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 double[] realC = new double[0];
90
91
92 private double[][] realA = new double[0][];
93
94
95 private double[] realB = new double[0];
96
97
98 private final double exp;
99
100
101 private T safety;
102
103
104 private T minReduction;
105
106
107 private T maxGrowth;
108
109
110 private boolean usingFieldCoefficients;
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
127 final double minStep, final double maxStep,
128 final double scalAbsoluteTolerance,
129 final double scalRelativeTolerance) {
130
131 super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
132
133 this.fsal = fsal;
134 this.c = getC();
135 this.a = getA();
136 this.b = getB();
137
138 exp = -1.0 / getOrder();
139
140
141 setSafety(field.getZero().add(0.9));
142 setMinReduction(field.getZero().add(0.2));
143 setMaxGrowth(field.getZero().add(10.0));
144
145 }
146
147
148
149
150
151
152
153
154
155
156
157
158
159 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
160 final double minStep, final double maxStep,
161 final double[] vecAbsoluteTolerance,
162 final double[] vecRelativeTolerance) {
163
164 super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
165 this.usingFieldCoefficients = false;
166
167 this.fsal = fsal;
168 this.c = getC();
169 this.a = getA();
170 this.b = getB();
171
172 exp = -1.0 / getOrder();
173
174
175 setSafety(field.getZero().add(0.9));
176 setMinReduction(field.getZero().add(0.2));
177 setMaxGrowth(field.getZero().add(10.0));
178
179 }
180
181
182
183
184
185
186
187
188
189 protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
190 FieldODEStateAndDerivative<T> globalPreviousState,
191 FieldODEStateAndDerivative<T> globalCurrentState,
192 FieldEquationsMapper<T> mapper);
193
194
195
196
197 public abstract int getOrder();
198
199
200
201
202 public T getSafety() {
203 return safety;
204 }
205
206
207
208
209 public void setSafety(final T safety) {
210 this.safety = safety;
211 }
212
213
214
215
216
217
218 public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
219 this.usingFieldCoefficients = usingFieldCoefficients;
220 }
221
222
223 @Override
224 public boolean isUsingFieldCoefficients() {
225 return usingFieldCoefficients;
226 }
227
228
229 @Override
230 public int getNumberOfStages() {
231 return b.length;
232 }
233
234
235 @Override
236 protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
237 if (!isUsingFieldCoefficients()) {
238 realA = getRealA();
239 realB = getRealB();
240 realC = getRealC();
241 }
242 return super.initIntegration(eqn, s0, t);
243 }
244
245
246 @Override
247 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
248 final FieldODEState<T> initialState, final T finalTime)
249 throws MathIllegalArgumentException, MathIllegalStateException {
250
251 sanityChecks(initialState, finalTime);
252 setStepStart(initIntegration(equations, initialState, finalTime));
253 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
254
255
256 final int stages = getNumberOfStages();
257 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
258 T[] yTmp = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
259
260
261 T hNew = getField().getZero();
262 boolean firstTime = true;
263
264
265 setIsLastStep(false);
266 do {
267
268
269 double error = 10.0;
270 while (error >= 1.0) {
271
272
273 final T[] y = getStepStart().getCompleteState();
274 yDotK[0] = getStepStart().getCompleteDerivative();
275
276 if (firstTime) {
277 final StepsizeHelper helper = getStepSizeHelper();
278 final T[] scale = MathArrays.buildArray(getField(), helper.getMainSetDimension());
279 for (int i = 0; i < scale.length; ++i) {
280 scale[i] = helper.getTolerance(i, y[i].abs());
281 }
282 hNew = getField().getZero().add(initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper()));
283 firstTime = false;
284 }
285
286 setStepSize(hNew);
287 if (forward) {
288 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
289 setStepSize(finalTime.subtract(getStepStart().getTime()));
290 }
291 } else {
292 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
293 setStepSize(finalTime.subtract(getStepStart().getTime()));
294 }
295 }
296
297
298 if (isUsingFieldCoefficients()) {
299 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
300 getStepStart().getTime(), y, getStepSize(), a, c, yDotK);
301 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
302 } else {
303 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
304 getStepStart().getTime(), y, getStepSize(), realA, realC, yDotK);
305 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
306 }
307
308 incrementEvaluations(stages - 1);
309
310
311 error = estimateError(yDotK, y, yTmp, getStepSize());
312 if (error >= 1.0) {
313
314 final T factor = MathUtils.min(maxGrowth,
315 MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
316 hNew = getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false);
317 }
318
319 }
320 final T stepEnd = getStepStart().getTime().add(getStepSize());
321 final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
322 final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
323
324
325 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
326 finalTime));
327
328 if (!isLastStep()) {
329
330
331 final T factor = MathUtils.min(maxGrowth,
332 MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
333 final T scaledH = getStepSize().multiply(factor);
334 final T nextT = getStepStart().getTime().add(scaledH);
335 final boolean nextIsLast = forward ?
336 nextT.subtract(finalTime).getReal() >= 0 :
337 nextT.subtract(finalTime).getReal() <= 0;
338 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
339
340 final T filteredNextT = getStepStart().getTime().add(hNew);
341 final boolean filteredNextIsLast = forward ?
342 filteredNextT.subtract(finalTime).getReal() >= 0 :
343 filteredNextT.subtract(finalTime).getReal() <= 0;
344 if (filteredNextIsLast) {
345 hNew = finalTime.subtract(getStepStart().getTime());
346 }
347
348 }
349
350 } while (!isLastStep());
351
352 final FieldODEStateAndDerivative<T> finalState = getStepStart();
353 resetInternalState();
354 return finalState;
355
356 }
357
358
359
360
361 public T getMinReduction() {
362 return minReduction;
363 }
364
365
366
367
368 public void setMinReduction(final T minReduction) {
369 this.minReduction = minReduction;
370 }
371
372
373
374
375 public T getMaxGrowth() {
376 return maxGrowth;
377 }
378
379
380
381
382 public void setMaxGrowth(final T maxGrowth) {
383 this.maxGrowth = maxGrowth;
384 }
385
386
387
388
389
390
391
392
393 protected abstract double estimateError(T[][] yDotK, T[] y0, T[] y1, T h);
394
395 }