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