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.events;
19
20 import org.apache.commons.math4.legacy.core.RealFieldElement;
21 import org.apache.commons.math4.legacy.analysis.RealFieldUnivariateFunction;
22 import org.apache.commons.math4.legacy.analysis.solvers.AllowedSolution;
23 import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver;
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.exception.NoBracketingException;
26 import org.apache.commons.math4.legacy.ode.FieldODEState;
27 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
28 import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
29 import org.apache.commons.math4.core.jdkmath.JdkMath;
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44 public class FieldEventState<T extends RealFieldElement<T>> {
45
46
47 private final FieldEventHandler<T> handler;
48
49
50 private final double maxCheckInterval;
51
52
53 private final T convergence;
54
55
56 private final int maxIterationCount;
57
58
59 private T t0;
60
61
62 private T g0;
63
64
65 private boolean g0Positive;
66
67
68 private boolean pendingEvent;
69
70
71 private T pendingEventTime;
72
73
74 private T previousEventTime;
75
76
77 private boolean forward;
78
79
80
81
82 private boolean increasing;
83
84
85 private Action nextAction;
86
87
88 private final BracketedRealFieldUnivariateSolver<T> solver;
89
90
91
92
93
94
95
96
97
98
99
100 public FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval,
101 final T convergence, final int maxIterationCount,
102 final BracketedRealFieldUnivariateSolver<T> solver) {
103 this.handler = handler;
104 this.maxCheckInterval = maxCheckInterval;
105 this.convergence = convergence.abs();
106 this.maxIterationCount = maxIterationCount;
107 this.solver = solver;
108
109
110 t0 = null;
111 g0 = null;
112 g0Positive = true;
113 pendingEvent = false;
114 pendingEventTime = null;
115 previousEventTime = null;
116 increasing = true;
117 nextAction = Action.CONTINUE;
118 }
119
120
121
122
123 public FieldEventHandler<T> getEventHandler() {
124 return handler;
125 }
126
127
128
129
130 public double getMaxCheckInterval() {
131 return maxCheckInterval;
132 }
133
134
135
136
137 public T getConvergence() {
138 return convergence;
139 }
140
141
142
143
144 public int getMaxIterationCount() {
145 return maxIterationCount;
146 }
147
148
149
150
151
152
153 public void reinitializeBegin(final FieldStepInterpolator<T> interpolator)
154 throws MaxCountExceededException {
155
156 final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
157 t0 = s0.getTime();
158 g0 = handler.g(s0);
159 if (g0.getReal() == 0) {
160
161
162
163
164
165
166
167
168
169
170
171
172
173 final double epsilon = JdkMath.max(solver.getAbsoluteAccuracy().getReal(),
174 JdkMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
175 final T tStart = t0.add(0.5 * epsilon);
176 g0 = handler.g(interpolator.getInterpolatedState(tStart));
177 }
178 g0Positive = g0.getReal() >= 0;
179 }
180
181
182
183
184
185
186
187
188
189 public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
190 throws MaxCountExceededException, NoBracketingException {
191
192 forward = interpolator.isForward();
193 final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
194 final T t1 = s1.getTime();
195 final T dt = t1.subtract(t0);
196 if (dt.abs().subtract(convergence).getReal() < 0) {
197
198 return false;
199 }
200 final int n = JdkMath.max(1, (int) JdkMath.ceil(JdkMath.abs(dt.getReal()) / maxCheckInterval));
201 final T h = dt.divide(n);
202
203 final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() {
204
205 @Override
206 public T value(final T t) {
207 return handler.g(interpolator.getInterpolatedState(t));
208 }
209 };
210
211 T ta = t0;
212 T ga = g0;
213 for (int i = 0; i < n; ++i) {
214
215
216 final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
217 final T gb = handler.g(interpolator.getInterpolatedState(tb));
218
219
220 if (g0Positive ^ (gb.getReal() >= 0)) {
221
222
223
224 increasing = gb.subtract(ga).getReal() >= 0;
225
226
227 final T root = forward ?
228 solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
229 solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
230
231 if (previousEventTime != null &&
232 root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
233 root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
234
235
236
237
238 do {
239 ta = forward ? ta.add(convergence) : ta.subtract(convergence);
240 ga = f.value(ta);
241 } while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
242
243 if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
244
245 --i;
246 } else {
247
248
249
250 pendingEventTime = root;
251 pendingEvent = true;
252 return true;
253 }
254 } else if (previousEventTime == null ||
255 previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
256 pendingEventTime = root;
257 pendingEvent = true;
258 return true;
259 } else {
260
261 ta = tb;
262 ga = gb;
263 }
264 } else {
265
266 ta = tb;
267 ga = gb;
268 }
269 }
270
271
272 pendingEvent = false;
273 pendingEventTime = null;
274 return false;
275 }
276
277
278
279
280
281 public T getEventTime() {
282 return pendingEvent ?
283 pendingEventTime :
284 t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
285 }
286
287
288
289
290 public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
291
292 t0 = state.getTime();
293 g0 = handler.g(state);
294
295 if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
296
297 previousEventTime = state.getTime();
298 g0Positive = increasing;
299 nextAction = handler.eventOccurred(state, increasing == forward);
300 } else {
301 g0Positive = g0.getReal() >= 0;
302 nextAction = Action.CONTINUE;
303 }
304 }
305
306
307
308
309
310 public boolean stop() {
311 return nextAction == Action.STOP;
312 }
313
314
315
316
317
318
319 public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
320
321 if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
322 return null;
323 }
324
325 final FieldODEState<T> newState;
326 if (nextAction == Action.RESET_STATE) {
327 newState = handler.resetState(state);
328 } else if (nextAction == Action.RESET_DERIVATIVES) {
329 newState = state;
330 } else {
331 newState = null;
332 }
333 pendingEvent = false;
334 pendingEventTime = null;
335
336 return newState;
337 }
338 }