1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.ode.events;
23
24 import org.hipparchus.analysis.UnivariateFunction;
25 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
26 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.ode.ODEIntegrator;
30 import org.hipparchus.ode.ODEState;
31 import org.hipparchus.ode.ODEStateAndDerivative;
32 import org.hipparchus.ode.OrdinaryDifferentialEquation;
33 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
34 import org.junit.jupiter.api.Test;
35
36 import java.util.ArrayList;
37 import java.util.List;
38
39 import static org.junit.jupiter.api.Assertions.assertEquals;
40
41 /** Tests for overlapping state events. Also tests an event function that does
42 * not converge to zero, but does have values of opposite sign around its root.
43 */
44 public class OverlappingEventsTest implements OrdinaryDifferentialEquation {
45
46 /** Expected event times for first event. */
47 private static final double[] EVENT_TIMES1 = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0,
48 7.0, 8.0, 9.0};
49
50 /** Expected event times for second event. */
51 private static final double[] EVENT_TIMES2 = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0,
52 3.5, 4.0, 4.5, 5.0, 5.5, 6.0,
53 6.5, 7.0, 7.5, 8.0, 8.5, 9.0,
54 9.5};
55
56 /** Test for events that occur at the exact same time, but due to numerical
57 * calculations occur very close together instead. Uses event type 0. See
58 * {@link org.hipparchus.ode.events.ODEEventDetector#g(org.hipparchus.ode.ODEStateAndDerivative)}
59 * ODEEventDetector.g(stateAndDerivative)}.
60 */
61 @Test
62 void testOverlappingEvents0()
63 throws MathIllegalArgumentException, MathIllegalStateException {
64 test(0);
65 }
66
67 /** Test for events that occur at the exact same time, but due to numerical
68 * calculations occur very close together instead. Uses event type 1. See
69 * {@link org.hipparchus.ode.events.ODEEventDetector#g(org.hipparchus.ode.ODEStateAndDerivative)}
70 * * ODEEventDetector.g(stateAndDerivative)}.
71 */
72 @Test
73 void testOverlappingEvents1()
74 throws MathIllegalArgumentException, MathIllegalStateException {
75 test(1);
76 }
77
78 /** Test for events that occur at the exact same time, but due to numerical
79 * calculations occur very close together instead.
80 * @param eventType the type of events to use. See
81 * {@link org.hipparchus.ode.events.ODEEventDetector#g(org.hipparchus.ode.ODEStateAndDerivative)}
82 ODEEventDetector.g(stateAndDerivative)}.
83 */
84 public void test(int eventType)
85 throws MathIllegalArgumentException, MathIllegalStateException {
86 double e = 1e-15;
87 ODEIntegrator integrator = new DormandPrince853Integrator(e, 100.0, 1e-7, 1e-7);
88 ODEEventDetector evt1 = new Event(0.1, e, 999, 0, eventType);
89 ODEEventDetector evt2 = new Event(0.1, e, 999, 1, eventType);
90 integrator.addEventDetector(evt1);
91 integrator.addEventDetector(evt2);
92 double t = 0.0;
93 double tEnd = 9.75;
94 double[] y = {0.0, 0.0};
95 List<Double> events1 = new ArrayList<Double>();
96 List<Double> events2 = new ArrayList<Double>();
97 while (t < tEnd) {
98 final ODEStateAndDerivative finalState =
99 integrator.integrate(this, new ODEState(t, y), tEnd);
100 t = finalState.getTime();
101 y = finalState.getPrimaryState();
102 //System.out.println("t=" + t + ",\t\ty=[" + y[0] + "," + y[1] + "]");
103 if (y[0] >= 1.0) {
104 y[0] = 0.0;
105 events1.add(t);
106 //System.out.println("Event 1 @ t=" + t);
107 }
108 if (y[1] >= 1.0) {
109 y[1] = 0.0;
110 events2.add(t);
111 //System.out.println("Event 2 @ t=" + t);
112 }
113 }
114 assertEquals(EVENT_TIMES1.length, events1.size());
115 assertEquals(EVENT_TIMES2.length, events2.size());
116 for(int i = 0; i < EVENT_TIMES1.length; i++) {
117 assertEquals(EVENT_TIMES1[i], events1.get(i), 1e-7);
118 }
119 for(int i = 0; i < EVENT_TIMES2.length; i++) {
120 assertEquals(EVENT_TIMES2[i], events2.get(i), 1e-7);
121 }
122 //System.out.println();
123 }
124
125 /** {@inheritDoc} */
126 public int getDimension() {
127 return 2;
128 }
129
130 /** {@inheritDoc} */
131 public double[] computeDerivatives(double t, double[] y) {
132 return new double[] { 1.0, 2.0 };
133 }
134
135 /** State events for this unit test. */
136 private class Event implements ODEEventDetector {
137
138 private final AdaptableInterval maxCheck;
139 private final int maxIter;
140 private final BracketingNthOrderBrentSolver solver;
141 private final int idx;
142 private final int eventType;
143
144 /** Constructor for the {@link Event} class.
145 * @param maxCheck maximum checking interval, must be strictly positive (s)
146 * @param threshold convergence threshold (s)
147 * @param maxIter maximum number of iterations in the event time search
148 * @param idx the index of the continuous variable to use
149 * @param eventType the type of event to use. See {@link #g}
150 */
151 public Event(final double maxCheck, final double threshold, final int maxIter,
152 final int idx, final int eventType) {
153 this.maxCheck = (s, isForward) -> maxCheck;
154 this.maxIter = maxIter;
155 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
156 this.idx = idx;
157 this.eventType = eventType;
158 }
159
160 public AdaptableInterval getMaxCheckInterval() {
161 return maxCheck;
162 }
163
164 public int getMaxIterationCount() {
165 return maxIter;
166 }
167
168 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
169 return solver;
170 }
171
172 public ODEEventHandler getHandler() {
173 return (state, detector, increasing) -> Action.STOP;
174 }
175
176 /** {@inheritDoc} */
177 public double g(ODEStateAndDerivative s) {
178 return (eventType == 0) ? s.getPrimaryState()[idx] >= 1.0 ? 1.0 : -1.0
179 : s.getPrimaryState()[idx] - 1.0;
180 }
181
182 }
183
184 }