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;
19
20 import java.util.ArrayList;
21 import java.util.Collection;
22 import java.util.Collections;
23 import java.util.Comparator;
24 import java.util.Iterator;
25 import java.util.List;
26 import java.util.SortedSet;
27 import java.util.TreeSet;
28
29 import org.apache.commons.math4.legacy.core.Field;
30 import org.apache.commons.math4.legacy.core.RealFieldElement;
31 import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver;
32 import org.apache.commons.math4.legacy.analysis.solvers.FieldBracketingNthOrderBrentSolver;
33 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
34 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
35 import org.apache.commons.math4.legacy.exception.NoBracketingException;
36 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
37 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
38 import org.apache.commons.math4.legacy.ode.events.FieldEventHandler;
39 import org.apache.commons.math4.legacy.ode.events.FieldEventState;
40 import org.apache.commons.math4.legacy.ode.sampling.AbstractFieldStepInterpolator;
41 import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
42 import org.apache.commons.math4.core.jdkmath.JdkMath;
43 import org.apache.commons.math4.legacy.core.IntegerSequence;
44
45
46
47
48
49
50 public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FirstOrderFieldIntegrator<T> {
51
52
53 private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
54
55
56 private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
57
58
59 private Collection<FieldStepHandler<T>> stepHandlers;
60
61
62 private FieldODEStateAndDerivative<T> stepStart;
63
64
65 private T stepSize;
66
67
68 private boolean isLastStep;
69
70
71 private boolean resetOccurred;
72
73
74 private final Field<T> field;
75
76
77 private Collection<FieldEventState<T>> eventsStates;
78
79
80 private boolean statesInitialized;
81
82
83 private final String name;
84
85
86 private IntegerSequence.Incrementor evaluations;
87
88
89 private transient FieldExpandableODE<T> equations;
90
91
92
93
94
95 protected AbstractFieldIntegrator(final Field<T> field, final String name) {
96 this.field = field;
97 this.name = name;
98 stepHandlers = new ArrayList<>();
99 stepStart = null;
100 stepSize = null;
101 eventsStates = new ArrayList<>();
102 statesInitialized = false;
103 evaluations = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
104 }
105
106
107
108
109 public Field<T> getField() {
110 return field;
111 }
112
113
114 @Override
115 public String getName() {
116 return name;
117 }
118
119
120 @Override
121 public void addStepHandler(final FieldStepHandler<T> handler) {
122 stepHandlers.add(handler);
123 }
124
125
126 @Override
127 public Collection<FieldStepHandler<T>> getStepHandlers() {
128 return Collections.unmodifiableCollection(stepHandlers);
129 }
130
131
132 @Override
133 public void clearStepHandlers() {
134 stepHandlers.clear();
135 }
136
137
138 @Override
139 public void addEventHandler(final FieldEventHandler<T> handler,
140 final double maxCheckInterval,
141 final double convergence,
142 final int maxIterationCount) {
143 addEventHandler(handler, maxCheckInterval, convergence,
144 maxIterationCount,
145 new FieldBracketingNthOrderBrentSolver<>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY),
146 field.getZero().add(convergence),
147 field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY),
148 5));
149 }
150
151
152 @Override
153 public void addEventHandler(final FieldEventHandler<T> handler,
154 final double maxCheckInterval,
155 final double convergence,
156 final int maxIterationCount,
157 final BracketedRealFieldUnivariateSolver<T> solver) {
158 eventsStates.add(new FieldEventState<>(handler, maxCheckInterval, field.getZero().add(convergence),
159 maxIterationCount, solver));
160 }
161
162
163 @Override
164 public Collection<FieldEventHandler<T>> getEventHandlers() {
165 final List<FieldEventHandler<T>> list = new ArrayList<>(eventsStates.size());
166 for (FieldEventState<T> state : eventsStates) {
167 list.add(state.getEventHandler());
168 }
169 return Collections.unmodifiableCollection(list);
170 }
171
172
173 @Override
174 public void clearEventHandlers() {
175 eventsStates.clear();
176 }
177
178
179 @Override
180 public FieldODEStateAndDerivative<T> getCurrentStepStart() {
181 return stepStart;
182 }
183
184
185 @Override
186 public T getCurrentSignedStepsize() {
187 return stepSize;
188 }
189
190
191 @Override
192 public void setMaxEvaluations(int maxEvaluations) {
193 evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
194 }
195
196
197 @Override
198 public int getMaxEvaluations() {
199 return evaluations.getMaximalCount();
200 }
201
202
203 @Override
204 public int getEvaluations() {
205 return evaluations.getCount();
206 }
207
208
209
210
211
212
213
214
215 protected FieldODEStateAndDerivative<T> initIntegration(final FieldExpandableODE<T> eqn,
216 final T t0, final T[] y0, final T t) {
217
218 this.equations = eqn;
219 evaluations = evaluations.withStart(0);
220
221
222 eqn.init(t0, y0, t);
223
224
225 final T[] y0Dot = computeDerivatives(t0, y0);
226 final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<>(t0, y0, y0Dot);
227
228
229 for (final FieldEventState<T> state : eventsStates) {
230 state.getEventHandler().init(state0, t);
231 }
232
233
234 for (FieldStepHandler<T> handler : stepHandlers) {
235 handler.init(state0, t);
236 }
237
238 setStateInitialized(false);
239
240 return state0;
241 }
242
243
244
245
246 protected FieldExpandableODE<T> getEquations() {
247 return equations;
248 }
249
250
251
252
253 protected IntegerSequence.Incrementor getEvaluationsCounter() {
254 return evaluations;
255 }
256
257
258
259
260
261
262
263
264
265
266
267 public T[] computeDerivatives(final T t, final T[] y)
268 throws DimensionMismatchException, MaxCountExceededException, NullPointerException {
269 evaluations.increment();
270 return equations.computeDerivatives(t, y);
271 }
272
273
274
275
276
277
278
279 protected void setStateInitialized(final boolean stateInitialized) {
280 this.statesInitialized = stateInitialized;
281 }
282
283
284
285
286
287
288
289
290
291
292 protected FieldODEStateAndDerivative<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator,
293 final T tEnd)
294 throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
295
296 FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState();
297 final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState();
298
299
300 if (! statesInitialized) {
301 for (FieldEventState<T> state : eventsStates) {
302 state.reinitializeBegin(interpolator);
303 }
304 statesInitialized = true;
305 }
306
307
308 final int orderingSign = interpolator.isForward() ? +1 : -1;
309 SortedSet<FieldEventState<T>> occurringEvents = new TreeSet<>(new Comparator<FieldEventState<T>>() {
310
311
312 @Override
313 public int compare(FieldEventState<T> es0, FieldEventState<T> es1) {
314 return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal());
315 }
316 });
317
318 for (final FieldEventState<T> state : eventsStates) {
319 if (state.evaluateStep(interpolator)) {
320
321 occurringEvents.add(state);
322 }
323 }
324
325 AbstractFieldStepInterpolator<T> restricted = interpolator;
326 while (!occurringEvents.isEmpty()) {
327
328
329 final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator();
330 final FieldEventState<T> currentEvent = iterator.next();
331 iterator.remove();
332
333
334 final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
335
336
337 restricted = restricted.restrictStep(previousState, eventState);
338
339
340 for (final FieldEventState<T> state : eventsStates) {
341 state.stepAccepted(eventState);
342 isLastStep = isLastStep || state.stop();
343 }
344
345
346 for (final FieldStepHandler<T> handler : stepHandlers) {
347 handler.handleStep(restricted, isLastStep);
348 }
349
350 if (isLastStep) {
351
352 return eventState;
353 }
354
355 FieldODEState<T> newState = null;
356 resetOccurred = false;
357 for (final FieldEventState<T> state : eventsStates) {
358 newState = state.reset(eventState);
359 if (newState != null) {
360
361
362 final T[] y = equations.getMapper().mapState(newState);
363 final T[] yDot = computeDerivatives(newState.getTime(), y);
364 resetOccurred = true;
365 return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot);
366 }
367 }
368
369
370 previousState = eventState;
371 restricted = restricted.restrictStep(eventState, currentState);
372
373
374 if (currentEvent.evaluateStep(restricted)) {
375
376 occurringEvents.add(currentEvent);
377 }
378 }
379
380
381 for (final FieldEventState<T> state : eventsStates) {
382 state.stepAccepted(currentState);
383 isLastStep = isLastStep || state.stop();
384 }
385 isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= JdkMath.ulp(tEnd.getReal());
386
387
388 for (FieldStepHandler<T> handler : stepHandlers) {
389 handler.handleStep(restricted, isLastStep);
390 }
391
392 return currentState;
393 }
394
395
396
397
398
399
400
401
402 protected void sanityChecks(final FieldODEState<T> eqn, final T t)
403 throws NumberIsTooSmallException, DimensionMismatchException {
404
405 final double threshold = 1000 * JdkMath.ulp(JdkMath.max(JdkMath.abs(eqn.getTime().getReal()),
406 JdkMath.abs(t.getReal())));
407 final double dt = eqn.getTime().subtract(t).abs().getReal();
408 if (dt <= threshold) {
409 throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
410 dt, threshold, false);
411 }
412 }
413
414
415
416
417 protected boolean resetOccurred() {
418 return resetOccurred;
419 }
420
421
422
423
424 protected void setStepSize(final T stepSize) {
425 this.stepSize = stepSize;
426 }
427
428
429
430
431 protected T getStepSize() {
432 return stepSize;
433 }
434
435
436
437 protected void setStepStart(final FieldODEStateAndDerivative<T> stepStart) {
438 this.stepStart = stepStart;
439 }
440
441
442
443
444 protected FieldODEStateAndDerivative<T> getStepStart() {
445 return stepStart;
446 }
447
448
449
450
451 protected void setIsLastStep(final boolean isLastStep) {
452 this.isLastStep = isLastStep;
453 }
454
455
456
457
458 protected boolean isLastStep() {
459 return isLastStep;
460 }
461 }