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 java.util.ArrayList;
25 import java.util.List;
26
27 import org.hipparchus.analysis.UnivariateFunction;
28 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
29 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
30 import org.hipparchus.exception.MathIllegalArgumentException;
31 import org.hipparchus.exception.MathIllegalStateException;
32 import org.hipparchus.ode.ODEIntegrator;
33 import org.hipparchus.ode.ODEState;
34 import org.hipparchus.ode.ODEStateAndDerivative;
35 import org.hipparchus.ode.OrdinaryDifferentialEquation;
36 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
37 import org.junit.Assert;
38 import org.junit.Test;
39
40 /** Tests for overlapping state events. Also tests an event function that does
41 * not converge to zero, but does have values of opposite sign around its root.
42 */
43 public class OverlappingEventsTest implements OrdinaryDifferentialEquation {
44
45 /** Expected event times for first event. */
46 private static final double[] EVENT_TIMES1 = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0,
47 7.0, 8.0, 9.0};
48
49 /** Expected event times for second event. */
50 private static final double[] EVENT_TIMES2 = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0,
51 3.5, 4.0, 4.5, 5.0, 5.5, 6.0,
52 6.5, 7.0, 7.5, 8.0, 8.5, 9.0,
53 9.5};
54
55 /** Test for events that occur at the exact same time, but due to numerical
56 * calculations occur very close together instead. Uses event type 0. See
57 * {@link org.hipparchus.ode.events.ODEEventDetector#g(org.hipparchus.ode.ODEStateAndDerivative)}
58 * ODEEventDetector.g(stateAndDerivative)}.
59 */
60 @Test
61 public void testOverlappingEvents0()
62 throws MathIllegalArgumentException, MathIllegalStateException {
63 test(0);
64 }
65
66 /** Test for events that occur at the exact same time, but due to numerical
67 * calculations occur very close together instead. Uses event type 1. See
68 * {@link org.hipparchus.ode.events.ODEEventDetector#g(org.hipparchus.ode.ODEStateAndDerivative)}
69 * * ODEEventDetector.g(stateAndDerivative)}.
70 */
71 @Test
72 public void testOverlappingEvents1()
73 throws MathIllegalArgumentException, MathIllegalStateException {
74 test(1);
75 }
76
77 /** Test for events that occur at the exact same time, but due to numerical
78 * calculations occur very close together instead.
79 * @param eventType the type of events to use. See
80 * {@link org.hipparchus.ode.events.ODEEventDetector#g(org.hipparchus.ode.ODEStateAndDerivative)}
81 ODEEventDetector.g(stateAndDerivative)}.
82 */
83 public void test(int eventType)
84 throws MathIllegalArgumentException, MathIllegalStateException {
85 double e = 1e-15;
86 ODEIntegrator integrator = new DormandPrince853Integrator(e, 100.0, 1e-7, 1e-7);
87 ODEEventDetector evt1 = new Event(0.1, e, 999, 0, eventType);
88 ODEEventDetector evt2 = new Event(0.1, e, 999, 1, eventType);
89 integrator.addEventDetector(evt1);
90 integrator.addEventDetector(evt2);
91 double t = 0.0;
92 double tEnd = 9.75;
93 double[] y = {0.0, 0.0};
94 List<Double> events1 = new ArrayList<Double>();
95 List<Double> events2 = new ArrayList<Double>();
96 while (t < tEnd) {
97 final ODEStateAndDerivative finalState =
98 integrator.integrate(this, new ODEState(t, y), tEnd);
99 t = finalState.getTime();
100 y = finalState.getPrimaryState();
101 //System.out.println("t=" + t + ",\t\ty=[" + y[0] + "," + y[1] + "]");
102 if (y[0] >= 1.0) {
103 y[0] = 0.0;
104 events1.add(t);
105 //System.out.println("Event 1 @ t=" + t);
106 }
107 if (y[1] >= 1.0) {
108 y[1] = 0.0;
109 events2.add(t);
110 //System.out.println("Event 2 @ t=" + t);
111 }
112 }
113 Assert.assertEquals(EVENT_TIMES1.length, events1.size());
114 Assert.assertEquals(EVENT_TIMES2.length, events2.size());
115 for(int i = 0; i < EVENT_TIMES1.length; i++) {
116 Assert.assertEquals(EVENT_TIMES1[i], events1.get(i), 1e-7);
117 }
118 for(int i = 0; i < EVENT_TIMES2.length; i++) {
119 Assert.assertEquals(EVENT_TIMES2[i], events2.get(i), 1e-7);
120 }
121 //System.out.println();
122 }
123
124 /** {@inheritDoc} */
125 public int getDimension() {
126 return 2;
127 }
128
129 /** {@inheritDoc} */
130 public double[] computeDerivatives(double t, double[] y) {
131 return new double[] { 1.0, 2.0 };
132 }
133
134 /** State events for this unit test. */
135 private class Event implements ODEEventDetector {
136
137 private final AdaptableInterval maxCheck;
138 private final int maxIter;
139 private final BracketingNthOrderBrentSolver solver;
140 private final int idx;
141 private final int eventType;
142
143 /** Constructor for the {@link Event} class.
144 * @param maxCheck maximum checking interval, must be strictly positive (s)
145 * @param threshold convergence threshold (s)
146 * @param maxIter maximum number of iterations in the event time search
147 * @param idx the index of the continuous variable to use
148 * @param eventType the type of event to use. See {@link #g}
149 */
150 public Event(final double maxCheck, final double threshold, final int maxIter,
151 final int idx, final int eventType) {
152 this.maxCheck = s -> maxCheck;
153 this.maxIter = maxIter;
154 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
155 this.idx = idx;
156 this.eventType = eventType;
157 }
158
159 public AdaptableInterval getMaxCheckInterval() {
160 return maxCheck;
161 }
162
163 public int getMaxIterationCount() {
164 return maxIter;
165 }
166
167 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
168 return solver;
169 }
170
171 public ODEEventHandler getHandler() {
172 return (state, detector, increasing) -> Action.STOP;
173 }
174
175 /** {@inheritDoc} */
176 public double g(ODEStateAndDerivative s) {
177 return (eventType == 0) ? s.getPrimaryState()[idx] >= 1.0 ? 1.0 : -1.0
178 : s.getPrimaryState()[idx] - 1.0;
179 }
180
181 }
182
183 }