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