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