1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.FunctionEvaluationException;
22 import org.apache.commons.math.ode.DerivativeException;
23 import org.apache.commons.math.ode.FirstOrderIntegrator;
24 import org.apache.commons.math.ode.HighamHall54Integrator;
25 import org.apache.commons.math.ode.IntegratorException;
26 import org.apache.commons.math.ode.StepHandler;
27 import org.apache.commons.math.ode.StepInterpolator;
28 import org.apache.commons.math.ode.SwitchingFunction;
29
30 import junit.framework.*;
31
32 public class HighamHall54IntegratorTest
33 extends TestCase {
34
35 public HighamHall54IntegratorTest(String name) {
36 super(name);
37 }
38
39 public void testWrongDerivative() {
40 try {
41 HighamHall54Integrator integrator =
42 new HighamHall54Integrator(0.0, 1.0, 1.0e-10, 1.0e-10);
43 FirstOrderDifferentialEquations equations =
44 new FirstOrderDifferentialEquations() {
45 public void computeDerivatives(double t, double[] y, double[] dot)
46 throws DerivativeException {
47 if (t < -0.5) {
48 throw new DerivativeException("{0}", new String[] { "oops" });
49 } else {
50 throw new DerivativeException(new RuntimeException("oops"));
51 }
52 }
53 public int getDimension() {
54 return 1;
55 }
56 };
57
58 try {
59 integrator.integrate(equations, -1.0, new double[1], 0.0, new double[1]);
60 fail("an exception should have been thrown");
61 } catch(DerivativeException de) {
62
63 }
64
65 try {
66 integrator.integrate(equations, 0.0, new double[1], 1.0, new double[1]);
67 fail("an exception should have been thrown");
68 } catch(DerivativeException de) {
69
70 }
71
72 } catch (Exception e) {
73 fail("wrong exception caught: " + e.getMessage());
74 }
75 }
76
77 public void testMinStep()
78 throws DerivativeException, IntegratorException {
79
80 try {
81 TestProblem1 pb = new TestProblem1();
82 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
83 double maxStep = pb.getFinalTime() - pb.getInitialTime();
84 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
85 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
86
87 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
88 vecAbsoluteTolerance,
89 vecRelativeTolerance);
90 TestProblemHandler handler = new TestProblemHandler(pb, integ);
91 integ.setStepHandler(handler);
92 integ.integrate(pb,
93 pb.getInitialTime(), pb.getInitialState(),
94 pb.getFinalTime(), new double[pb.getDimension()]);
95 fail("an exception should have been thrown");
96 } catch(DerivativeException de) {
97 fail("wrong exception caught");
98 } catch(IntegratorException ie) {
99 }
100
101 }
102
103 public void testIncreasingTolerance()
104 throws DerivativeException, IntegratorException {
105
106 int previousCalls = Integer.MAX_VALUE;
107 for (int i = -12; i < -2; ++i) {
108 TestProblem1 pb = new TestProblem1();
109 double minStep = 0;
110 double maxStep = pb.getFinalTime() - pb.getInitialTime();
111 double scalAbsoluteTolerance = Math.pow(10.0, i);
112 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
113
114 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
115 scalAbsoluteTolerance,
116 scalRelativeTolerance);
117 TestProblemHandler handler = new TestProblemHandler(pb, integ);
118 integ.setStepHandler(handler);
119 integ.integrate(pb,
120 pb.getInitialTime(), pb.getInitialState(),
121 pb.getFinalTime(), new double[pb.getDimension()]);
122
123
124
125
126 assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance));
127 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
128
129 int calls = pb.getCalls();
130 assertTrue(calls <= previousCalls);
131 previousCalls = calls;
132
133 }
134
135 }
136
137 public void testSwitchingFunctions()
138 throws DerivativeException, IntegratorException {
139
140 TestProblem4 pb = new TestProblem4();
141 double minStep = 0;
142 double maxStep = pb.getFinalTime() - pb.getInitialTime();
143 double scalAbsoluteTolerance = 1.0e-8;
144 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
145
146 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
147 scalAbsoluteTolerance,
148 scalRelativeTolerance);
149 TestProblemHandler handler = new TestProblemHandler(pb, integ);
150 integ.setStepHandler(handler);
151 SwitchingFunction[] functions = pb.getSwitchingFunctions();
152 for (int l = 0; l < functions.length; ++l) {
153 integ.addSwitchingFunction(functions[l],
154 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
155 }
156 integ.integrate(pb,
157 pb.getInitialTime(), pb.getInitialState(),
158 pb.getFinalTime(), new double[pb.getDimension()]);
159
160 assertTrue(handler.getMaximalValueError() < 1.0e-7);
161 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
162 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
163
164 }
165
166 public void testSwitchingFunctionsError()
167 throws DerivativeException, IntegratorException {
168
169 final TestProblem1 pb = new TestProblem1();
170 double minStep = 0;
171 double maxStep = pb.getFinalTime() - pb.getInitialTime();
172 double scalAbsoluteTolerance = 1.0e-8;
173 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
174
175 FirstOrderIntegrator integ =
176 new HighamHall54Integrator(minStep, maxStep,
177 scalAbsoluteTolerance, scalRelativeTolerance);
178 TestProblemHandler handler = new TestProblemHandler(pb, integ);
179 integ.setStepHandler(handler);
180
181 integ.addSwitchingFunction(new SwitchingFunction() {
182 public int eventOccurred(double t, double[] y) {
183 return SwitchingFunction.CONTINUE;
184 }
185 public double g(double t, double[] y) throws FunctionEvaluationException {
186 double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
187 double offset = t - middle;
188 if (offset > 0) {
189 throw new FunctionEvaluationException(t);
190 }
191 return offset;
192 }
193 public void resetState(double t, double[] y) {
194 }
195 private static final long serialVersionUID = 935652725339916361L;
196 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
197
198 try {
199 integ.integrate(pb,
200 pb.getInitialTime(), pb.getInitialState(),
201 pb.getFinalTime(), new double[pb.getDimension()]);
202 fail("an exception should have been thrown");
203 } catch (IntegratorException ie) {
204
205 } catch (Exception e) {
206 fail("wrong exception type caught");
207 }
208
209 }
210
211 public void testSwitchingFunctionsNoConvergence()
212 throws DerivativeException, IntegratorException {
213
214 final TestProblem1 pb = new TestProblem1();
215 double minStep = 0;
216 double maxStep = pb.getFinalTime() - pb.getInitialTime();
217 double scalAbsoluteTolerance = 1.0e-8;
218 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
219
220 FirstOrderIntegrator integ =
221 new HighamHall54Integrator(minStep, maxStep,
222 scalAbsoluteTolerance, scalRelativeTolerance);
223 TestProblemHandler handler = new TestProblemHandler(pb, integ);
224 integ.setStepHandler(handler);
225
226 integ.addSwitchingFunction(new SwitchingFunction() {
227 public int eventOccurred(double t, double[] y) {
228 return SwitchingFunction.CONTINUE;
229 }
230 public double g(double t, double[] y) {
231 double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2;
232 double offset = t - middle;
233 return (offset > 0) ? (offset + 0.5) : (offset - 0.5);
234 }
235 public void resetState(double t, double[] y) {
236 }
237 private static final long serialVersionUID = 935652725339916361L;
238 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
239
240 try {
241 integ.integrate(pb,
242 pb.getInitialTime(), pb.getInitialState(),
243 pb.getFinalTime(), new double[pb.getDimension()]);
244 fail("an exception should have been thrown");
245 } catch (IntegratorException ie) {
246 assertTrue(ie.getCause() != null);
247 assertTrue(ie.getCause() instanceof ConvergenceException);
248 } catch (Exception e) {
249 fail("wrong exception type caught");
250 }
251
252 }
253
254 public void testSanityChecks() {
255 try {
256 final TestProblem3 pb = new TestProblem3(0.9);
257 double minStep = 0;
258 double maxStep = pb.getFinalTime() - pb.getInitialTime();
259
260 try {
261 FirstOrderIntegrator integ =
262 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
263 integ.integrate(pb, pb.getInitialTime(), new double[6],
264 pb.getFinalTime(), new double[pb.getDimension()]);
265 fail("an exception should have been thrown");
266 } catch (IntegratorException ie) {
267
268 }
269
270 try {
271 FirstOrderIntegrator integ =
272 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
273 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
274 pb.getFinalTime(), new double[6]);
275 fail("an exception should have been thrown");
276 } catch (IntegratorException ie) {
277
278 }
279
280 try {
281 FirstOrderIntegrator integ =
282 new HighamHall54Integrator(minStep, maxStep, new double[2], new double[4]);
283 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
284 pb.getFinalTime(), new double[pb.getDimension()]);
285 fail("an exception should have been thrown");
286 } catch (IntegratorException ie) {
287
288 }
289
290 try {
291 FirstOrderIntegrator integ =
292 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[2]);
293 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
294 pb.getFinalTime(), new double[pb.getDimension()]);
295 fail("an exception should have been thrown");
296 } catch (IntegratorException ie) {
297
298 }
299
300 try {
301 FirstOrderIntegrator integ =
302 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]);
303 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
304 pb.getInitialTime(), new double[pb.getDimension()]);
305 fail("an exception should have been thrown");
306 } catch (IntegratorException ie) {
307
308 }
309
310 } catch (Exception e) {
311 fail("wrong exception caught: " + e.getMessage());
312 }
313 }
314
315 public void testKepler()
316 throws DerivativeException, IntegratorException {
317
318 final TestProblem3 pb = new TestProblem3(0.9);
319 double minStep = 0;
320 double maxStep = pb.getFinalTime() - pb.getInitialTime();
321 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
322 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
323
324 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep,
325 vecAbsoluteTolerance,
326 vecRelativeTolerance);
327 integ.setStepHandler(new KeplerHandler(pb));
328 integ.integrate(pb,
329 pb.getInitialTime(), pb.getInitialState(),
330 pb.getFinalTime(), new double[pb.getDimension()]);
331 assertEquals("Higham-Hall 5(4)", integ.getName());
332 }
333
334 private static class KeplerHandler implements StepHandler {
335 public KeplerHandler(TestProblem3 pb) {
336 this.pb = pb;
337 nbSteps = 0;
338 maxError = 0;
339 }
340 public boolean requiresDenseOutput() {
341 return false;
342 }
343 public void reset() {
344 nbSteps = 0;
345 maxError = 0;
346 }
347 public void handleStep(StepInterpolator interpolator,
348 boolean isLast) {
349
350 ++nbSteps;
351 double[] interpolatedY = interpolator.getInterpolatedState ();
352 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime());
353 double dx = interpolatedY[0] - theoreticalY[0];
354 double dy = interpolatedY[1] - theoreticalY[1];
355 double error = dx * dx + dy * dy;
356 if (error > maxError) {
357 maxError = error;
358 }
359 if (isLast) {
360 assertTrue(maxError < 4e-11);
361 assertTrue(nbSteps < 670);
362 }
363 }
364 private TestProblem3 pb;
365 private int nbSteps;
366 private double maxError;
367 }
368
369 public static Test suite() {
370 return new TestSuite(HighamHall54IntegratorTest.class);
371 }
372
373 }