1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.ode.events;
24
25
26 import org.hipparchus.analysis.UnivariateFunction;
27 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
28 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
29 import org.hipparchus.exception.MathIllegalArgumentException;
30 import org.hipparchus.exception.MathIllegalStateException;
31 import org.hipparchus.ode.EquationsMapper;
32 import org.hipparchus.ode.ExpandableODE;
33 import org.hipparchus.ode.ODEState;
34 import org.hipparchus.ode.ODEStateAndDerivative;
35 import org.hipparchus.ode.OrdinaryDifferentialEquation;
36 import org.hipparchus.ode.SecondaryODE;
37 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
38 import org.hipparchus.ode.nonstiff.LutherIntegrator;
39 import org.hipparchus.ode.sampling.DummyStepInterpolator;
40 import org.hipparchus.ode.sampling.ODEStateInterpolator;
41 import org.junit.Assert;
42 import org.junit.Test;
43
44 public class DetectorBasedEventStateTest {
45
46
47 @Test
48 public void closeEvents() throws MathIllegalArgumentException, MathIllegalStateException {
49
50 final double r1 = 90.0;
51 final double r2 = 135.0;
52 final double gap = r2 - r1;
53
54 final double tolerance = 0.1;
55 DetectorBasedEventState es = new DetectorBasedEventState(new CloseEventsGenerator(r1, r2, 1.5 * gap, tolerance, 100));
56 EquationsMapper mapper = new ExpandableODE(new OrdinaryDifferentialEquation() {
57 @Override
58 public int getDimension() {
59 return 0;
60 }
61 @Override
62 public double[] computeDerivatives(double t, double[] y) {
63 return new double[0];
64 }
65 }).getMapper();
66
67 double[] a = new double[0];
68 ODEStateAndDerivative osdLongBefore = new ODEStateAndDerivative(r1 - 1.5 * gap, a, a);
69 ODEStateAndDerivative osBefore = new ODEStateAndDerivative(r1 - 0.5 * gap, a, a);
70 ODEStateInterpolator interpolatorA = new DummyStepInterpolator(true,
71 osdLongBefore, osBefore,
72 osdLongBefore, osBefore,
73 mapper);
74 es.reinitializeBegin(interpolatorA);
75 Assert.assertFalse(es.evaluateStep(interpolatorA));
76
77 ODEStateAndDerivative osdBetween = new ODEStateAndDerivative(0.5 * (r1 + r2), a, a);
78 ODEStateInterpolator interpolatorB = new DummyStepInterpolator(true,
79 osBefore, osdBetween,
80 osBefore, osdBetween,
81 mapper);
82 Assert.assertTrue(es.evaluateStep(interpolatorB));
83 Assert.assertEquals(r1, es.getEventTime(), tolerance);
84 ODEStateAndDerivative osdAtEvent = new ODEStateAndDerivative(es.getEventTime(), a, a);
85 es.doEvent(osdAtEvent);
86
87 ODEStateAndDerivative osdAfterSecond = new ODEStateAndDerivative(r2 + 0.4 * gap, a, a);
88 ODEStateInterpolator interpolatorC = new DummyStepInterpolator(true,
89 osdAtEvent, osdAfterSecond,
90 osdAtEvent, osdAfterSecond,
91 mapper);
92 Assert.assertTrue(es.evaluateStep(interpolatorC));
93 Assert.assertEquals(r2, es.getEventTime(), tolerance);
94
95 }
96
97
98 @Test
99 public void testIssue695()
100 throws MathIllegalArgumentException, MathIllegalStateException {
101
102 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
103 public int getDimension() {
104 return 1;
105 }
106 public double[] computeDerivatives(double t, double[] y) {
107 return new double[] { 1.0 };
108 }
109 };
110
111 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
112 integrator.addEventDetector(new ResettingEvent(10.99, 0.1, 1.0e-9, 1000));
113 integrator.addEventDetector(new ResettingEvent(11.01, 0.1, 1.0e-9, 1000));
114 integrator.setInitialStepSize(3.0);
115
116 double target = 30.0;
117 ODEStateAndDerivative finalState =
118 integrator.integrate(equation, new ODEState(0.0, new double[1]), target);
119 Assert.assertEquals(target, finalState.getTime(), 1.0e-10);
120 Assert.assertEquals(32.0, finalState.getPrimaryState()[0], 1.0e-10);
121
122 }
123
124 private static class ResettingEvent implements ODEEventDetector {
125
126 private static double lastTriggerTime = Double.NEGATIVE_INFINITY;
127 private final AdaptableInterval maxCheck;
128 private final int maxIter;
129 private final BracketingNthOrderBrentSolver solver;
130 private final double tEvent;
131
132 public ResettingEvent(final double tEvent,
133 final double maxCheck, final double threshold, final int maxIter) {
134 this.maxCheck = s -> maxCheck;
135 this.maxIter = maxIter;
136 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
137 this.tEvent = tEvent;
138 }
139
140 public AdaptableInterval getMaxCheckInterval() {
141 return maxCheck;
142 }
143
144 public int getMaxIterationCount() {
145 return maxIter;
146 }
147
148 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
149 return solver;
150 }
151
152 public double g(ODEStateAndDerivative s) {
153
154
155
156
157 Assert.assertTrue("going backard in time! (" + s.getTime() + " < " + lastTriggerTime + ")",
158 s.getTime() >= lastTriggerTime);
159 return s.getTime() - tEvent;
160 }
161
162 public ODEEventHandler getHandler() {
163 return new ODEEventHandler() {
164 public Action eventOccurred(ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) {
165
166 lastTriggerTime = s.getTime();
167 return Action.RESET_STATE;
168 }
169
170 public ODEStateAndDerivative resetState(ODEEventDetector detector, ODEStateAndDerivative s) {
171 double[] y = s.getPrimaryState();
172 y[0] += 1.0;
173 return new ODEStateAndDerivative(s.getTime(), y, s.getPrimaryDerivative());
174 }
175 };
176 }
177
178 }
179
180
181 @Test
182 public void testIssue965()
183 throws MathIllegalArgumentException, MathIllegalStateException {
184
185 ExpandableODE equation = new ExpandableODE(new OrdinaryDifferentialEquation() {
186 public int getDimension() {
187 return 1;
188 }
189 public double[] computeDerivatives(double t, double[] y) {
190 return new double[] { 2.0 };
191 }
192 });
193 int index = equation.addSecondaryEquations(new SecondaryODE() {
194 public int getDimension() {
195 return 1;
196 }
197 public double[] computeDerivatives(double t, double[] primary,
198 double[] primaryDot, double[] secondary) {
199 return new double[] { -3.0 };
200 }
201 });
202 Assert.assertEquals(1, index);
203
204 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
205 integrator.addEventDetector(new SecondaryStateEvent(index, -3.0, 0.1, 1.0e-9, 1000));
206 integrator.setInitialStepSize(3.0);
207
208 ODEState initialState = new ODEState(0.0,
209 new double[] { 0.0 },
210 new double[][] { { 0.0 } });
211 ODEStateAndDerivative finalState = integrator.integrate(equation, initialState, 30.0);
212 Assert.assertEquals( 1.0, finalState.getTime(), 1.0e-10);
213 Assert.assertEquals( 2.0, finalState.getPrimaryState()[0], 1.0e-10);
214 Assert.assertEquals(-3.0, finalState.getSecondaryState(index)[0], 1.0e-10);
215
216 }
217
218 private static class SecondaryStateEvent implements ODEEventDetector {
219
220 private final AdaptableInterval maxCheck;
221 private final int maxIter;
222 private final BracketingNthOrderBrentSolver solver;
223 private int index;
224 private final double target;
225
226 public SecondaryStateEvent(final int index, final double target,
227 final double maxCheck, final double threshold, final int maxIter) {
228 this.maxCheck = s -> maxCheck;
229 this.maxIter = maxIter;
230 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
231 this.index = index;
232 this.target = target;
233 }
234
235 public AdaptableInterval getMaxCheckInterval() {
236 return maxCheck;
237 }
238
239 public int getMaxIterationCount() {
240 return maxIter;
241 }
242
243 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
244 return solver;
245 }
246
247
248 public ODEEventHandler getHandler() {
249 return (state, detector, increasing) -> Action.STOP;
250 }
251
252 public double g(ODEStateAndDerivative s) {
253 return s.getSecondaryState(index)[0] - target;
254 }
255
256 }
257
258 @Test
259 public void testEventsCloserThanThreshold()
260 throws MathIllegalArgumentException, MathIllegalStateException {
261
262 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
263
264 public int getDimension() {
265 return 1;
266 }
267
268 public double[] computeDerivatives(double t, double[] y) {
269 return new double[] { 1.0 };
270 }
271 };
272
273 LutherIntegrator integrator = new LutherIntegrator(20.0);
274 CloseEventsGenerator eventsGenerator =
275 new CloseEventsGenerator(9.0 - 1.0 / 128, 9.0 + 1.0 / 128, 1.0, 0.02, 1000);
276 integrator.addEventDetector(eventsGenerator);
277 double tEnd = integrator.integrate(equation, new ODEState(0.0, new double[1]), 100.0).getTime();
278 Assert.assertEquals( 2, eventsGenerator.getCount());
279 Assert.assertEquals( 9.0 + 1.0 / 128, tEnd, 1.0 / 32.0);
280
281 }
282
283 private class CloseEventsGenerator implements ODEEventDetector {
284
285 private final AdaptableInterval maxCheck;
286 private final int maxIter;
287 private final BracketingNthOrderBrentSolver solver;
288 final double r1;
289 final double r2;
290 int count;
291
292 public CloseEventsGenerator(final double r1, final double r2,
293 final double maxCheck, final double threshold, final int maxIter) {
294 this.maxCheck = s -> maxCheck;
295 this.maxIter = maxIter;
296 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
297 this.r1 = r1;
298 this.r2 = r2;
299 this.count = 0;
300 }
301
302 public AdaptableInterval getMaxCheckInterval() {
303 return maxCheck;
304 }
305
306 public int getMaxIterationCount() {
307 return maxIter;
308 }
309
310 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
311 return solver;
312 }
313
314 public double g(ODEStateAndDerivative s) {
315 return (s.getTime() - r1) * (r2 - s.getTime());
316 }
317
318 public ODEEventHandler getHandler() {
319 return (state, detector, increasing) -> ++count < 2 ? Action.CONTINUE : Action.STOP;
320 }
321
322 public int getCount() {
323 return count;
324 }
325
326 }
327
328 }