1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.hipparchus.ode.events;
18
19 import java.util.ArrayList;
20 import java.util.Arrays;
21 import java.util.List;
22
23 import org.hipparchus.analysis.UnivariateFunction;
24 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
25 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
26 import org.hipparchus.ode.ODEIntegrator;
27 import org.hipparchus.ode.ODEState;
28 import org.hipparchus.ode.ODEStateAndDerivative;
29 import org.hipparchus.ode.OrdinaryDifferentialEquation;
30 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
31 import org.hipparchus.ode.nonstiff.LutherIntegrator;
32 import org.hipparchus.ode.sampling.ODEStateInterpolator;
33 import org.hipparchus.ode.sampling.ODEStepHandler;
34 import org.hipparchus.util.FastMath;
35 import org.junit.Assert;
36 import org.junit.Test;
37
38
39
40
41
42
43 public class CloseEventsTest {
44
45 @Test
46 public void testCloseEventsFirstOneIsReset() {
47
48
49
50
51
52
53
54
55
56
57 ODEIntegrator integrator =
58 new DormandPrince853Integrator(10, 100.0, 1e-7, 1e-7);
59
60 TimeDetector detector1 = new TimeDetector(10, 1e-9, 100, Action.RESET_DERIVATIVES, 9);
61 integrator.addEventDetector(detector1);
62 TimeDetector detector2 = new TimeDetector(11, 1e-9, 100, 9 + 1e-15, 9 + 4.9);
63 integrator.addEventDetector(detector2);
64
65
66 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 20);
67
68
69 List<Event> events1 = detector1.getEvents();
70 Assert.assertEquals(1, events1.size());
71 Assert.assertEquals(9, events1.get(0).getT(), 0.0);
72 List<Event> events2 = detector2.getEvents();
73 Assert.assertEquals(0, events2.size());
74 }
75
76 @Test
77 public void testCloseEvents() {
78
79 double e = 1e-15;
80 ODEIntegrator integrator =
81 new DormandPrince853Integrator(e, 100.0, 1e-7, 1e-7);
82
83 TimeDetector detector1 = new TimeDetector(10, 1, 100, 5);
84 integrator.addEventDetector(detector1);
85 TimeDetector detector2 = new TimeDetector(10, 1, 100, 5.5);
86 integrator.addEventDetector(detector2);
87
88
89 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 20);
90
91
92 List<Event> events1 = detector1.getEvents();
93 Assert.assertEquals(1, events1.size());
94 Assert.assertEquals(5, events1.get(0).getT(), 0.0);
95 List<Event> events2 = detector2.getEvents();
96 Assert.assertEquals(1, events2.size());
97 Assert.assertEquals(5.5, events2.get(0).getT(), 0.0);
98 }
99
100 @Test
101 public void testSimultaneousEvents() {
102
103 ODEIntegrator integrator =
104 new DormandPrince853Integrator(10, 100.0, 1e-7, 1e-7);
105
106 TimeDetector detector1 = new TimeDetector(10, 1, 100, 5);
107 integrator.addEventDetector(detector1);
108 TimeDetector detector2 = new TimeDetector(10, 1, 100, 5);
109 integrator.addEventDetector(detector2);
110
111
112 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 20);
113
114
115 List<Event> events1 = detector1.getEvents();
116 Assert.assertEquals(1, events1.size());
117 Assert.assertEquals(5, events1.get(0).getT(), 0.0);
118 List<Event> events2 = detector2.getEvents();
119 Assert.assertEquals(1, events2.size());
120 Assert.assertEquals(5, events2.get(0).getT(), 0.0);
121 }
122
123
124
125
126
127
128 @Test
129 public void testSimultaneousEventsReset() {
130
131 double tol = 1e-10;
132 ODEIntegrator integrator =
133 new DormandPrince853Integrator(10, 100.0, 1e-7, 1e-7);
134 boolean[] firstEventOccurred = {false};
135 List<Event> events = new ArrayList<>();
136
137 TimeDetector detector1 = new TimeDetector(10, tol, 100, events, 5) {
138 @Override
139 public ODEEventHandler getHandler() {
140 return (state, detector, increasing) -> {
141 firstEventOccurred[0] = true;
142 super.getHandler().eventOccurred(state, detector, increasing);
143 return Action.RESET_STATE;
144 };
145 }
146 };
147 integrator.addEventDetector(detector1);
148
149 TimeDetector detector2 = new TimeDetector(1, tol, 100, events, 1, 3, 5) {
150 @Override
151 public double g(final ODEStateAndDerivative state) {
152 if (firstEventOccurred[0]) {
153 return super.g(state);
154 }
155 return new TimeDetector(1, tol, 100, 5).g(state);
156 }
157 };
158 integrator.addEventDetector(detector2);
159
160
161 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 20);
162
163
164
165 Assert.assertEquals(5, events.get(0).getT(), 0.0);
166 Assert.assertTrue(events.get(0).isIncreasing());
167 Assert.assertEquals(detector1, events.get(0).getDetector());
168 Assert.assertEquals(5, events.get(1).getT(), 0.0);
169 Assert.assertTrue(events.get(1).isIncreasing());
170 Assert.assertEquals(detector2, events.get(1).getDetector());
171 Assert.assertEquals(2, events.size());
172 }
173
174
175
176
177
178
179
180 @Test
181 public void testSimultaneousDiscontinuousEventsAfterReset() {
182
183 double t = FastMath.PI;
184 double tol = 1e-10;
185 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
186 List<Event> events = new ArrayList<>();
187
188 TimeDetector resetDetector =
189 new ResetDetector(10, tol, 100, events, new ODEState(t, new double[]{1e100, 0}), t);
190 integrator.addEventDetector(resetDetector);
191 List<BaseDetector> detectors = new ArrayList<>();
192 for (int i = 0; i < 2; i++) {
193 BaseDetector detector1 = new StateDetector(10, tol, 100, events, 0.0);
194 integrator.addEventDetector(detector1);
195 detectors.add(detector1);
196 }
197
198
199 integrator.integrate(new Equation(), new ODEState(0, new double[]{-1e100, 0}), 10);
200
201
202 Assert.assertEquals(t, events.get(0).getT(), tol);
203 Assert.assertTrue(events.get(0).isIncreasing());
204 Assert.assertEquals(resetDetector, events.get(0).getDetector());
205
206 Assert.assertEquals(t, events.get(1).getT(), tol);
207 Assert.assertTrue(events.get(1).isIncreasing());
208 Assert.assertEquals(detectors.get(0), events.get(1).getDetector());
209 Assert.assertEquals(t, events.get(2).getT(), tol);
210 Assert.assertTrue(events.get(2).isIncreasing());
211 Assert.assertEquals(detectors.get(1), events.get(2).getDetector());
212 Assert.assertEquals(events.size(), 3);
213 }
214
215
216
217
218
219
220 @Test
221 public void testFastSwitching() {
222
223
224 ODEIntegrator integrator =
225 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
226
227 TimeDetector detector1 = new TimeDetector(10, 0.2, 100, 9.9, 10.1, 12);
228 integrator.addEventDetector(detector1);
229
230
231 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 20);
232
233
234
235 List<Event> events1 = detector1.getEvents();
236 Assert.assertEquals(1, events1.size());
237 Assert.assertEquals(9.9, events1.get(0).getT(), 0.1);
238 Assert.assertEquals(true, events1.get(0).isIncreasing());
239 }
240
241
242 @Test
243 public void testTrickyCaseLower() {
244
245 double maxCheck = 10;
246 double tolerance = 1e-6;
247 double t1 = 1.0, t2 = 15, t3 = 16, t4 = 17, t5 = 18;
248
249 List<Event> events = new ArrayList<>();
250 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t3);
251 TimeDetector detectorB = new TimeDetector(maxCheck, tolerance, 100, events, -10, t1, t2, t5);
252 TimeDetector detectorC = new TimeDetector(maxCheck, tolerance, 100, Action.RESET_DERIVATIVES, events, t4);
253
254 ODEIntegrator integrator =
255 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
256 integrator.addEventDetector(detectorA);
257 integrator.addEventDetector(detectorB);
258 integrator.addEventDetector(detectorC);
259
260
261 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
262
263
264
265
266 Assert.assertEquals(5, events.size());
267 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
268 Assert.assertEquals(false, events.get(0).isIncreasing());
269 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
270 Assert.assertEquals(true, events.get(1).isIncreasing());
271 Assert.assertEquals(t3, events.get(2).getT(), tolerance);
272 Assert.assertEquals(true, events.get(2).isIncreasing());
273 Assert.assertEquals(t4, events.get(3).getT(), tolerance);
274 Assert.assertEquals(true, events.get(3).isIncreasing());
275 Assert.assertEquals(t5, events.get(4).getT(), tolerance);
276 Assert.assertEquals(false, events.get(4).isIncreasing());
277 }
278
279
280
281
282
283
284 @Test
285 public void testRootFindingTolerance() {
286
287 double maxCheck = 10;
288 double t2 = 11, t3 = t2 + 1e-5;
289 List<Event> events = new ArrayList<>();
290 TimeDetector detectorA = new TimeDetector(maxCheck, 1e-6, 100, events, t2);
291 TimeDetector detectorB = new FlatDetector(maxCheck, 0.5, 100, events, t3);
292 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
293 integrator.addEventDetector(detectorA);
294 integrator.addEventDetector(detectorB);
295
296
297 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
298
299
300
301
302 Assert.assertSame(detectorB, events.get(0).getDetector());
303 Assert.assertSame(detectorA, events.get(1).getDetector());
304 Assert.assertTrue(events.get(0).getT() < events.get(1).getT());
305
306
307 Assert.assertEquals(2, events.size());
308 Assert.assertEquals(t3, events.get(0).getT(), 0.5);
309 Assert.assertEquals(true, events.get(0).isIncreasing());
310 Assert.assertEquals(t2, events.get(1).getT(), 1e-6);
311 Assert.assertEquals(true, events.get(1).isIncreasing());
312 }
313
314
315 @Test
316 public void testRootPlusToleranceHasWrongSign() {
317
318 double maxCheck = 10;
319 double tolerance = 1e-6;
320 final double toleranceB = 0.3;
321 double t1 = 11, t2 = 11.1, t3 = 11.2;
322
323 List<Event> events = new ArrayList<>();
324 TimeDetector detectorA = new TimeDetector(maxCheck, 1e-6, 100, events, t2);
325 TimeDetector detectorB = new TimeDetector(maxCheck, toleranceB, 100, events, t1, t3);
326 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
327 integrator.addEventDetector(detectorA);
328 integrator.addEventDetector(detectorB);
329
330
331 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
332
333
334
335 Assert.assertEquals(3, events.size());
336 Assert.assertEquals(t1, events.get(0).getT(), toleranceB);
337 Assert.assertEquals(true, events.get(0).isIncreasing());
338 Assert.assertSame(detectorB, events.get(0).getDetector());
339 Assert.assertEquals(t3, events.get(1).getT(), toleranceB);
340 Assert.assertEquals(false, events.get(1).isIncreasing());
341 Assert.assertSame(detectorB, events.get(1).getDetector());
342 Assert.assertEquals(t2, events.get(2).getT(), tolerance);
343 Assert.assertEquals(true, events.get(2).isIncreasing());
344 Assert.assertSame(detectorA, events.get(2).getDetector());
345
346 for (int i = 1; i < events.size(); i++) {
347 Assert.assertTrue(events.get(i).getT() >= events.get(i - 1).getT());
348 }
349 }
350
351
352 @Test
353 public void testRootPlusToleranceHasWrongSignAndLessThanTb() {
354
355
356 double maxCheck = 10;
357 double tolerance = 0.5;
358 double t1 = 11, t2 = 11.4, t3 = 12.0;
359
360 List<Event> events = new ArrayList<>();
361 TimeDetector detectorB = new FlatDetector(maxCheck, tolerance, 100, events, t1, t2, t3);
362 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
363 integrator.addEventDetector(detectorB);
364
365
366 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
367
368
369
370 Assert.assertEquals(1, events.size());
371 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
372 Assert.assertEquals(true, events.get(0).isIncreasing());
373 Assert.assertSame(detectorB, events.get(0).getDetector());
374 }
375
376
377
378
379
380 @Test
381 public void testDoubleRoot() {
382
383 double maxCheck = 10;
384 double tolerance = 1e-6;
385 double t1 = 11;
386
387 List<Event> events = new ArrayList<>();
388 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t1);
389 TimeDetector detectorB = new TimeDetector(maxCheck, tolerance, 100, events, t1, t1);
390 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
391 integrator.addEventDetector(detectorA);
392 integrator.addEventDetector(detectorB);
393
394
395 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
396
397
398 Assert.assertEquals(1, events.size());
399 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
400 Assert.assertEquals(true, events.get(0).isIncreasing());
401 Assert.assertSame(detectorA, events.get(0).getDetector());
402
403 Assert.assertTrue(detectorB.g(state(t1)) == 0.0);
404 Assert.assertTrue(detectorB.g(state(t1 - 1e-6)) < 0);
405 Assert.assertTrue(detectorB.g(state(t1 + 1e-6)) < 0);
406 }
407
408
409
410
411
412 @Test
413 public void testDoubleRootOppositeSign() {
414
415 double maxCheck = 10;
416 double tolerance = 1e-6;
417 double t1 = 11;
418
419 List<Event> events = new ArrayList<>();
420 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t1);
421 TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, -20, t1, t1);
422 detectorB.g(state(t1));
423 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
424 integrator.addEventDetector(detectorA);
425 integrator.addEventDetector(detectorB);
426
427
428 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
429
430
431 Assert.assertEquals(1, events.size());
432 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
433 Assert.assertEquals(true, events.get(0).isIncreasing());
434 Assert.assertSame(detectorA, events.get(0).getDetector());
435
436 Assert.assertEquals(0.0, detectorB.g(state(t1)), 0.0);
437 Assert.assertTrue(detectorB.g(state(t1 - 1e-6)) > 0);
438 Assert.assertTrue(detectorB.g(state(t1 + 1e-6)) > 0);
439 }
440
441
442 @Test
443 public void testZeroAtBeginningAndEndOfInterval() {
444
445 double maxCheck = 10;
446 double tolerance = 1e-6;
447 double t1 = 10, t2 = 20;
448
449 List<Event> events = new ArrayList<>();
450 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t1, t2);
451 ODEIntegrator integrator =
452 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
453 integrator.addEventDetector(detectorA);
454
455
456 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
457
458
459 Assert.assertEquals(2, events.size());
460 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
461 Assert.assertEquals(true, events.get(0).isIncreasing());
462 Assert.assertSame(detectorA, events.get(0).getDetector());
463 Assert.assertEquals(t2, events.get(1).getT(), 0.0);
464 Assert.assertEquals(false, events.get(1).isIncreasing());
465 Assert.assertSame(detectorA, events.get(1).getDetector());
466 }
467
468
469 @Test
470 public void testZeroAtBeginningAndEndOfIntervalOppositeSign() {
471
472 double maxCheck = 10;
473 double tolerance = 1e-6;
474 double t1 = 10, t2 = 20;
475
476 List<Event> events = new ArrayList<>();
477 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, -10, t1, t2);
478 ODEIntegrator integrator =
479 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
480 integrator.addEventDetector(detectorA);
481
482
483 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
484
485
486 Assert.assertEquals(2, events.size());
487 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
488 Assert.assertEquals(false, events.get(0).isIncreasing());
489 Assert.assertSame(detectorA, events.get(0).getDetector());
490 Assert.assertEquals(t2, events.get(1).getT(), 0.0);
491 Assert.assertEquals(true, events.get(1).isIncreasing());
492 Assert.assertSame(detectorA, events.get(1).getDetector());
493 }
494
495
496 @Test
497 public void testMultipleBackups() {
498
499 double maxCheck = 10;
500 double tolerance = 1e-6;
501 double t1 = 11.0, t2 = 12, t3 = 13, t4 = 14, t5 = 15, t6 = 16, t7 = 17;
502
503 List<Event> events = new ArrayList<>();
504 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t6);
505 TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t1, t3, t4, t7);
506 TimeDetector detectorC = new ContinuousDetector(maxCheck, tolerance, 100, events, t2, t5);
507
508 ODEIntegrator integrator =
509 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
510 integrator.addEventDetector(detectorA);
511 integrator.addEventDetector(detectorB);
512 integrator.addEventDetector(detectorC);
513
514
515 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
516
517
518
519 Assert.assertEquals(5, events.size());
520 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
521 Assert.assertEquals(true, events.get(0).isIncreasing());
522 Assert.assertEquals(detectorB, events.get(0).getDetector());
523 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
524 Assert.assertEquals(true, events.get(1).isIncreasing());
525 Assert.assertEquals(detectorC, events.get(1).getDetector());
526
527
528
529
530
531
532
533
534
535
536 Assert.assertEquals(t5, events.get(2).getT(), tolerance);
537 Assert.assertEquals(false, events.get(2).isIncreasing());
538 Assert.assertEquals(detectorC, events.get(2).getDetector());
539 Assert.assertEquals(t6, events.get(3).getT(), tolerance);
540 Assert.assertEquals(true, events.get(3).isIncreasing());
541 Assert.assertEquals(detectorA, events.get(3).getDetector());
542 Assert.assertEquals(t7, events.get(4).getT(), tolerance);
543 Assert.assertEquals(false, events.get(4).isIncreasing());
544 Assert.assertEquals(detectorB, events.get(4).getDetector());
545 }
546
547
548 @Test
549 public void testEventCausedByStateReset() {
550
551 double maxCheck = 10;
552 double tolerance = 1e-6;
553 double t1 = 15.0;
554 final ODEState newState = new ODEState(t1, new double[]{-20, 0});
555
556 List<Event> events = new ArrayList<>();
557 TimeDetector detectorA = new ResetDetector(maxCheck, tolerance, 100, events, newState, t1);
558 BaseDetector detectorB = new StateDetector(maxCheck, tolerance, 100, events, -1);
559
560 ODEIntegrator integrator =
561 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
562 integrator.addEventDetector(detectorA);
563 integrator.addEventDetector(detectorB);
564
565
566 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 40.0);
567
568
569
570 Assert.assertEquals(3, events.size());
571 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
572 Assert.assertEquals(true, events.get(0).isIncreasing());
573 Assert.assertEquals(detectorA, events.get(0).getDetector());
574 Assert.assertEquals(t1, events.get(1).getT(), tolerance);
575 Assert.assertEquals(false, events.get(1).isIncreasing());
576 Assert.assertEquals(detectorB, events.get(1).getDetector());
577 Assert.assertEquals(t1 + 19, events.get(2).getT(), tolerance);
578 Assert.assertEquals(true, events.get(2).isIncreasing());
579 Assert.assertEquals(detectorB, events.get(2).getDetector());
580 }
581
582
583 @Test
584 public void testConvergenceTooTight() {
585
586 double maxCheck = 10;
587 double tolerance = 1e-18;
588 double t1 = 15;
589
590 List<Event> events = new ArrayList<>();
591 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t1);
592 ODEIntegrator integrator =
593 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
594 integrator.addEventDetector(detectorA);
595
596
597 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
598
599
600 Assert.assertEquals(1, events.size());
601 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
602 Assert.assertEquals(true, events.get(0).isIncreasing());
603 Assert.assertSame(detectorA, events.get(0).getDetector());
604 }
605
606
607 @Test
608 public void testToleranceMismatch() {
609
610 double maxCheck = 10;
611 double tolerance = 1e-18;
612 double t1 = 15.1;
613
614 List<Event> events = new ArrayList<>();
615 TimeDetector detectorA = new FlatDetector(maxCheck, tolerance, 100, events, t1);
616 ODEIntegrator integrator =
617 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
618 integrator.addEventDetector(detectorA);
619
620
621 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
622
623
624 Assert.assertEquals(1, events.size());
625
626 Assert.assertEquals(t1, events.get(0).getT(), 1e-3);
627 Assert.assertEquals(true, events.get(0).isIncreasing());
628 Assert.assertSame(detectorA, events.get(0).getDetector());
629 }
630
631
632
633
634
635
636
637
638 @Test
639 public void testEventChangesGFunctionDefinition() {
640
641 double maxCheck = 5;
642 double tolerance = 1e-6;
643 double t1 = 11, t2 = 19;
644
645 List<Event> events = new ArrayList<>();
646
647 boolean[] swap = new boolean[1];
648 final TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t1) {
649 @Override
650 public ODEEventHandler getHandler() {
651 return (state, detector, increasing) -> {
652 swap[0] = true;
653 return super.getHandler().eventOccurred(state, detector, increasing);
654 };
655 }
656 };
657 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
658 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
659
660 @Override
661 public double g(ODEStateAndDerivative state) {
662 if (swap[0]) {
663 return detectorB.g(state);
664 } else {
665 return -1;
666 }
667 }
668
669 };
670 ODEIntegrator integrator =
671 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
672 integrator.addEventDetector(detectorA);
673 integrator.addEventDetector(detectorC);
674
675
676 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
677
678
679 Assert.assertEquals(2, events.size());
680 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
681 Assert.assertEquals(true, events.get(0).isIncreasing());
682 Assert.assertSame(detectorA, events.get(0).getDetector());
683 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
684 Assert.assertEquals(true, events.get(1).isIncreasing());
685 Assert.assertSame(detectorC, events.get(1).getDetector());
686 }
687
688
689
690
691
692
693 @Test
694 public void testEventChangesGFunctionDefinitionCancel() {
695
696 double maxCheck = 5;
697 double tolerance = 1e-6;
698 double t1 = 11, t2 = 11.1;
699
700 List<Event> events = new ArrayList<>();
701
702 boolean[] swap = new boolean[1];
703 final TimeDetector detectorA =
704 new ContinuousDetector(maxCheck, tolerance, 100, Action.RESET_EVENTS, events, t1) {
705 @Override
706 public ODEEventHandler getHandler() {
707 return (state, detector, increasing) -> {
708 swap[0] = true;
709 return super.getHandler().eventOccurred(state, detector, increasing);
710 };
711 }
712 };
713 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
714 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
715
716 @Override
717 public double g(ODEStateAndDerivative state) {
718 if (!swap[0]) {
719 return detectorB.g(state);
720 } else {
721 return -1;
722 }
723 }
724
725 };
726 ODEIntegrator integrator =
727 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
728 integrator.addEventDetector(detectorA);
729 integrator.addEventDetector(detectorC);
730
731
732 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
733
734
735 Assert.assertEquals(1, events.size());
736 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
737 Assert.assertEquals(true, events.get(0).isIncreasing());
738 Assert.assertSame(detectorA, events.get(0).getDetector());
739 }
740
741
742
743
744
745
746 @Test
747 public void testEventChangesGFunctionDefinitionDelay() {
748
749 double maxCheck = 5;
750 double tolerance = 1e-6;
751 double t1 = 11, t2 = 11.1, t3 = 11.2;
752
753 List<Event> events = new ArrayList<>();
754
755 boolean[] swap = new boolean[1];
756 final TimeDetector detectorA =
757 new ContinuousDetector(maxCheck, tolerance, 100, Action.RESET_EVENTS, events, t1) {
758 @Override
759 public ODEEventHandler getHandler() {
760 return (state, detector, increasing) -> {
761 swap[0] = true;
762 return super.getHandler().eventOccurred(state, detector, increasing);
763 };
764 }
765 };
766 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
767 final TimeDetector detectorD = new ContinuousDetector(maxCheck, tolerance, 100, events, t3);
768 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
769
770 @Override
771 public double g(ODEStateAndDerivative state) {
772 if (!swap[0]) {
773 return detectorB.g(state);
774 } else {
775 return detectorD.g(state);
776 }
777 }
778
779 };
780 ODEIntegrator integrator =
781 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
782 integrator.addEventDetector(detectorA);
783 integrator.addEventDetector(detectorC);
784
785
786 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
787
788
789 Assert.assertEquals(2, events.size());
790 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
791 Assert.assertEquals(true, events.get(0).isIncreasing());
792 Assert.assertSame(detectorA, events.get(0).getDetector());
793 Assert.assertEquals(t3, events.get(1).getT(), tolerance);
794 Assert.assertEquals(true, events.get(1).isIncreasing());
795 Assert.assertSame(detectorC, events.get(1).getDetector());
796 }
797
798
799
800
801
802
803 @Test
804 public void testEventChangesGFunctionDefinitionAccelerate() {
805
806 double maxCheck = 5;
807 double tolerance = 1e-6;
808 double t1 = 11, t2 = 11.1, t3 = 11.2;
809
810 List<Event> events = new ArrayList<>();
811
812 boolean[] swap = new boolean[1];
813 final TimeDetector detectorA =
814 new ContinuousDetector(maxCheck, tolerance, 100, Action.RESET_EVENTS, events, t1) {
815 @Override
816 public ODEEventHandler getHandler() {
817 return (state, detector, increasing) -> {
818 swap[0] = true;
819 return super.getHandler().eventOccurred(state, detector, increasing);
820 };
821 }
822 };
823 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
824 final TimeDetector detectorD = new ContinuousDetector(maxCheck, tolerance, 100, events, t3);
825 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
826
827 @Override
828 public double g(ODEStateAndDerivative state) {
829 if (swap[0]) {
830 return detectorB.g(state);
831 } else {
832 return detectorD.g(state);
833 }
834 }
835
836 };
837 ODEIntegrator integrator =
838 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
839 integrator.addEventDetector(detectorA);
840 integrator.addEventDetector(detectorC);
841
842
843 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
844
845
846 Assert.assertEquals(2, events.size());
847 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
848 Assert.assertEquals(true, events.get(0).isIncreasing());
849 Assert.assertSame(detectorA, events.get(0).getDetector());
850 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
851 Assert.assertEquals(true, events.get(1).isIncreasing());
852 Assert.assertSame(detectorC, events.get(1).getDetector());
853 }
854
855
856 @Test
857 public void testToleranceStop() {
858
859 double maxCheck = 10;
860 double tolerance = 1e-18;
861 double t1 = 15.1;
862
863 List<Event> events = new ArrayList<>();
864 TimeDetector detectorA = new FlatDetector(maxCheck, tolerance, 100, Action.STOP, events, t1);
865 ODEIntegrator integrator =
866 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
867 integrator.addEventDetector(detectorA);
868
869
870 ODEStateAndDerivative finalState =
871 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
872
873
874 Assert.assertEquals(1, events.size());
875
876 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
877 Assert.assertEquals(true, events.get(0).isIncreasing());
878 Assert.assertSame(detectorA, events.get(0).getDetector());
879 Assert.assertEquals(t1, finalState.getTime(), tolerance);
880
881
882 finalState = integrator.integrate(new Equation(), finalState, 30.0);
883
884
885 Assert.assertEquals(30.0, finalState.getTime(), 0.0);
886 }
887
888
889
890
891
892 @Test
893 public void testLongInitialZero() {
894
895 double maxCheck = 10;
896 double tolerance = 1;
897
898 List<Event> events = new ArrayList<>();
899 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, Action.STOP, events, -50) {
900 @Override
901 public double g(ODEStateAndDerivative state) {
902 if (state.getTime() < 2) {
903 return 0;
904 } else {
905 return super.g(state);
906 }
907 }
908 };
909 ODEIntegrator integrator =
910 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
911 integrator.addEventDetector(detectorA);
912
913
914 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
915
916
917 Assert.assertEquals(0, events.size());
918 }
919
920
921
922
923
924
925 @Test
926 public void testShortBracketingInterval() {
927
928 double maxCheck = 10;
929 double tolerance = 1e-6;
930 final double t1 = FastMath.nextUp(10.0), t2 = 10.5;
931
932 List<Event> events = new ArrayList<>();
933
934 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events) {
935 @Override
936 public double g(ODEStateAndDerivative state) {
937 final double t = state.getTime();
938 if (t < t1) {
939 return -1;
940 } else if (t < t2) {
941 return 1;
942 } else {
943 return -1;
944 }
945 }
946 };
947 TimeDetector detectorB = new TimeDetector(maxCheck, tolerance, 100, events, t1);
948 ODEIntegrator integrator =
949 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
950 integrator.addEventDetector(detectorA);
951 integrator.addEventDetector(detectorB);
952
953
954 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 30.0);
955
956
957 Assert.assertEquals(3, events.size());
958 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
959 Assert.assertEquals(true, events.get(0).isIncreasing());
960 Assert.assertSame(detectorA, events.get(0).getDetector());
961 Assert.assertEquals(t1, events.get(1).getT(), tolerance);
962 Assert.assertEquals(true, events.get(1).isIncreasing());
963 Assert.assertSame(detectorB, events.get(1).getDetector());
964 Assert.assertEquals(t2, events.get(2).getT(), tolerance);
965 Assert.assertEquals(false, events.get(2).isIncreasing());
966 Assert.assertSame(detectorA, events.get(2).getDetector());
967 }
968
969
970 @Test
971 public void testEventStepHandler() {
972
973 double tolerance = 1e-18;
974 ODEIntegrator integrator =
975 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
976 integrator.addEventDetector(new TimeDetector(100, tolerance, 100, 5));
977 StepHandler stepHandler = new StepHandler();
978 integrator.addStepHandler(stepHandler);
979
980
981 ODEStateAndDerivative finalState = integrator
982 .integrate(new Equation(), new ODEState(0, new double[2]), 10);
983
984
985 Assert.assertEquals(10.0, finalState.getTime(), tolerance);
986 Assert.assertEquals(0.0,
987 stepHandler.initialState.getTime(), tolerance);
988 Assert.assertEquals(10.0, stepHandler.finalTime, tolerance);
989 Assert.assertEquals(10.0,
990 stepHandler.finalState.getTime(), tolerance);
991 ODEStateInterpolator interpolator = stepHandler.interpolators.get(0);
992 Assert.assertEquals(0.0,
993 interpolator.getPreviousState().getTime(), tolerance);
994 Assert.assertEquals(5.0,
995 interpolator.getCurrentState().getTime(), tolerance);
996 interpolator = stepHandler.interpolators.get(1);
997 Assert.assertEquals(5.0,
998 interpolator.getPreviousState().getTime(), tolerance);
999 Assert.assertEquals(10.0,
1000 interpolator.getCurrentState().getTime(), tolerance);
1001 Assert.assertEquals(2, stepHandler.interpolators.size());
1002 }
1003
1004
1005 @Test
1006 public void testEventCausedByDerivativesReset() {
1007
1008 TimeDetector detectorA = new TimeDetector(10, 1e-6, 100, Action.RESET_STATE, 15.0) {
1009 @Override
1010 public ODEEventHandler getHandler() {
1011 return new ODEEventHandler() {
1012 @Override
1013 public Action eventOccurred(ODEStateAndDerivative state, ODEEventDetector detector, boolean increasing) {
1014 return Action.RESET_STATE;
1015 }
1016 @Override
1017 public ODEState resetState(ODEEventDetector detector, ODEStateAndDerivative state) {
1018 return null;
1019 }
1020 };
1021 }
1022 };
1023 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1024 integrator.addEventDetector(detectorA);
1025
1026 try {
1027
1028 integrator.integrate(new Equation(), new ODEState(0, new double[2]), 20.0);
1029 Assert.fail("Expected Exception");
1030 } catch (NullPointerException e) {
1031
1032 }
1033 }
1034
1035 @Test
1036 public void testResetChangesSign() {
1037 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
1038 public int getDimension() { return 1; }
1039 public double[] computeDerivatives(double t, double[] y) { return new double[] { 1.0 }; }
1040 };
1041
1042 LutherIntegrator integrator = new LutherIntegrator(20.0);
1043 final double small = 1.0e-10;
1044 ResetChangesSignGenerator eventsGenerator = new ResetChangesSignGenerator(6.0, 9.0, -0.5 * small, 8.0, small, 1000);
1045 integrator.addEventDetector(eventsGenerator);
1046 final ODEStateAndDerivative end = integrator.integrate(equation, new ODEState(0.0, new double[1]), 100.0);
1047 Assert.assertEquals(2, eventsGenerator.getCount());
1048 Assert.assertEquals(9.0, end.getCompleteState()[0], 1.0e-12);
1049 Assert.assertEquals(9.0 + 0.5 * small, end.getTime(), 1.0e-12);
1050 }
1051
1052
1053
1054
1055
1056
1057 @Test
1058 public void testCloseEventsFirstOneIsResetReverse() {
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069 ODEIntegrator integrator =
1070 new DormandPrince853Integrator(10, 100.0, 1e-7, 1e-7);
1071
1072
1073 double t1 = -1;
1074 TimeDetector detector1 = new TimeDetector(10, 1e-9, 100, Action.RESET_DERIVATIVES, t1);
1075 integrator.addEventDetector(detector1);
1076 TimeDetector detector2 = new TimeDetector(11, 1e-9, 100, t1 - 1e-15, t1 - 4.9);
1077 integrator.addEventDetector(detector2);
1078
1079
1080 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -20);
1081
1082
1083 List<Event> events1 = detector1.getEvents();
1084 Assert.assertEquals(1, events1.size());
1085 Assert.assertEquals(t1, events1.get(0).getT(), 0.0);
1086 List<Event> events2 = detector2.getEvents();
1087 Assert.assertEquals(0, events2.size());
1088 }
1089
1090 @Test
1091 public void testCloseEventsReverse() {
1092
1093 double e = 1e-15;
1094 ODEIntegrator integrator =
1095 new DormandPrince853Integrator(e, 100.0, 1e-7, 1e-7);
1096
1097 TimeDetector detector1 = new TimeDetector(10, 1, 100, -5);
1098 integrator.addEventDetector(detector1);
1099 TimeDetector detector2 = new TimeDetector(10, 1, 100, -5.5);
1100 integrator.addEventDetector(detector2);
1101
1102
1103 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -20);
1104
1105
1106 List<Event> events1 = detector1.getEvents();
1107 Assert.assertEquals(1, events1.size());
1108 Assert.assertEquals(-5, events1.get(0).getT(), 0.0);
1109 List<Event> events2 = detector2.getEvents();
1110 Assert.assertEquals(1, events2.size());
1111 Assert.assertEquals(-5.5, events2.get(0).getT(), 0.0);
1112 }
1113
1114 @Test
1115 public void testSimultaneousEventsReverse() {
1116
1117 ODEIntegrator integrator =
1118 new DormandPrince853Integrator(10, 100.0, 1e-7, 1e-7);
1119
1120 TimeDetector detector1 = new TimeDetector(10, 1, 100, -5);
1121 integrator.addEventDetector(detector1);
1122 TimeDetector detector2 = new TimeDetector(10, 1, 100, -5);
1123 integrator.addEventDetector(detector2);
1124
1125
1126 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -20);
1127
1128
1129 List<Event> events1 = detector1.getEvents();
1130 Assert.assertEquals(1, events1.size());
1131 Assert.assertEquals(-5, events1.get(0).getT(), 0.0);
1132 List<Event> events2 = detector2.getEvents();
1133 Assert.assertEquals(1, events2.size());
1134 Assert.assertEquals(-5, events2.get(0).getT(), 0.0);
1135 }
1136
1137
1138
1139
1140
1141
1142 @Test
1143 public void testSimultaneousEventsResetReverse() {
1144
1145 double tol = 1e-10;
1146 ODEIntegrator integrator =
1147 new DormandPrince853Integrator(10, 100.0, 1e-7, 1e-7);
1148 boolean[] firstEventOccurred = {false};
1149 List<Event> events = new ArrayList<>();
1150
1151 TimeDetector detector1 = new TimeDetector(10, tol, 100, events, -5) {
1152 @Override
1153 public ODEEventHandler getHandler() {
1154 return (state, detector, increasing) -> {
1155 firstEventOccurred[0] = true;
1156 super.getHandler().eventOccurred(state, detector, increasing);
1157 return Action.RESET_STATE;
1158 };
1159 }
1160 };
1161 integrator.addEventDetector(detector1);
1162
1163 TimeDetector detector2 = new TimeDetector(1, tol, 100, events, -1, -3, -5) {
1164 @Override
1165 public double g(final ODEStateAndDerivative state) {
1166 if (firstEventOccurred[0]) {
1167 return super.g(state);
1168 }
1169 return new TimeDetector(1, tol, 100, -5).g(state);
1170 }
1171 };
1172 integrator.addEventDetector(detector2);
1173
1174
1175 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -20);
1176
1177
1178
1179 Assert.assertEquals(-5, events.get(0).getT(), 0.0);
1180 Assert.assertTrue(events.get(0).isIncreasing());
1181 Assert.assertEquals(detector1, events.get(0).getDetector());
1182 Assert.assertEquals(-5, events.get(1).getT(), 0.0);
1183 Assert.assertTrue(events.get(1).isIncreasing());
1184 Assert.assertEquals(detector2, events.get(1).getDetector());
1185 Assert.assertEquals(2, events.size());
1186 }
1187
1188
1189
1190
1191
1192
1193
1194 @Test
1195 public void testSimultaneousDiscontinuousEventsAfterResetReverse() {
1196
1197 double t = -FastMath.PI;
1198 double tol = 1e-10;
1199 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1200 List<Event> events = new ArrayList<>();
1201
1202 TimeDetector resetDetector =
1203 new ResetDetector(10, tol, 100, events, new ODEState(t, new double[]{1e100, 0}), t);
1204 integrator.addEventDetector(resetDetector);
1205 List<BaseDetector> detectors = new ArrayList<>();
1206 for (int i = 0; i < 2; i++) {
1207 BaseDetector detector1 = new StateDetector(10, tol, 100, events, 0.0);
1208 integrator.addEventDetector(detector1);
1209 detectors.add(detector1);
1210 }
1211
1212
1213 integrator.integrate(new Equation(), new ODEState(0, new double[]{-1e100, 0}), -10);
1214
1215
1216 Assert.assertEquals(t, events.get(0).getT(), tol);
1217 Assert.assertTrue(events.get(0).isIncreasing());
1218 Assert.assertEquals(resetDetector, events.get(0).getDetector());
1219
1220 Assert.assertEquals(t, events.get(1).getT(), tol);
1221 Assert.assertFalse(events.get(1).isIncreasing());
1222 Assert.assertEquals(detectors.get(0), events.get(1).getDetector());
1223 Assert.assertEquals(t, events.get(2).getT(), tol);
1224 Assert.assertFalse(events.get(2).isIncreasing());
1225 Assert.assertEquals(detectors.get(1), events.get(2).getDetector());
1226 Assert.assertEquals(events.size(), 3);
1227 }
1228
1229
1230
1231
1232
1233
1234 @Test
1235 public void testFastSwitchingReverse() {
1236
1237
1238 ODEIntegrator integrator =
1239 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1240
1241 TimeDetector detector1 = new TimeDetector(10, 0.2, 100, -9.9, -10.1, -12);
1242 integrator.addEventDetector(detector1);
1243
1244
1245 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -20);
1246
1247
1248
1249 List<Event> events1 = detector1.getEvents();
1250 Assert.assertEquals(1, events1.size());
1251 Assert.assertEquals(-9.9, events1.get(0).getT(), 0.2);
1252 Assert.assertEquals(true, events1.get(0).isIncreasing());
1253 }
1254
1255
1256 @Test
1257 public void testTrickyCaseLowerReverse() {
1258
1259 double maxCheck = 10;
1260 double tolerance = 1e-6;
1261 double t1 = -1.0, t2 = -15, t3 = -16, t4 = -17, t5 = -18;
1262
1263 List<Event> events = new ArrayList<>();
1264 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t3);
1265 TimeDetector detectorB = new TimeDetector(maxCheck, tolerance, 100, events, -50, t1, t2, t5);
1266 TimeDetector detectorC = new TimeDetector(maxCheck, tolerance, 100, Action.RESET_DERIVATIVES, events, t4);
1267
1268 ODEIntegrator integrator =
1269 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1270 integrator.addEventDetector(detectorA);
1271 integrator.addEventDetector(detectorB);
1272 integrator.addEventDetector(detectorC);
1273
1274
1275 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1276
1277
1278
1279
1280 Assert.assertEquals(5, events.size());
1281 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1282 Assert.assertEquals(false, events.get(0).isIncreasing());
1283 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
1284 Assert.assertEquals(true, events.get(1).isIncreasing());
1285 Assert.assertEquals(t3, events.get(2).getT(), tolerance);
1286 Assert.assertEquals(true, events.get(2).isIncreasing());
1287 Assert.assertEquals(t4, events.get(3).getT(), tolerance);
1288 Assert.assertEquals(true, events.get(3).isIncreasing());
1289 Assert.assertEquals(t5, events.get(4).getT(), tolerance);
1290 Assert.assertEquals(false, events.get(4).isIncreasing());
1291 }
1292
1293
1294
1295
1296
1297
1298 @Test
1299 public void testRootFindingToleranceReverse() {
1300
1301 double maxCheck = 10;
1302 double t2 = -11, t3 = t2 - 1e-5;
1303 List<Event> events = new ArrayList<>();
1304 TimeDetector detectorA = new TimeDetector(maxCheck, 1e-6, 100, events, t2);
1305 FlatDetector detectorB = new FlatDetector(maxCheck, 0.5, 100, events, t3);
1306 ODEIntegrator integrator =
1307 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1308 integrator.addEventDetector(detectorA);
1309 integrator.addEventDetector(detectorB);
1310
1311
1312 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1313
1314
1315
1316
1317 Assert.assertSame(detectorB, events.get(0).getDetector());
1318 Assert.assertSame(detectorA, events.get(1).getDetector());
1319 Assert.assertTrue(events.get(0).getT() > events.get(1).getT());
1320
1321
1322 Assert.assertEquals(2, events.size());
1323 Assert.assertEquals(t3, events.get(0).getT(), 0.5);
1324 Assert.assertEquals(true, events.get(0).isIncreasing());
1325 Assert.assertEquals(t2, events.get(1).getT(), 1e-6);
1326 Assert.assertEquals(true, events.get(1).isIncreasing());
1327 }
1328
1329
1330 @Test
1331 public void testRootPlusToleranceHasWrongSignReverse() {
1332
1333 double maxCheck = 10;
1334 double tolerance = 1e-6;
1335 final double toleranceB = 0.3;
1336 double t1 = -11, t2 = -11.1, t3 = -11.2;
1337
1338 List<Event> events = new ArrayList<>();
1339 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t2);
1340 TimeDetector detectorB = new TimeDetector(maxCheck, toleranceB, 100, events, -50, t1, t3);
1341 ODEIntegrator integrator =
1342 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1343 integrator.addEventDetector(detectorA);
1344 integrator.addEventDetector(detectorB);
1345
1346
1347 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1348
1349
1350
1351 Assert.assertEquals(3, events.size());
1352 Assert.assertEquals(t1, events.get(0).getT(), toleranceB);
1353 Assert.assertEquals(true, events.get(0).isIncreasing());
1354 Assert.assertSame(detectorB, events.get(0).getDetector());
1355 Assert.assertEquals(t3, events.get(1).getT(), toleranceB);
1356 Assert.assertEquals(false, events.get(1).isIncreasing());
1357 Assert.assertSame(detectorB, events.get(1).getDetector());
1358 Assert.assertEquals(t2, events.get(2).getT(), tolerance);
1359 Assert.assertEquals(true, events.get(2).isIncreasing());
1360 Assert.assertSame(detectorA, events.get(2).getDetector());
1361
1362 Assert.assertTrue(events.get(0).getT() >= events.get(1).getT());
1363 Assert.assertTrue(events.get(1).getT() >= events.get(2).getT());
1364 }
1365
1366
1367 @Test
1368 public void testRootPlusToleranceHasWrongSignAndLessThanTbReverse() {
1369
1370
1371 double maxCheck = 10;
1372 double tolerance = 0.5;
1373 double t1 = -11, t2 = -11.4, t3 = -12.0;
1374
1375 List<Event> events = new ArrayList<>();
1376 TimeDetector detectorB = new FlatDetector(maxCheck, tolerance, 100, events, t1, t2, t3);
1377 ODEIntegrator integrator =
1378 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1379 integrator.addEventDetector(detectorB);
1380
1381
1382 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1383
1384
1385
1386 Assert.assertEquals(1, events.size());
1387 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1388 Assert.assertEquals(true, events.get(0).isIncreasing());
1389 Assert.assertSame(detectorB, events.get(0).getDetector());
1390 }
1391
1392
1393
1394
1395
1396 @Test
1397 public void testDoubleRootReverse() {
1398
1399 double maxCheck = 10;
1400 double tolerance = 1e-6;
1401 double t1 = -11;
1402
1403 List<Event> events = new ArrayList<>();
1404 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t1);
1405 TimeDetector detectorB = new TimeDetector(maxCheck, tolerance, 100, events, t1, t1);
1406 ODEIntegrator integrator =
1407 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1408 integrator.addEventDetector(detectorA);
1409 integrator.addEventDetector(detectorB);
1410
1411
1412 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1413
1414
1415 Assert.assertEquals(1, events.size());
1416 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
1417 Assert.assertEquals(true, events.get(0).isIncreasing());
1418 Assert.assertSame(detectorA, events.get(0).getDetector());
1419
1420 Assert.assertTrue(detectorB.g(state(t1)) == 0.0);
1421 Assert.assertTrue(detectorB.g(state(t1 + 1e-6)) < 0);
1422 Assert.assertTrue(detectorB.g(state(t1 - 1e-6)) < 0);
1423 }
1424
1425
1426
1427
1428
1429 @Test
1430 public void testDoubleRootOppositeSignReverse() {
1431
1432 double maxCheck = 10;
1433 double tolerance = 1e-6;
1434 double t1 = -11;
1435
1436 List<Event> events = new ArrayList<>();
1437 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events, t1);
1438 TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, -50, t1, t1);
1439 detectorB.g(state(t1));
1440 ODEIntegrator integrator =
1441 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1442 integrator.addEventDetector(detectorA);
1443 integrator.addEventDetector(detectorB);
1444
1445
1446 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1447
1448
1449 Assert.assertEquals(1, events.size());
1450 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
1451 Assert.assertEquals(true, events.get(0).isIncreasing());
1452 Assert.assertSame(detectorA, events.get(0).getDetector());
1453
1454 Assert.assertEquals(0.0, detectorB.g(state(t1)), 0.0);
1455 Assert.assertTrue(detectorB.g(state(t1 + 1e-6)) > 0);
1456 Assert.assertTrue(detectorB.g(state(t1 - 1e-6)) > 0);
1457 }
1458
1459
1460 @Test
1461 public void testZeroAtBeginningAndEndOfIntervalReverse() {
1462
1463 double maxCheck = 10;
1464 double tolerance = 1e-6;
1465 double t1 = -10, t2 = -20;
1466
1467 List<Event> events = new ArrayList<>();
1468 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, -50, t1, t2);
1469 ODEIntegrator integrator =
1470 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1471 integrator.addEventDetector(detectorA);
1472
1473
1474 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1475
1476
1477 Assert.assertEquals(2, events.size());
1478 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
1479 Assert.assertEquals(true, events.get(0).isIncreasing());
1480 Assert.assertSame(detectorA, events.get(0).getDetector());
1481 Assert.assertEquals(t2, events.get(1).getT(), 0.0);
1482 Assert.assertEquals(false, events.get(1).isIncreasing());
1483 Assert.assertSame(detectorA, events.get(1).getDetector());
1484 }
1485
1486
1487 @Test
1488 public void testZeroAtBeginningAndEndOfIntervalOppositeSignReverse() {
1489
1490 double maxCheck = 10;
1491 double tolerance = 1e-6;
1492 double t1 = -10, t2 = -20;
1493
1494 List<Event> events = new ArrayList<>();
1495 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t1, t2);
1496 ODEIntegrator integrator =
1497 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1498 integrator.addEventDetector(detectorA);
1499
1500
1501 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1502
1503
1504 Assert.assertEquals(2, events.size());
1505 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
1506 Assert.assertEquals(false, events.get(0).isIncreasing());
1507 Assert.assertSame(detectorA, events.get(0).getDetector());
1508 Assert.assertEquals(t2, events.get(1).getT(), 0.0);
1509 Assert.assertEquals(true, events.get(1).isIncreasing());
1510 Assert.assertSame(detectorA, events.get(1).getDetector());
1511 }
1512
1513
1514 @Test
1515 public void testMultipleBackupsReverse() {
1516
1517 double maxCheck = 10;
1518 double tolerance = 1e-6;
1519 double t1 = -11.0, t2 = -12, t3 = -13, t4 = -14, t5 = -15, t6 = -16, t7 = -17;
1520
1521 List<Event> events = new ArrayList<>();
1522 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t6);
1523 TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, -50, t1, t3, t4, t7);
1524 TimeDetector detectorC = new ContinuousDetector(maxCheck, tolerance, 100, events, -50, t2, t5);
1525
1526 ODEIntegrator integrator =
1527 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1528 integrator.addEventDetector(detectorA);
1529 integrator.addEventDetector(detectorB);
1530 integrator.addEventDetector(detectorC);
1531
1532
1533 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1534
1535
1536
1537 Assert.assertEquals(5, events.size());
1538 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1539 Assert.assertEquals(true, events.get(0).isIncreasing());
1540 Assert.assertEquals(detectorB, events.get(0).getDetector());
1541 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
1542 Assert.assertEquals(true, events.get(1).isIncreasing());
1543 Assert.assertEquals(detectorC, events.get(1).getDetector());
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554 Assert.assertEquals(t5, events.get(2).getT(), tolerance);
1555 Assert.assertEquals(false, events.get(2).isIncreasing());
1556 Assert.assertEquals(detectorC, events.get(2).getDetector());
1557 Assert.assertEquals(t6, events.get(3).getT(), tolerance);
1558 Assert.assertEquals(true, events.get(3).isIncreasing());
1559 Assert.assertEquals(detectorA, events.get(3).getDetector());
1560 Assert.assertEquals(t7, events.get(4).getT(), tolerance);
1561 Assert.assertEquals(false, events.get(4).isIncreasing());
1562 Assert.assertEquals(detectorB, events.get(4).getDetector());
1563 }
1564
1565
1566 @Test
1567 public void testEventCausedByStateResetReverse() {
1568
1569 double maxCheck = 10;
1570 double tolerance = 1e-6;
1571 double t1 = -15.0;
1572 final ODEState newState = new ODEState(t1, new double[]{20, 0});
1573
1574 List<Event> events = new ArrayList<>();
1575 TimeDetector detectorA = new ResetDetector(maxCheck, tolerance, 100, events, newState, t1);
1576 BaseDetector detectorB = new StateDetector(maxCheck, tolerance, 100, events, 1);
1577
1578 ODEIntegrator integrator =
1579 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1580 integrator.addEventDetector(detectorA);
1581 integrator.addEventDetector(detectorB);
1582
1583
1584 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -40.0);
1585
1586
1587
1588 Assert.assertEquals(3, events.size());
1589 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1590 Assert.assertEquals(true, events.get(0).isIncreasing());
1591 Assert.assertEquals(detectorA, events.get(0).getDetector());
1592 Assert.assertEquals(t1, events.get(1).getT(), tolerance);
1593 Assert.assertEquals(false, events.get(1).isIncreasing());
1594 Assert.assertEquals(detectorB, events.get(1).getDetector());
1595 Assert.assertEquals(t1 - 19, events.get(2).getT(), tolerance);
1596 Assert.assertEquals(true, events.get(2).isIncreasing());
1597 Assert.assertEquals(detectorB, events.get(2).getDetector());
1598 }
1599
1600
1601 @Test
1602 public void testConvergenceTooTightReverse() {
1603
1604 double maxCheck = 10;
1605 double tolerance = 1e-18;
1606 double t1 = -15;
1607
1608 List<Event> events = new ArrayList<>();
1609 TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t1);
1610 ODEIntegrator integrator =
1611 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1612 integrator.addEventDetector(detectorA);
1613
1614
1615 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1616
1617
1618 Assert.assertEquals(1, events.size());
1619 Assert.assertEquals(t1, events.get(0).getT(), 0.0);
1620 Assert.assertEquals(true, events.get(0).isIncreasing());
1621 Assert.assertSame(detectorA, events.get(0).getDetector());
1622 }
1623
1624
1625 @Test
1626 public void testToleranceMismatchReverse() {
1627
1628 double maxCheck = 10;
1629 double tolerance = 1e-18;
1630 double t1 = -15.1;
1631
1632 List<Event> events = new ArrayList<>();
1633 TimeDetector detectorA = new FlatDetector(maxCheck, tolerance, 100, events, t1);
1634 ODEIntegrator integrator =
1635 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1636 integrator.addEventDetector(detectorA);
1637
1638
1639 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1640
1641
1642 Assert.assertEquals(1, events.size());
1643
1644 Assert.assertEquals(t1, events.get(0).getT(), 1e-3);
1645 Assert.assertEquals(true, events.get(0).isIncreasing());
1646 Assert.assertSame(detectorA, events.get(0).getDetector());
1647 }
1648
1649
1650
1651
1652
1653
1654
1655 @Test
1656 public void testEventChangesGFunctionDefinitionReverse() {
1657
1658 double maxCheck = 5;
1659 double tolerance = 1e-6;
1660 double t1 = -11, t2 = -19;
1661
1662 List<Event> events = new ArrayList<>();
1663
1664 boolean[] swap = new boolean[1];
1665 final TimeDetector detectorA = new ContinuousDetector(maxCheck, tolerance, 100, events, t1) {
1666 @Override
1667 public ODEEventHandler getHandler() {
1668 return (state, detector, increasing) -> {
1669 swap[0] = true;
1670 return super.getHandler().eventOccurred(state, detector, increasing);
1671 };
1672 }
1673 };
1674 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
1675 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
1676
1677 @Override
1678 public double g(ODEStateAndDerivative state) {
1679 if (swap[0]) {
1680 return detectorB.g(state);
1681 } else {
1682 return 1;
1683 }
1684 }
1685
1686 };
1687 ODEIntegrator integrator =
1688 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1689 integrator.addEventDetector(detectorA);
1690 integrator.addEventDetector(detectorC);
1691
1692
1693 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1694
1695
1696 Assert.assertEquals(2, events.size());
1697 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1698 Assert.assertEquals(true, events.get(0).isIncreasing());
1699 Assert.assertSame(detectorA, events.get(0).getDetector());
1700 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
1701 Assert.assertEquals(true, events.get(1).isIncreasing());
1702 Assert.assertSame(detectorC, events.get(1).getDetector());
1703 }
1704
1705
1706
1707
1708
1709
1710 @Test
1711 public void testEventChangesGFunctionDefinitionCancelReverse() {
1712
1713 double maxCheck = 5;
1714 double tolerance = 1e-6;
1715 double t1 = -11, t2 = -11.1;
1716
1717 List<Event> events = new ArrayList<>();
1718
1719 boolean[] swap = new boolean[1];
1720 final TimeDetector detectorA =
1721 new ContinuousDetector(maxCheck, tolerance, 100, Action.RESET_EVENTS, events, t1) {
1722 @Override
1723 public ODEEventHandler getHandler() {
1724 return (state, detector, increasing) -> {
1725 swap[0] = true;
1726 return super.getHandler().eventOccurred(state, detector, increasing);
1727 };
1728 }
1729 };
1730 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
1731 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
1732
1733 @Override
1734 public double g(ODEStateAndDerivative state) {
1735 if (!swap[0]) {
1736 return detectorB.g(state);
1737 } else {
1738 return 1;
1739 }
1740 }
1741
1742 };
1743 ODEIntegrator integrator =
1744 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1745 integrator.addEventDetector(detectorA);
1746 integrator.addEventDetector(detectorC);
1747
1748
1749 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1750
1751
1752 Assert.assertEquals(1, events.size());
1753 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1754 Assert.assertEquals(true, events.get(0).isIncreasing());
1755 Assert.assertSame(detectorA, events.get(0).getDetector());
1756 }
1757
1758
1759
1760
1761
1762
1763
1764 @Test
1765 public void testEventChangesGFunctionDefinitionDelayReverse() {
1766
1767 double maxCheck = 5;
1768 double tolerance = 1e-6;
1769 double t1 = -11, t2 = -11.1, t3 = -11.2;
1770
1771 List<Event> events = new ArrayList<>();
1772
1773 boolean[] swap = new boolean[1];
1774 final TimeDetector detectorA =
1775 new ContinuousDetector(maxCheck, tolerance, 100, Action.RESET_EVENTS, events, t1) {
1776 @Override
1777 public ODEEventHandler getHandler() {
1778 return (state, detector, increasing) -> {
1779 swap[0] = true;
1780 return super.getHandler().eventOccurred(state, detector, increasing);
1781 };
1782 }
1783 };
1784 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
1785 final TimeDetector detectorD = new ContinuousDetector(maxCheck, tolerance, 100, events, t3);
1786 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
1787
1788 @Override
1789 public double g(ODEStateAndDerivative state) {
1790 if (!swap[0]) {
1791 return detectorB.g(state);
1792 } else {
1793 return detectorD.g(state);
1794 }
1795 }
1796
1797 };
1798 ODEIntegrator integrator =
1799 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1800 integrator.addEventDetector(detectorA);
1801 integrator.addEventDetector(detectorC);
1802
1803
1804 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1805
1806
1807 Assert.assertEquals(2, events.size());
1808 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1809 Assert.assertEquals(true, events.get(0).isIncreasing());
1810 Assert.assertSame(detectorA, events.get(0).getDetector());
1811 Assert.assertEquals(t3, events.get(1).getT(), tolerance);
1812 Assert.assertEquals(true, events.get(1).isIncreasing());
1813 Assert.assertSame(detectorC, events.get(1).getDetector());
1814 }
1815
1816
1817
1818
1819
1820
1821 @Test
1822 public void testEventChangesGFunctionDefinitionAccelerateReverse() {
1823
1824 double maxCheck = 5;
1825 double tolerance = 1e-6;
1826 double t1 = -11, t2 = -11.1, t3 = -11.2;
1827
1828 List<Event> events = new ArrayList<>();
1829
1830 boolean[] swap = new boolean[1];
1831 final TimeDetector detectorA =
1832 new ContinuousDetector(maxCheck, tolerance, 100, Action.RESET_EVENTS, events, t1) {
1833 @Override
1834 public ODEEventHandler getHandler() {
1835 return (state, detector, increasing) -> {
1836 swap[0] = true;
1837 return super.getHandler().eventOccurred(state, detector, increasing);
1838 };
1839 }
1840 };
1841 final TimeDetector detectorB = new ContinuousDetector(maxCheck, tolerance, 100, events, t2);
1842 final TimeDetector detectorD = new ContinuousDetector(maxCheck, tolerance, 100, events, t3);
1843 BaseDetector detectorC = new BaseDetector(maxCheck, tolerance, 100, Action.CONTINUE, events) {
1844
1845 @Override
1846 public double g(ODEStateAndDerivative state) {
1847 if (swap[0]) {
1848 return detectorB.g(state);
1849 } else {
1850 return detectorD.g(state);
1851 }
1852 }
1853
1854 };
1855 ODEIntegrator integrator =
1856 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1857 integrator.addEventDetector(detectorA);
1858 integrator.addEventDetector(detectorC);
1859
1860
1861 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1862
1863
1864 Assert.assertEquals(2, events.size());
1865 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1866 Assert.assertEquals(true, events.get(0).isIncreasing());
1867 Assert.assertSame(detectorA, events.get(0).getDetector());
1868 Assert.assertEquals(t2, events.get(1).getT(), tolerance);
1869 Assert.assertEquals(true, events.get(1).isIncreasing());
1870 Assert.assertSame(detectorC, events.get(1).getDetector());
1871 }
1872
1873
1874 @Test
1875 public void testToleranceStopReverse() {
1876
1877 double maxCheck = 10;
1878 double tolerance = 1e-18;
1879 double t1 = -15.1;
1880
1881 List<Event> events = new ArrayList<>();
1882 TimeDetector detectorA = new FlatDetector(maxCheck, tolerance, 100, Action.STOP, events, t1);
1883 ODEIntegrator integrator =
1884 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1885 integrator.addEventDetector(detectorA);
1886
1887
1888 ODEStateAndDerivative finalState =
1889 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1890
1891
1892 Assert.assertEquals(1, events.size());
1893
1894 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1895 Assert.assertEquals(true, events.get(0).isIncreasing());
1896 Assert.assertSame(detectorA, events.get(0).getDetector());
1897 Assert.assertEquals(t1, finalState.getTime(), tolerance);
1898
1899
1900 finalState = integrator.integrate(new Equation(), finalState, -30.0);
1901
1902
1903 Assert.assertEquals(-30.0, finalState.getTime(), 0.0);
1904 }
1905
1906
1907
1908
1909
1910 @Test
1911 public void testLongInitialZeroReverse() {
1912
1913 double maxCheck = 10;
1914 double tolerance = 1;
1915
1916 List<Event> events = new ArrayList<>();
1917 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, Action.STOP, events, 50) {
1918 @Override
1919 public double g(ODEStateAndDerivative state) {
1920 if (state.getTime() > -2) {
1921 return 0;
1922 } else {
1923 return super.g(state);
1924 }
1925 }
1926 };
1927 ODEIntegrator integrator =
1928 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1929 integrator.addEventDetector(detectorA);
1930
1931
1932 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1933
1934
1935 Assert.assertEquals(0, events.size());
1936 }
1937
1938
1939
1940
1941
1942
1943 @Test
1944 public void testShortBracketingIntervalReverse() {
1945
1946 double maxCheck = 10;
1947 double tolerance = 1e-6;
1948 final double t1 = FastMath.nextDown(-10.0), t2 = -10.5;
1949
1950 List<Event> events = new ArrayList<>();
1951
1952 TimeDetector detectorA = new TimeDetector(maxCheck, tolerance, 100, events) {
1953 @Override
1954 public double g(ODEStateAndDerivative state) {
1955 final double t = state.getTime();
1956 if (t > t1) {
1957 return -1;
1958 } else if (t > t2) {
1959 return 1;
1960 } else {
1961 return -1;
1962 }
1963 }
1964 };
1965 TimeDetector detectorB = new TimeDetector(maxCheck, tolerance, 100, events, t1);
1966 ODEIntegrator integrator =
1967 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1968 integrator.addEventDetector(detectorA);
1969 integrator.addEventDetector(detectorB);
1970
1971
1972 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -30.0);
1973
1974
1975 Assert.assertEquals(3, events.size());
1976 Assert.assertEquals(t1, events.get(0).getT(), tolerance);
1977 Assert.assertEquals(false, events.get(0).isIncreasing());
1978 Assert.assertSame(detectorA, events.get(0).getDetector());
1979 Assert.assertEquals(t1, events.get(1).getT(), tolerance);
1980 Assert.assertEquals(true, events.get(1).isIncreasing());
1981 Assert.assertSame(detectorB, events.get(1).getDetector());
1982 Assert.assertEquals(t2, events.get(2).getT(), tolerance);
1983 Assert.assertEquals(true, events.get(2).isIncreasing());
1984 Assert.assertSame(detectorA, events.get(2).getDetector());
1985 }
1986
1987
1988 @Test
1989 public void testEventStepHandlerReverse() {
1990
1991 double tolerance = 1e-18;
1992 ODEIntegrator integrator =
1993 new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
1994 integrator.addEventDetector(new TimeDetector(100, tolerance, 100, -5));
1995 StepHandler stepHandler = new StepHandler();
1996 integrator.addStepHandler(stepHandler);
1997
1998
1999 ODEStateAndDerivative finalState = integrator
2000 .integrate(new Equation(), new ODEState(0, new double[2]), -10);
2001
2002
2003 Assert.assertEquals(-10.0, finalState.getTime(), tolerance);
2004 Assert.assertEquals(0.0,
2005 stepHandler.initialState.getTime(), tolerance);
2006 Assert.assertEquals(-10.0, stepHandler.finalTime, tolerance);
2007 Assert.assertEquals(-10.0,
2008 stepHandler.finalState.getTime(), tolerance);
2009 ODEStateInterpolator interpolator = stepHandler.interpolators.get(0);
2010 Assert.assertEquals(0.0,
2011 interpolator.getPreviousState().getTime(), tolerance);
2012 Assert.assertEquals(-5.0,
2013 interpolator.getCurrentState().getTime(), tolerance);
2014 interpolator = stepHandler.interpolators.get(1);
2015 Assert.assertEquals(-5.0,
2016 interpolator.getPreviousState().getTime(), tolerance);
2017 Assert.assertEquals(-10.0,
2018 interpolator.getCurrentState().getTime(), tolerance);
2019 Assert.assertEquals(2, stepHandler.interpolators.size());
2020 }
2021
2022
2023 @Test
2024 public void testEventCausedByDerivativesResetReverse() {
2025
2026 TimeDetector detectorA = new TimeDetector(10, 1e-6, 100, Action.RESET_STATE, -15.0) {
2027 @Override
2028 public ODEEventHandler getHandler() {
2029 return new ODEEventHandler() {
2030 @Override
2031 public Action eventOccurred(ODEStateAndDerivative state,
2032 ODEEventDetector detector, boolean increasing) {
2033 return Action.RESET_STATE;
2034 }
2035 @Override
2036 public ODEState resetState(ODEEventDetector detector, ODEStateAndDerivative state) {
2037 return null;
2038 }
2039 };
2040 }
2041 };
2042 ODEIntegrator integrator = new DormandPrince853Integrator(10, 10, 1e-7, 1e-7);
2043 integrator.addEventDetector(detectorA);
2044
2045 try {
2046
2047 integrator.integrate(new Equation(), new ODEState(0, new double[2]), -20.0);
2048 Assert.fail("Expected Exception");
2049 } catch (NullPointerException e) {
2050
2051 }
2052 }
2053
2054 @Test
2055 public void testResetChangesSignReverse() {
2056 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
2057 public int getDimension() { return 1; }
2058 public double[] computeDerivatives(double t, double[] y) { return new double[] { 1.0 }; }
2059 };
2060
2061 LutherIntegrator integrator = new LutherIntegrator(20.0);
2062 final double small = 1.0e-10;
2063 ResetChangesSignGenerator eventsGenerator = new ResetChangesSignGenerator(-6.0, -9.0, +0.5 * small, 8.0, small, 1000);
2064 integrator.addEventDetector(eventsGenerator);
2065 final ODEStateAndDerivative end = integrator.integrate(equation, new ODEState(0.0, new double[1]), -100.0);
2066 Assert.assertEquals(2, eventsGenerator.getCount());
2067 Assert.assertEquals(-9.0, end.getCompleteState()[0], 1.0e-12);
2068 Assert.assertEquals(-9.0 - 0.5 * small, end.getTime(), 1.0e-12);
2069 }
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079 private ODEStateAndDerivative state(double t) {
2080 return new ODEStateAndDerivative(t, new double[0], new double[0]);
2081 }
2082
2083
2084 private static abstract class BaseDetector implements ODEEventDetector {
2085
2086 private final AdaptableInterval maxCheck;
2087 private final int maxIter;
2088 private final BracketingNthOrderBrentSolver solver;
2089 protected final Action action;
2090
2091
2092 private final List<Event> events;
2093
2094 public BaseDetector(final double maxCheck, final double threshold, final int maxIter,
2095 Action action, List<Event> events) {
2096 this.maxCheck = s -> maxCheck;
2097 this.maxIter = maxIter;
2098 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
2099 this.action = action;
2100 this.events = events;
2101 }
2102
2103 public AdaptableInterval getMaxCheckInterval() {
2104 return maxCheck;
2105 }
2106
2107 public int getMaxIterationCount() {
2108 return maxIter;
2109 }
2110
2111 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
2112 return solver;
2113 }
2114
2115
2116
2117
2118
2119
2120 public List<Event> getEvents() {
2121 return events;
2122 }
2123
2124 @Override
2125 public ODEEventHandler getHandler() {
2126 return new ODEEventHandler() {
2127 @Override
2128 public Action eventOccurred(ODEStateAndDerivative state,
2129 ODEEventDetector detector, boolean increasing) {
2130 events.add(new Event(state, detector, increasing));
2131 return action;
2132 }
2133 };
2134 }
2135
2136 }
2137
2138
2139 private static class TimeDetector extends BaseDetector implements ODEEventDetector {
2140
2141
2142 protected final double[] eventTs;
2143
2144
2145
2146
2147
2148
2149 public TimeDetector(final double maxCheck, final double threshold, final int maxIter,
2150 double... eventTs) {
2151 this(maxCheck, threshold, maxIter, Action.CONTINUE, eventTs);
2152 }
2153
2154 public TimeDetector(final double maxCheck, final double threshold, final int maxIter,
2155 List<Event> events, double... eventTs) {
2156 this(maxCheck, threshold, maxIter, Action.CONTINUE, events, eventTs);
2157 }
2158
2159 public TimeDetector(final double maxCheck, final double threshold, final int maxIter,
2160 Action action, double... eventTs) {
2161 this(maxCheck, threshold, maxIter, action, new ArrayList<Event>(), eventTs);
2162 }
2163
2164 public TimeDetector(final double maxCheck, final double threshold, final int maxIter,
2165 Action action, List<Event> events, double... eventTs) {
2166 super(maxCheck, threshold, maxIter, action, events);
2167 this.eventTs = eventTs.clone();
2168 Arrays.sort(this.eventTs);
2169 }
2170
2171 @Override
2172 public double g(final ODEStateAndDerivative state) {
2173 final double t = state.getTime();
2174 int i = 0;
2175 while (i < eventTs.length && t > eventTs[i]) {
2176 i++;
2177 }
2178 i--;
2179 if (i < 0) {
2180 return t - eventTs[0];
2181 } else {
2182 int sign = (i % 2) * 2 - 1;
2183 return -sign * (t - eventTs[i]);
2184 }
2185 }
2186
2187 }
2188
2189 private static class Event {
2190
2191 private final ODEStateAndDerivative state;
2192 private final boolean increasing;
2193 private final ODEEventDetector detector;
2194
2195 public Event(ODEStateAndDerivative state, ODEEventDetector detector, boolean increasing) {
2196 this.increasing = increasing;
2197 this.state = state;
2198 this.detector = detector;
2199 }
2200
2201 public boolean isIncreasing() {
2202 return increasing;
2203 }
2204
2205 public double getT() {
2206 return state.getTime();
2207 }
2208
2209 public ODEEventDetector getDetector() {
2210 return detector;
2211 }
2212
2213 @Override
2214 public String toString() {
2215 return "Event{" +
2216 "increasing=" + increasing +
2217 ", time=" + state.getTime() +
2218 ", state=" + Arrays.toString(state.getCompleteState()) +
2219 '}';
2220 }
2221 }
2222
2223
2224
2225
2226
2227 private static class FlatDetector extends TimeDetector {
2228
2229 public FlatDetector(double maxCheck, double threshold, int maxIter,
2230 double... eventTs) {
2231 super(maxCheck, threshold, maxIter, eventTs);
2232 }
2233
2234 public FlatDetector(double maxCheck, double threshold, int maxIter,
2235 List<Event> events, double... eventTs) {
2236 super(maxCheck, threshold, maxIter, events, eventTs);
2237 }
2238
2239 public FlatDetector(double maxCheck, double threshold, int maxIter,
2240 Action action, List<Event> events, double... eventTs) {
2241 super(maxCheck, threshold, maxIter, action, events, eventTs);
2242 }
2243
2244 @Override
2245 public double g(ODEStateAndDerivative state) {
2246 final double g = super.g(state);
2247 return FastMath.signum(g);
2248 }
2249
2250 }
2251
2252
2253 private static class ContinuousDetector extends TimeDetector {
2254
2255 public ContinuousDetector(double maxCheck, double threshold, int maxIter,
2256 List<Event> events, double... eventTs) {
2257 super(maxCheck, threshold, maxIter, events, eventTs);
2258 }
2259
2260 public ContinuousDetector(double maxCheck, double threshold, int maxIter,
2261 Action action, List<Event> events, double... eventTs) {
2262 super(maxCheck, threshold, maxIter, action, events, eventTs);
2263 }
2264
2265 @Override
2266 public double g(ODEStateAndDerivative state) {
2267 final double t = state.getTime();
2268 int i = 0;
2269 while (i < eventTs.length && t > eventTs[i]) {
2270 i++;
2271 }
2272 i--;
2273 if (i < 0) {
2274 return t - eventTs[0];
2275 } else if (i < eventTs.length - 1) {
2276 int sign = (i % 2) * 2 - 1;
2277 return -sign * (t - eventTs[i]) * (eventTs[i + 1] - t);
2278 } else {
2279 int sign = (i % 2) * 2 - 1;
2280 return -sign * (t - eventTs[i]);
2281 }
2282 }
2283 }
2284
2285
2286 private static class ResetDetector extends TimeDetector {
2287
2288 private final ODEState resetState;
2289
2290 public ResetDetector(double maxCheck, double threshold, int maxIter,
2291 List<Event> events, ODEState state, double eventT) {
2292 super(maxCheck, threshold, maxIter, Action.RESET_STATE, events, eventT);
2293 this.resetState = state;
2294 }
2295
2296 @Override
2297 public ODEEventHandler getHandler() {
2298 return new ODEEventHandler() {
2299 @Override
2300 public Action eventOccurred(ODEStateAndDerivative state,
2301 ODEEventDetector detector, boolean increasing) {
2302 return ResetDetector.super.getHandler().eventOccurred(state, detector, increasing);
2303 }
2304 @Override
2305 public ODEState resetState(ODEEventDetector detector, ODEStateAndDerivative state) {
2306 Assert.assertEquals(eventTs[0], state.getTime(), 0);
2307 return resetState;
2308 }
2309 };
2310 }
2311
2312 }
2313
2314
2315 private static class StateDetector extends BaseDetector {
2316
2317 private final double triggerState;
2318
2319 public StateDetector(double maxCheck, double threshold, int maxIter,
2320 List<Event> events, double triggerState) {
2321 super(maxCheck, threshold, maxIter, Action.CONTINUE, events);
2322 this.triggerState = triggerState;
2323 }
2324
2325 @Override
2326 public double g(ODEStateAndDerivative state) {
2327 return state.getPrimaryState()[0] - this.triggerState;
2328 }
2329 }
2330
2331
2332 public static class Equation implements OrdinaryDifferentialEquation {
2333
2334 public int getDimension() {
2335 return 2;
2336 }
2337
2338 public double[] computeDerivatives(double t, double[] y) {
2339 return new double[]{1.0, 2.0};
2340 }
2341
2342 }
2343
2344 private static class StepHandler implements ODEStepHandler {
2345
2346 private ODEStateAndDerivative initialState;
2347 private double finalTime;
2348 private List<ODEStateInterpolator> interpolators = new ArrayList<>();
2349 private ODEStateAndDerivative finalState;
2350
2351 @Override
2352 public void init(ODEStateAndDerivative initialState, double finalTime) {
2353 this.initialState = initialState;
2354 this.finalTime = finalTime;
2355 }
2356
2357 @Override
2358 public void handleStep(ODEStateInterpolator interpolator) {
2359 this.interpolators.add(interpolator);
2360 }
2361
2362 @Override
2363 public void finish(ODEStateAndDerivative finalState) {
2364 this.finalState = finalState;
2365 }
2366 }
2367
2368 private class ResetChangesSignGenerator implements ODEEventDetector {
2369
2370 private final AdaptableInterval maxCheck;
2371 private final int maxIter;
2372 private final BracketingNthOrderBrentSolver solver;
2373 final double y1;
2374 final double y2;
2375 final double change;
2376 int count;
2377
2378 public ResetChangesSignGenerator(final double y1, final double y2, final double change,
2379 final double maxCheck, final double threshold, final int maxIter) {
2380 this.maxCheck = s -> maxCheck;
2381 this.maxIter = maxIter;
2382 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
2383 this.y1 = y1;
2384 this.y2 = y2;
2385 this.change = change;
2386 this.count = 0;
2387 }
2388
2389 public AdaptableInterval getMaxCheckInterval() {
2390 return maxCheck;
2391 }
2392
2393 public int getMaxIterationCount() {
2394 return maxIter;
2395 }
2396
2397 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
2398 return solver;
2399 }
2400
2401 public ODEEventHandler getHandler() {
2402 return new ODEEventHandler() {
2403 public Action eventOccurred(ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) {
2404 return ++count < 2 ? Action.RESET_STATE : Action.STOP;
2405 }
2406
2407 public ODEState resetState(ODEEventDetector detector, ODEStateAndDerivative s) {
2408 return new ODEState(s.getTime(), new double[] { s.getCompleteState()[0] + change });
2409 }
2410 };
2411 }
2412
2413 public double g(ODEStateAndDerivative s) {
2414 return (s.getCompleteState()[0] - y1) * (s.getCompleteState()[0] - y2);
2415 }
2416
2417 public int getCount() {
2418 return count;
2419 }
2420
2421 }
2422
2423 }