001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.ode.nonstiff; 019 020 import junit.framework.Test; 021 import junit.framework.TestCase; 022 import junit.framework.TestSuite; 023 024 import org.apache.commons.math.ConvergenceException; 025 import org.apache.commons.math.ode.DerivativeException; 026 import org.apache.commons.math.ode.FirstOrderDifferentialEquations; 027 import org.apache.commons.math.ode.FirstOrderIntegrator; 028 import org.apache.commons.math.ode.IntegratorException; 029 import org.apache.commons.math.ode.TestProblem1; 030 import org.apache.commons.math.ode.TestProblem3; 031 import org.apache.commons.math.ode.TestProblem4; 032 import org.apache.commons.math.ode.TestProblem5; 033 import org.apache.commons.math.ode.TestProblemHandler; 034 import org.apache.commons.math.ode.events.EventException; 035 import org.apache.commons.math.ode.events.EventHandler; 036 import org.apache.commons.math.ode.nonstiff.HighamHall54Integrator; 037 import org.apache.commons.math.ode.sampling.StepHandler; 038 import org.apache.commons.math.ode.sampling.StepInterpolator; 039 040 public class HighamHall54IntegratorTest 041 extends TestCase { 042 043 public HighamHall54IntegratorTest(String name) { 044 super(name); 045 } 046 047 public void testWrongDerivative() { 048 try { 049 HighamHall54Integrator integrator = 050 new HighamHall54Integrator(0.0, 1.0, 1.0e-10, 1.0e-10); 051 FirstOrderDifferentialEquations equations = 052 new FirstOrderDifferentialEquations() { 053 private static final long serialVersionUID = -1157081786301178032L; 054 public void computeDerivatives(double t, double[] y, double[] dot) 055 throws DerivativeException { 056 if (t < -0.5) { 057 throw new DerivativeException("{0}", "oops"); 058 } else { 059 throw new DerivativeException(new RuntimeException("oops")); 060 } 061 } 062 public int getDimension() { 063 return 1; 064 } 065 }; 066 067 try { 068 integrator.integrate(equations, -1.0, new double[1], 0.0, new double[1]); 069 fail("an exception should have been thrown"); 070 } catch(DerivativeException de) { 071 // expected behavior 072 } 073 074 try { 075 integrator.integrate(equations, 0.0, new double[1], 1.0, new double[1]); 076 fail("an exception should have been thrown"); 077 } catch(DerivativeException de) { 078 // expected behavior 079 } 080 081 } catch (Exception e) { 082 fail("wrong exception caught: " + e.getMessage()); 083 } 084 } 085 086 public void testMinStep() { 087 088 try { 089 TestProblem1 pb = new TestProblem1(); 090 double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime()); 091 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 092 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 }; 093 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 }; 094 095 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep, 096 vecAbsoluteTolerance, 097 vecRelativeTolerance); 098 TestProblemHandler handler = new TestProblemHandler(pb, integ); 099 integ.addStepHandler(handler); 100 integ.integrate(pb, 101 pb.getInitialTime(), pb.getInitialState(), 102 pb.getFinalTime(), new double[pb.getDimension()]); 103 fail("an exception should have been thrown"); 104 } catch(DerivativeException de) { 105 fail("wrong exception caught"); 106 } catch(IntegratorException ie) { 107 } 108 109 } 110 111 public void testIncreasingTolerance() 112 throws DerivativeException, IntegratorException { 113 114 int previousCalls = Integer.MAX_VALUE; 115 for (int i = -12; i < -2; ++i) { 116 TestProblem1 pb = new TestProblem1(); 117 double minStep = 0; 118 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 119 double scalAbsoluteTolerance = Math.pow(10.0, i); 120 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 121 122 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep, 123 scalAbsoluteTolerance, 124 scalRelativeTolerance); 125 TestProblemHandler handler = new TestProblemHandler(pb, integ); 126 integ.addStepHandler(handler); 127 integ.integrate(pb, 128 pb.getInitialTime(), pb.getInitialState(), 129 pb.getFinalTime(), new double[pb.getDimension()]); 130 131 // the 1.3 factor is only valid for this test 132 // and has been obtained from trial and error 133 // there is no general relation between local and global errors 134 assertTrue(handler.getMaximalValueError() < (1.3 * scalAbsoluteTolerance)); 135 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 136 137 int calls = pb.getCalls(); 138 assertEquals(integ.getEvaluations(), calls); 139 assertTrue(calls <= previousCalls); 140 previousCalls = calls; 141 142 } 143 144 } 145 146 public void testBackward() 147 throws DerivativeException, IntegratorException { 148 149 TestProblem5 pb = new TestProblem5(); 150 double minStep = 0; 151 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 152 double scalAbsoluteTolerance = 1.0e-8; 153 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 154 155 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep, 156 scalAbsoluteTolerance, 157 scalRelativeTolerance); 158 TestProblemHandler handler = new TestProblemHandler(pb, integ); 159 integ.addStepHandler(handler); 160 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 161 pb.getFinalTime(), new double[pb.getDimension()]); 162 163 assertTrue(handler.getLastError() < 5.0e-7); 164 assertTrue(handler.getMaximalValueError() < 5.0e-7); 165 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 166 assertEquals("Higham-Hall 5(4)", integ.getName()); 167 } 168 169 public void testEvents() 170 throws DerivativeException, IntegratorException { 171 172 TestProblem4 pb = new TestProblem4(); 173 double minStep = 0; 174 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 175 double scalAbsoluteTolerance = 1.0e-8; 176 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 177 178 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep, 179 scalAbsoluteTolerance, 180 scalRelativeTolerance); 181 TestProblemHandler handler = new TestProblemHandler(pb, integ); 182 integ.addStepHandler(handler); 183 EventHandler[] functions = pb.getEventsHandlers(); 184 for (int l = 0; l < functions.length; ++l) { 185 integ.addEventHandler(functions[l], 186 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000); 187 } 188 assertEquals(functions.length, integ.getEventHandlers().size()); 189 integ.integrate(pb, 190 pb.getInitialTime(), pb.getInitialState(), 191 pb.getFinalTime(), new double[pb.getDimension()]); 192 193 assertTrue(handler.getMaximalValueError() < 1.0e-7); 194 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 195 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep); 196 integ.clearEventHandlers(); 197 assertEquals(0, integ.getEventHandlers().size()); 198 199 } 200 201 public void testEventsErrors() { 202 203 final TestProblem1 pb = new TestProblem1(); 204 double minStep = 0; 205 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 206 double scalAbsoluteTolerance = 1.0e-8; 207 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 208 209 FirstOrderIntegrator integ = 210 new HighamHall54Integrator(minStep, maxStep, 211 scalAbsoluteTolerance, scalRelativeTolerance); 212 TestProblemHandler handler = new TestProblemHandler(pb, integ); 213 integ.addStepHandler(handler); 214 215 integ.addEventHandler(new EventHandler() { 216 public int eventOccurred(double t, double[] y, boolean increasing) { 217 return EventHandler.CONTINUE; 218 } 219 public double g(double t, double[] y) throws EventException { 220 double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2; 221 double offset = t - middle; 222 if (offset > 0) { 223 throw new EventException("Evaluation failed for argument = {0}", t); 224 } 225 return offset; 226 } 227 public void resetState(double t, double[] y) { 228 } 229 private static final long serialVersionUID = 935652725339916361L; 230 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000); 231 232 try { 233 integ.integrate(pb, 234 pb.getInitialTime(), pb.getInitialState(), 235 pb.getFinalTime(), new double[pb.getDimension()]); 236 fail("an exception should have been thrown"); 237 } catch (IntegratorException ie) { 238 // expected behavior 239 } catch (Exception e) { 240 fail("wrong exception type caught"); 241 } 242 243 } 244 245 public void testEventsNoConvergence() { 246 247 final TestProblem1 pb = new TestProblem1(); 248 double minStep = 0; 249 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 250 double scalAbsoluteTolerance = 1.0e-8; 251 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 252 253 FirstOrderIntegrator integ = 254 new HighamHall54Integrator(minStep, maxStep, 255 scalAbsoluteTolerance, scalRelativeTolerance); 256 TestProblemHandler handler = new TestProblemHandler(pb, integ); 257 integ.addStepHandler(handler); 258 259 integ.addEventHandler(new EventHandler() { 260 public int eventOccurred(double t, double[] y, boolean increasing) { 261 return EventHandler.CONTINUE; 262 } 263 public double g(double t, double[] y) { 264 double middle = (pb.getInitialTime() + pb.getFinalTime()) / 2; 265 double offset = t - middle; 266 return (offset > 0) ? (offset + 0.5) : (offset - 0.5); 267 } 268 public void resetState(double t, double[] y) { 269 } 270 private static final long serialVersionUID = 935652725339916361L; 271 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3); 272 273 try { 274 integ.integrate(pb, 275 pb.getInitialTime(), pb.getInitialState(), 276 pb.getFinalTime(), new double[pb.getDimension()]); 277 fail("an exception should have been thrown"); 278 } catch (IntegratorException ie) { 279 assertTrue(ie.getCause() != null); 280 assertTrue(ie.getCause() instanceof ConvergenceException); 281 } catch (Exception e) { 282 fail("wrong exception type caught"); 283 } 284 285 } 286 287 public void testSanityChecks() { 288 try { 289 final TestProblem3 pb = new TestProblem3(0.9); 290 double minStep = 0; 291 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 292 293 try { 294 FirstOrderIntegrator integ = 295 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]); 296 integ.integrate(pb, pb.getInitialTime(), new double[6], 297 pb.getFinalTime(), new double[pb.getDimension()]); 298 fail("an exception should have been thrown"); 299 } catch (IntegratorException ie) { 300 // expected behavior 301 } 302 303 try { 304 FirstOrderIntegrator integ = 305 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]); 306 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 307 pb.getFinalTime(), new double[6]); 308 fail("an exception should have been thrown"); 309 } catch (IntegratorException ie) { 310 // expected behavior 311 } 312 313 try { 314 FirstOrderIntegrator integ = 315 new HighamHall54Integrator(minStep, maxStep, new double[2], new double[4]); 316 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 317 pb.getFinalTime(), new double[pb.getDimension()]); 318 fail("an exception should have been thrown"); 319 } catch (IntegratorException ie) { 320 // expected behavior 321 } 322 323 try { 324 FirstOrderIntegrator integ = 325 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[2]); 326 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 327 pb.getFinalTime(), new double[pb.getDimension()]); 328 fail("an exception should have been thrown"); 329 } catch (IntegratorException ie) { 330 // expected behavior 331 } 332 333 try { 334 FirstOrderIntegrator integ = 335 new HighamHall54Integrator(minStep, maxStep, new double[4], new double[4]); 336 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 337 pb.getInitialTime(), new double[pb.getDimension()]); 338 fail("an exception should have been thrown"); 339 } catch (IntegratorException ie) { 340 // expected behavior 341 } 342 343 } catch (Exception e) { 344 fail("wrong exception caught: " + e.getMessage()); 345 } 346 } 347 348 public void testKepler() 349 throws DerivativeException, IntegratorException { 350 351 final TestProblem3 pb = new TestProblem3(0.9); 352 double minStep = 0; 353 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 354 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 }; 355 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 }; 356 357 FirstOrderIntegrator integ = new HighamHall54Integrator(minStep, maxStep, 358 vecAbsoluteTolerance, 359 vecRelativeTolerance); 360 integ.addStepHandler(new KeplerHandler(pb)); 361 integ.integrate(pb, 362 pb.getInitialTime(), pb.getInitialState(), 363 pb.getFinalTime(), new double[pb.getDimension()]); 364 assertEquals("Higham-Hall 5(4)", integ.getName()); 365 } 366 367 private static class KeplerHandler implements StepHandler { 368 public KeplerHandler(TestProblem3 pb) { 369 this.pb = pb; 370 nbSteps = 0; 371 maxError = 0; 372 } 373 public boolean requiresDenseOutput() { 374 return false; 375 } 376 public void reset() { 377 nbSteps = 0; 378 maxError = 0; 379 } 380 public void handleStep(StepInterpolator interpolator, 381 boolean isLast) throws DerivativeException { 382 383 ++nbSteps; 384 double[] interpolatedY = interpolator.getInterpolatedState(); 385 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime()); 386 double dx = interpolatedY[0] - theoreticalY[0]; 387 double dy = interpolatedY[1] - theoreticalY[1]; 388 double error = dx * dx + dy * dy; 389 if (error > maxError) { 390 maxError = error; 391 } 392 if (isLast) { 393 assertTrue(maxError < 4e-11); 394 assertTrue(nbSteps < 670); 395 } 396 } 397 private TestProblem3 pb; 398 private int nbSteps; 399 private double maxError; 400 } 401 402 public static Test suite() { 403 return new TestSuite(HighamHall54IntegratorTest.class); 404 } 405 406 }