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;
24
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.exception.MathIllegalStateException;
27 import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
28 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
29 import org.hipparchus.ode.nonstiff.EulerIntegrator;
30 import org.hipparchus.ode.sampling.DummyStepInterpolator;
31 import org.hipparchus.ode.sampling.ODEStateInterpolator;
32 import org.hipparchus.ode.sampling.ODEStepHandler;
33 import org.hipparchus.util.FastMath;
34 import org.junit.jupiter.api.AfterEach;
35 import org.junit.jupiter.api.BeforeEach;
36 import org.junit.jupiter.api.Test;
37
38 import java.io.ByteArrayInputStream;
39 import java.io.ByteArrayOutputStream;
40 import java.io.IOException;
41 import java.io.ObjectInputStream;
42 import java.io.ObjectOutputStream;
43 import java.util.Random;
44
45 import static org.junit.jupiter.api.Assertions.assertEquals;
46 import static org.junit.jupiter.api.Assertions.assertFalse;
47 import static org.junit.jupiter.api.Assertions.assertTrue;
48 import static org.junit.jupiter.api.Assertions.fail;
49
50 public class DenseOutputModelTest {
51
52 TestProblem3 pb;
53 ODEIntegrator integ;
54
55 @Test
56 void testBoundaries() throws MathIllegalArgumentException, MathIllegalStateException {
57 integ.addStepHandler(new DenseOutputModel());
58 integ.integrate(pb, pb.getInitialState(), pb.getFinalTime());
59 DenseOutputModel dom = (DenseOutputModel) integ.getStepHandlers().iterator().next();
60 double tBefore = 2.0 * pb.getInitialTime() - pb.getFinalTime();
61 assertEquals(tBefore, dom.getInterpolatedState(tBefore).getTime(), 1.0e-10);
62 double tAfter = 2.0 * pb.getFinalTime() - pb.getInitialTime();
63 assertEquals(tAfter, dom.getInterpolatedState(tAfter).getTime(), 1.0e-10);
64 double tMiddle = 2.0 * pb.getFinalTime() - pb.getInitialTime();
65 assertEquals(tMiddle, dom.getInterpolatedState(tMiddle).getTime(), 1.0e-10);
66 }
67
68 @Test
69 void testRandomAccess() throws MathIllegalArgumentException, MathIllegalStateException {
70
71 DenseOutputModel dom = new DenseOutputModel();
72 integ.addStepHandler(dom);
73 integ.integrate(pb, pb.getInitialState(), pb.getFinalTime());
74
75 Random random = new Random(347588535632l);
76 double maxError = 0.0;
77 double maxErrorDot = 0.0;
78 for (int i = 0; i < 1000; ++i) {
79 double r = random.nextDouble();
80 double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
81 ODEStateAndDerivative sd = dom.getInterpolatedState(time);
82 double[] interpolatedY = sd.getPrimaryState();
83 double[] interpolatedYDot = sd.getPrimaryDerivative();
84 double[] theoreticalY = pb.computeTheoreticalState(time);
85 double[] theoreticalYDot = pb.doComputeDerivatives(time, theoreticalY);
86 double dx = interpolatedY[0] - theoreticalY[0];
87 double dy = interpolatedY[1] - theoreticalY[1];
88 double error = dx * dx + dy * dy;
89 maxError = FastMath.max(maxError, error);
90 double dxDot = interpolatedYDot[0] - theoreticalYDot[0];
91 double dyDot = interpolatedYDot[1] - theoreticalYDot[1];
92 double errorDot = dxDot * dxDot + dyDot * dyDot;
93 maxErrorDot = FastMath.max(maxErrorDot, errorDot);
94 }
95
96 assertEquals(0.0, maxError, 1.0e-9);
97 assertEquals(0.0, maxErrorDot, 4.0e-7);
98
99 }
100
101 @Test
102 void testModelsMerging() throws MathIllegalArgumentException, MathIllegalStateException {
103
104
105 OrdinaryDifferentialEquation problem =
106 new OrdinaryDifferentialEquation() {
107 @Override
108 public double[] computeDerivatives(double t, double[] y) {
109 return new double[] { -y[1], y[0] };
110 }
111 @Override
112 public int getDimension() {
113 return 2;
114 }
115 };
116
117
118 DenseOutputModel dom1 = new DenseOutputModel();
119 ODEIntegrator integ1 = new DormandPrince853Integrator(0, 1.0, 1.0e-8, 1.0e-8);
120 integ1.addStepHandler(dom1);
121 integ1.integrate(problem, new ODEState(FastMath.PI, new double[] { -1.0, 0.0 }), 0);
122
123
124 DenseOutputModel dom2 = new DenseOutputModel();
125 ODEIntegrator integ2 = new DormandPrince853Integrator(0, 0.1, 1.0e-12, 1.0e-12);
126 integ2.addStepHandler(dom2);
127 integ2.integrate(problem, new ODEState(2.0 * FastMath.PI, new double[] { 1.0, 0.0 }), FastMath.PI);
128
129
130 DenseOutputModel dom = new DenseOutputModel();
131 dom.append(dom2);
132 dom.append(new DenseOutputModel());
133 dom.append(dom1);
134
135
136 assertEquals(2.0 * FastMath.PI, dom.getInitialTime(), 1.0e-12);
137 assertEquals(0, dom.getFinalTime(), 1.0e-12);
138 for (double t = 0; t < 2.0 * FastMath.PI; t += 0.1) {
139 final double[] y = dom.getInterpolatedState(t).getPrimaryState();
140 assertEquals(FastMath.cos(t), y[0], 1.0e-7);
141 assertEquals(FastMath.sin(t), y[1], 1.0e-7);
142 }
143
144 }
145
146 @Test
147 void testErrorConditions() throws MathIllegalArgumentException, MathIllegalStateException {
148
149 DenseOutputModel cm = new DenseOutputModel();
150 cm.handleStep(buildInterpolator(0, new double[] { 0.0, 1.0, -2.0 }, 1));
151
152
153 assertTrue(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0 }, 2.0));
154
155
156 assertTrue(checkAppendError(cm, 10.0, new double[] { 0.0, 1.0, -2.0 }, 20.0));
157
158
159 assertTrue(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0, -2.0 }, 0.0));
160
161
162 assertFalse(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0, -2.0 }, 2.0));
163
164 }
165
166 @Test
167 void testSerialization() {
168 try {
169 TestProblem1 pb = new TestProblem1();
170 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
171 EulerIntegrator integ = new EulerIntegrator(step);
172 integ.addStepHandler(new DenseOutputModel());
173 integ.integrate(pb, pb.getInitialState(), pb.getFinalTime());
174
175 ByteArrayOutputStream bos = new ByteArrayOutputStream();
176 ObjectOutputStream oos = new ObjectOutputStream(bos);
177 for (ODEStepHandler handler : integ.getStepHandlers()) {
178 oos.writeObject(handler);
179 }
180
181 int expectedSize = 131976;
182 assertTrue(bos.size () > 9 * expectedSize / 10, "size = " + bos.size());
183 assertTrue(bos.size () < 11 * expectedSize / 10, "size = " + bos.size());
184
185 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
186 ObjectInputStream ois = new ObjectInputStream(bis);
187 DenseOutputModel cm = (DenseOutputModel) ois.readObject();
188
189 Random random = new Random(347588535632l);
190 double maxError = 0.0;
191 for (int i = 0; i < 1000; ++i) {
192 double r = random.nextDouble();
193 double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
194 double[] interpolatedY = cm.getInterpolatedState(time).getPrimaryState();
195 double[] theoreticalY = pb.computeTheoreticalState(time);
196 double dx = interpolatedY[0] - theoreticalY[0];
197 double dy = interpolatedY[1] - theoreticalY[1];
198 double error = dx * dx + dy * dy;
199 if (error > maxError) {
200 maxError = error;
201 }
202 }
203 assertEquals(0.0, maxError, 5.5e-7);
204 } catch (ClassNotFoundException | IOException e) {
205 fail(e.getLocalizedMessage());
206 }
207
208 }
209 private boolean checkAppendError(DenseOutputModel cm,
210 double t0, double[] y0, double t1)
211 throws MathIllegalArgumentException, MathIllegalStateException {
212 try {
213 DenseOutputModel otherCm = new DenseOutputModel();
214 otherCm.handleStep(buildInterpolator(t0, y0, t1));
215 cm.append(otherCm);
216 } catch(MathIllegalArgumentException iae) {
217 return true;
218 }
219 return false;
220 }
221
222 private ODEStateInterpolator buildInterpolator(double t0, double[] y0, double t1) {
223 return new DummyStepInterpolator(t1 >= t0,
224 new ODEStateAndDerivative(t0, y0, new double[y0.length]),
225 new ODEStateAndDerivative(t1, y0, new double[y0.length]),
226 new ODEStateAndDerivative(t0, y0, new double[y0.length]),
227 new ODEStateAndDerivative(t1, y0, new double[y0.length]),
228 new EquationsMapper(null, y0.length));
229 }
230
231 public void checkValue(double value, double reference) {
232 assertTrue(FastMath.abs(value - reference) < 1.0e-10);
233 }
234
235 @BeforeEach
236 void setUp() {
237 pb = new TestProblem3(0.9);
238 double minStep = 0;
239 double maxStep = pb.getFinalTime() - pb.getInitialTime();
240 integ = new DormandPrince54Integrator(minStep, maxStep, 1.0e-8, 1.0e-8);
241 }
242
243 @AfterEach
244 void tearDown() {
245 pb = null;
246 integ = null;
247 }
248
249 }