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.events;
24
25
26 import org.hipparchus.analysis.UnivariateFunction;
27 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
28 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
29 import org.hipparchus.exception.MathIllegalArgumentException;
30 import org.hipparchus.exception.MathIllegalStateException;
31 import org.hipparchus.ode.EquationsMapper;
32 import org.hipparchus.ode.ExpandableODE;
33 import org.hipparchus.ode.ODEState;
34 import org.hipparchus.ode.ODEStateAndDerivative;
35 import org.hipparchus.ode.OrdinaryDifferentialEquation;
36 import org.hipparchus.ode.SecondaryODE;
37 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
38 import org.hipparchus.ode.nonstiff.LutherIntegrator;
39 import org.hipparchus.ode.sampling.DummyStepInterpolator;
40 import org.hipparchus.ode.sampling.ODEStateInterpolator;
41 import org.junit.jupiter.api.Assertions;
42 import org.junit.jupiter.api.Test;
43 import org.junit.jupiter.params.ParameterizedTest;
44 import org.junit.jupiter.params.provider.ValueSource;
45 import org.mockito.Mockito;
46
47 import static org.junit.jupiter.api.Assertions.assertEquals;
48 import static org.junit.jupiter.api.Assertions.assertFalse;
49 import static org.junit.jupiter.api.Assertions.assertTrue;
50
51 class DetectorBasedEventStateTest {
52
53
54 @Test
55 void closeEvents() throws MathIllegalArgumentException, MathIllegalStateException {
56
57 final double r1 = 90.0;
58 final double r2 = 135.0;
59 final double gap = r2 - r1;
60
61 final double tolerance = 0.1;
62 DetectorBasedEventState es = new DetectorBasedEventState(new CloseEventsGenerator(r1, r2, 1.5 * gap, tolerance, 100));
63 EquationsMapper mapper = new ExpandableODE(new OrdinaryDifferentialEquation() {
64 @Override
65 public int getDimension() {
66 return 0;
67 }
68 @Override
69 public double[] computeDerivatives(double t, double[] y) {
70 return new double[0];
71 }
72 }).getMapper();
73
74 double[] a = new double[0];
75 ODEStateAndDerivative osdLongBefore = new ODEStateAndDerivative(r1 - 1.5 * gap, a, a);
76 ODEStateAndDerivative osBefore = new ODEStateAndDerivative(r1 - 0.5 * gap, a, a);
77 ODEStateInterpolator interpolatorA = new DummyStepInterpolator(true,
78 osdLongBefore, osBefore,
79 osdLongBefore, osBefore,
80 mapper);
81 es.reinitializeBegin(interpolatorA);
82 assertFalse(es.evaluateStep(interpolatorA));
83
84 ODEStateAndDerivative osdBetween = new ODEStateAndDerivative(0.5 * (r1 + r2), a, a);
85 ODEStateInterpolator interpolatorB = new DummyStepInterpolator(true,
86 osBefore, osdBetween,
87 osBefore, osdBetween,
88 mapper);
89 assertTrue(es.evaluateStep(interpolatorB));
90 assertEquals(r1, es.getEventTime(), tolerance);
91 ODEStateAndDerivative osdAtEvent = new ODEStateAndDerivative(es.getEventTime(), a, a);
92 es.doEvent(osdAtEvent);
93
94 ODEStateAndDerivative osdAfterSecond = new ODEStateAndDerivative(r2 + 0.4 * gap, a, a);
95 ODEStateInterpolator interpolatorC = new DummyStepInterpolator(true,
96 osdAtEvent, osdAfterSecond,
97 osdAtEvent, osdAfterSecond,
98 mapper);
99 assertTrue(es.evaluateStep(interpolatorC));
100 assertEquals(r2, es.getEventTime(), tolerance);
101
102 }
103
104
105 @Test
106 void testIssue695()
107 throws MathIllegalArgumentException, MathIllegalStateException {
108
109 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
110 public int getDimension() {
111 return 1;
112 }
113 public double[] computeDerivatives(double t, double[] y) {
114 return new double[] { 1.0 };
115 }
116 };
117
118 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
119 integrator.addEventDetector(new ResettingEvent(10.99, 0.1, 1.0e-9, 1000));
120 integrator.addEventDetector(new ResettingEvent(11.01, 0.1, 1.0e-9, 1000));
121 integrator.setInitialStepSize(3.0);
122
123 double target = 30.0;
124 ODEStateAndDerivative finalState =
125 integrator.integrate(equation, new ODEState(0.0, new double[1]), target);
126 assertEquals(target, finalState.getTime(), 1.0e-10);
127 assertEquals(32.0, finalState.getPrimaryState()[0], 1.0e-10);
128
129 }
130
131 private static class ResettingEvent implements ODEEventDetector {
132
133 private static double lastTriggerTime = Double.NEGATIVE_INFINITY;
134 private final AdaptableInterval maxCheck;
135 private final int maxIter;
136 private final BracketingNthOrderBrentSolver solver;
137 private final double tEvent;
138
139 public ResettingEvent(final double tEvent,
140 final double maxCheck, final double threshold, final int maxIter) {
141 this.maxCheck = (s, isForward) -> maxCheck;
142 this.maxIter = maxIter;
143 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
144 this.tEvent = tEvent;
145 }
146
147 public AdaptableInterval getMaxCheckInterval() {
148 return maxCheck;
149 }
150
151 public int getMaxIterationCount() {
152 return maxIter;
153 }
154
155 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
156 return solver;
157 }
158
159 public double g(ODEStateAndDerivative s) {
160
161
162
163
164 assertTrue(s.getTime() >= lastTriggerTime,
165 "going backard in time! (" + s.getTime() + " < " + lastTriggerTime + ")");
166 return s.getTime() - tEvent;
167 }
168
169 public ODEEventHandler getHandler() {
170 return new ODEEventHandler() {
171 public Action eventOccurred(ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) {
172
173 lastTriggerTime = s.getTime();
174 return Action.RESET_STATE;
175 }
176
177 public ODEStateAndDerivative resetState(ODEEventDetector detector, ODEStateAndDerivative s) {
178 double[] y = s.getPrimaryState();
179 y[0] += 1.0;
180 return new ODEStateAndDerivative(s.getTime(), y, s.getPrimaryDerivative());
181 }
182 };
183 }
184
185 }
186
187
188 @Test
189 void testIssue965()
190 throws MathIllegalArgumentException, MathIllegalStateException {
191
192 ExpandableODE equation = new ExpandableODE(new OrdinaryDifferentialEquation() {
193 public int getDimension() {
194 return 1;
195 }
196 public double[] computeDerivatives(double t, double[] y) {
197 return new double[] { 2.0 };
198 }
199 });
200 int index = equation.addSecondaryEquations(new SecondaryODE() {
201 public int getDimension() {
202 return 1;
203 }
204 public double[] computeDerivatives(double t, double[] primary,
205 double[] primaryDot, double[] secondary) {
206 return new double[] { -3.0 };
207 }
208 });
209 assertEquals(1, index);
210
211 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
212 integrator.addEventDetector(new SecondaryStateEvent(index, -3.0, 0.1, 1.0e-9, 1000));
213 integrator.setInitialStepSize(3.0);
214
215 ODEState initialState = new ODEState(0.0,
216 new double[] { 0.0 },
217 new double[][] { { 0.0 } });
218 ODEStateAndDerivative finalState = integrator.integrate(equation, initialState, 30.0);
219 assertEquals( 1.0, finalState.getTime(), 1.0e-10);
220 assertEquals( 2.0, finalState.getPrimaryState()[0], 1.0e-10);
221 assertEquals(-3.0, finalState.getSecondaryState(index)[0], 1.0e-10);
222
223 }
224
225 private static class SecondaryStateEvent implements ODEEventDetector {
226
227 private final AdaptableInterval maxCheck;
228 private final int maxIter;
229 private final BracketingNthOrderBrentSolver solver;
230 private int index;
231 private final double target;
232
233 public SecondaryStateEvent(final int index, final double target,
234 final double maxCheck, final double threshold, final int maxIter) {
235 this.maxCheck = (s, isForward) -> maxCheck;
236 this.maxIter = maxIter;
237 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
238 this.index = index;
239 this.target = target;
240 }
241
242 public AdaptableInterval getMaxCheckInterval() {
243 return maxCheck;
244 }
245
246 public int getMaxIterationCount() {
247 return maxIter;
248 }
249
250 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
251 return solver;
252 }
253
254
255 public ODEEventHandler getHandler() {
256 return (state, detector, increasing) -> Action.STOP;
257 }
258
259 public double g(ODEStateAndDerivative s) {
260 return s.getSecondaryState(index)[0] - target;
261 }
262
263 }
264
265 @Test
266 void testEventsCloserThanThreshold()
267 throws MathIllegalArgumentException, MathIllegalStateException {
268
269 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
270
271 public int getDimension() {
272 return 1;
273 }
274
275 public double[] computeDerivatives(double t, double[] y) {
276 return new double[] { 1.0 };
277 }
278 };
279
280 LutherIntegrator integrator = new LutherIntegrator(20.0);
281 CloseEventsGenerator eventsGenerator =
282 new CloseEventsGenerator(9.0 - 1.0 / 128, 9.0 + 1.0 / 128, 1.0, 0.02, 1000);
283 integrator.addEventDetector(eventsGenerator);
284 double tEnd = integrator.integrate(equation, new ODEState(0.0, new double[1]), 100.0).getTime();
285 assertEquals( 2, eventsGenerator.getCount());
286 assertEquals( 9.0 + 1.0 / 128, tEnd, 1.0 / 32.0);
287
288 }
289
290 private class CloseEventsGenerator implements ODEEventDetector {
291
292 private final AdaptableInterval maxCheck;
293 private final int maxIter;
294 private final BracketingNthOrderBrentSolver solver;
295 final double r1;
296 final double r2;
297 int count;
298
299 public CloseEventsGenerator(final double r1, final double r2,
300 final double maxCheck, final double threshold, final int maxIter) {
301 this.maxCheck = (s, isForward) -> maxCheck;
302 this.maxIter = maxIter;
303 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
304 this.r1 = r1;
305 this.r2 = r2;
306 this.count = 0;
307 }
308
309 public AdaptableInterval getMaxCheckInterval() {
310 return maxCheck;
311 }
312
313 public int getMaxIterationCount() {
314 return maxIter;
315 }
316
317 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
318 return solver;
319 }
320
321 public double g(ODEStateAndDerivative s) {
322 return (s.getTime() - r1) * (r2 - s.getTime());
323 }
324
325 public ODEEventHandler getHandler() {
326 return (state, detector, increasing) -> ++count < 2 ? Action.CONTINUE : Action.STOP;
327 }
328
329 public int getCount() {
330 return count;
331 }
332
333 }
334
335 @ParameterizedTest
336 @ValueSource(booleans = {true, false})
337 void testAdaptableInterval(final boolean isForward) {
338
339 final TestDetector detector = new TestDetector();
340 final DetectorBasedEventState eventState = new DetectorBasedEventState(detector);
341 final ODEStateInterpolator mockedInterpolator = Mockito.mock(ODEStateInterpolator.class);
342 final ODEStateAndDerivative stateAndDerivative1 = getStateAndDerivative(1);
343 final ODEStateAndDerivative stateAndDerivative2 = getStateAndDerivative(-1);
344 if (isForward) {
345 Mockito.when(mockedInterpolator.getCurrentState()).thenReturn(stateAndDerivative1);
346 Mockito.when(mockedInterpolator.getPreviousState()).thenReturn(stateAndDerivative2);
347 } else {
348 Mockito.when(mockedInterpolator.getCurrentState()).thenReturn(stateAndDerivative2);
349 Mockito.when(mockedInterpolator.getPreviousState()).thenReturn(stateAndDerivative1);
350 }
351 Mockito.when(mockedInterpolator.isForward()).thenReturn(isForward);
352
353 eventState.evaluateStep(mockedInterpolator);
354
355 if (isForward) {
356 Assertions.assertEquals(1, detector.triggeredForward);
357 Assertions.assertEquals(0, detector.triggeredBackward);
358 } else {
359 Assertions.assertEquals(0, detector.triggeredForward);
360 Assertions.assertEquals(1, detector.triggeredBackward);
361 }
362 }
363
364 private static ODEStateAndDerivative getStateAndDerivative(final double time) {
365 return new ODEStateAndDerivative(time, new double[] {time}, new double[1]);
366 }
367
368 private static class TestDetector implements ODEEventDetector {
369
370 int triggeredForward = 0;
371 int triggeredBackward = 0;
372
373 @Override
374 public AdaptableInterval getMaxCheckInterval() {
375 return (state, isForward) -> {
376 if (isForward) {
377 triggeredForward++;
378 } else {
379 triggeredBackward++;
380 }
381 return 1.;
382 };
383 }
384
385 @Override
386 public int getMaxIterationCount() {
387 return 10;
388 }
389
390 @Override
391 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
392 return new BracketingNthOrderBrentSolver();
393 }
394
395 @Override
396 public ODEEventHandler getHandler() {
397 return null;
398 }
399
400 @Override
401 public double g(ODEStateAndDerivative state) {
402 return 0.;
403 }
404 }
405
406 }