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 org.apache.commons.math.ode.DerivativeException; 021 import org.apache.commons.math.ode.FirstOrderIntegrator; 022 import org.apache.commons.math.ode.IntegratorException; 023 import org.apache.commons.math.ode.TestProblem1; 024 import org.apache.commons.math.ode.TestProblem3; 025 import org.apache.commons.math.ode.TestProblem4; 026 import org.apache.commons.math.ode.TestProblem5; 027 import org.apache.commons.math.ode.TestProblemAbstract; 028 import org.apache.commons.math.ode.TestProblemHandler; 029 import org.apache.commons.math.ode.events.EventHandler; 030 import org.apache.commons.math.ode.nonstiff.GraggBulirschStoerIntegrator; 031 import org.apache.commons.math.ode.sampling.StepHandler; 032 import org.apache.commons.math.ode.sampling.StepInterpolator; 033 034 import junit.framework.*; 035 036 public class GraggBulirschStoerIntegratorTest 037 extends TestCase { 038 039 public GraggBulirschStoerIntegratorTest(String name) { 040 super(name); 041 } 042 043 public void testDimensionCheck() { 044 try { 045 TestProblem1 pb = new TestProblem1(); 046 AdaptiveStepsizeIntegrator integrator = 047 new GraggBulirschStoerIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10); 048 integrator.integrate(pb, 049 0.0, new double[pb.getDimension()+10], 050 1.0, new double[pb.getDimension()+10]); 051 fail("an exception should have been thrown"); 052 } catch(DerivativeException de) { 053 fail("wrong exception caught"); 054 } catch(IntegratorException ie) { 055 } 056 } 057 058 public void testNullIntervalCheck() { 059 try { 060 TestProblem1 pb = new TestProblem1(); 061 GraggBulirschStoerIntegrator integrator = 062 new GraggBulirschStoerIntegrator(0.0, 1.0, 1.0e-10, 1.0e-10); 063 integrator.integrate(pb, 064 0.0, new double[pb.getDimension()], 065 0.0, new double[pb.getDimension()]); 066 fail("an exception should have been thrown"); 067 } catch(DerivativeException de) { 068 fail("wrong exception caught"); 069 } catch(IntegratorException ie) { 070 } 071 } 072 073 public void testMinStep() { 074 075 try { 076 TestProblem5 pb = new TestProblem5(); 077 double minStep = 0.1 * Math.abs(pb.getFinalTime() - pb.getInitialTime()); 078 double maxStep = Math.abs(pb.getFinalTime() - pb.getInitialTime()); 079 double[] vecAbsoluteTolerance = { 1.0e-20, 1.0e-21 }; 080 double[] vecRelativeTolerance = { 1.0e-20, 1.0e-21 }; 081 082 FirstOrderIntegrator integ = 083 new GraggBulirschStoerIntegrator(minStep, maxStep, 084 vecAbsoluteTolerance, vecRelativeTolerance); 085 TestProblemHandler handler = new TestProblemHandler(pb, integ); 086 integ.addStepHandler(handler); 087 integ.integrate(pb, 088 pb.getInitialTime(), pb.getInitialState(), 089 pb.getFinalTime(), new double[pb.getDimension()]); 090 fail("an exception should have been thrown"); 091 } catch(DerivativeException de) { 092 fail("wrong exception caught"); 093 } catch(IntegratorException ie) { 094 } 095 096 } 097 098 public void testBackward() 099 throws DerivativeException, IntegratorException { 100 101 TestProblem5 pb = new TestProblem5(); 102 double minStep = 0; 103 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 104 double scalAbsoluteTolerance = 1.0e-8; 105 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 106 107 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(minStep, maxStep, 108 scalAbsoluteTolerance, 109 scalRelativeTolerance); 110 TestProblemHandler handler = new TestProblemHandler(pb, integ); 111 integ.addStepHandler(handler); 112 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 113 pb.getFinalTime(), new double[pb.getDimension()]); 114 115 assertTrue(handler.getLastError() < 9.0e-10); 116 assertTrue(handler.getMaximalValueError() < 9.0e-10); 117 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 118 assertEquals("Gragg-Bulirsch-Stoer", integ.getName()); 119 } 120 121 public void testIncreasingTolerance() 122 throws DerivativeException, IntegratorException { 123 124 int previousCalls = Integer.MAX_VALUE; 125 for (int i = -12; i < -4; ++i) { 126 TestProblem1 pb = new TestProblem1(); 127 double minStep = 0; 128 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 129 double absTolerance = Math.pow(10.0, i); 130 double relTolerance = absTolerance; 131 132 FirstOrderIntegrator integ = 133 new GraggBulirschStoerIntegrator(minStep, maxStep, 134 absTolerance, relTolerance); 135 TestProblemHandler handler = new TestProblemHandler(pb, integ); 136 integ.addStepHandler(handler); 137 integ.integrate(pb, 138 pb.getInitialTime(), pb.getInitialState(), 139 pb.getFinalTime(), new double[pb.getDimension()]); 140 141 // the coefficients are only valid for this test 142 // and have been obtained from trial and error 143 // there is no general relation between local and global errors 144 double ratio = handler.getMaximalValueError() / absTolerance; 145 assertTrue(ratio < 2.4); 146 assertTrue(ratio > 0.02); 147 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 148 149 int calls = pb.getCalls(); 150 assertEquals(integ.getEvaluations(), calls); 151 assertTrue(calls <= previousCalls); 152 previousCalls = calls; 153 154 } 155 156 } 157 158 public void testIntegratorControls() 159 throws DerivativeException, IntegratorException { 160 161 TestProblem3 pb = new TestProblem3(0.999); 162 GraggBulirschStoerIntegrator integ = 163 new GraggBulirschStoerIntegrator(0, pb.getFinalTime() - pb.getInitialTime(), 164 1.0e-8, 1.0e-10); 165 166 double errorWithDefaultSettings = getMaxError(integ, pb); 167 168 // stability control 169 integ.setStabilityCheck(true, 2, 1, 0.99); 170 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb)); 171 integ.setStabilityCheck(true, -1, -1, -1); 172 173 integ.setStepsizeControl(0.5, 0.99, 0.1, 2.5); 174 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb)); 175 integ.setStepsizeControl(-1, -1, -1, -1); 176 177 integ.setOrderControl(10, 0.7, 0.95); 178 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb)); 179 integ.setOrderControl(-1, -1, -1); 180 181 integ.setInterpolationControl(true, 3); 182 assertTrue(errorWithDefaultSettings < getMaxError(integ, pb)); 183 integ.setInterpolationControl(true, -1); 184 185 } 186 187 private double getMaxError(FirstOrderIntegrator integrator, TestProblemAbstract pb) 188 throws DerivativeException, IntegratorException { 189 TestProblemHandler handler = new TestProblemHandler(pb, integrator); 190 integrator.addStepHandler(handler); 191 integrator.integrate(pb, 192 pb.getInitialTime(), pb.getInitialState(), 193 pb.getFinalTime(), new double[pb.getDimension()]); 194 return handler.getMaximalValueError(); 195 } 196 197 public void testEvents() 198 throws DerivativeException, IntegratorException { 199 200 TestProblem4 pb = new TestProblem4(); 201 double minStep = 0; 202 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 203 double scalAbsoluteTolerance = 1.0e-10; 204 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance; 205 206 FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(minStep, maxStep, 207 scalAbsoluteTolerance, 208 scalRelativeTolerance); 209 TestProblemHandler handler = new TestProblemHandler(pb, integ); 210 integ.addStepHandler(handler); 211 EventHandler[] functions = pb.getEventsHandlers(); 212 for (int l = 0; l < functions.length; ++l) { 213 integ.addEventHandler(functions[l], 214 Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000); 215 } 216 assertEquals(functions.length, integ.getEventHandlers().size()); 217 integ.integrate(pb, 218 pb.getInitialTime(), pb.getInitialState(), 219 pb.getFinalTime(), new double[pb.getDimension()]); 220 221 assertTrue(handler.getMaximalValueError() < 5.0e-8); 222 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 223 assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep); 224 integ.clearEventHandlers(); 225 assertEquals(0, integ.getEventHandlers().size()); 226 227 } 228 229 public void testKepler() 230 throws DerivativeException, IntegratorException { 231 232 final TestProblem3 pb = new TestProblem3(0.9); 233 double minStep = 0; 234 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 235 double absTolerance = 1.0e-6; 236 double relTolerance = 1.0e-6; 237 238 FirstOrderIntegrator integ = 239 new GraggBulirschStoerIntegrator(minStep, maxStep, 240 absTolerance, relTolerance); 241 integ.addStepHandler(new KeplerStepHandler(pb)); 242 integ.integrate(pb, 243 pb.getInitialTime(), pb.getInitialState(), 244 pb.getFinalTime(), new double[pb.getDimension()]); 245 246 assertEquals(integ.getEvaluations(), pb.getCalls()); 247 assertTrue(pb.getCalls() < 2150); 248 249 } 250 251 public void testVariableSteps() 252 throws DerivativeException, IntegratorException { 253 254 final TestProblem3 pb = new TestProblem3(0.9); 255 double minStep = 0; 256 double maxStep = pb.getFinalTime() - pb.getInitialTime(); 257 double absTolerance = 1.0e-8; 258 double relTolerance = 1.0e-8; 259 FirstOrderIntegrator integ = 260 new GraggBulirschStoerIntegrator(minStep, maxStep, 261 absTolerance, relTolerance); 262 integ.addStepHandler(new VariableStepHandler()); 263 double stopTime = integ.integrate(pb, 264 pb.getInitialTime(), pb.getInitialState(), 265 pb.getFinalTime(), new double[pb.getDimension()]); 266 assertEquals(pb.getFinalTime(), stopTime, 1.0e-10); 267 assertEquals("Gragg-Bulirsch-Stoer", integ.getName()); 268 } 269 270 public void testUnstableDerivative() 271 throws DerivativeException, IntegratorException { 272 final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0); 273 FirstOrderIntegrator integ = 274 new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0); 275 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000); 276 double[] y = { Double.NaN }; 277 integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y); 278 assertEquals(8.0, y[0], 1.0e-12); 279 } 280 281 private static class KeplerStepHandler implements StepHandler { 282 public KeplerStepHandler(TestProblem3 pb) { 283 this.pb = pb; 284 reset(); 285 } 286 public boolean requiresDenseOutput() { 287 return true; 288 } 289 public void reset() { 290 nbSteps = 0; 291 maxError = 0; 292 } 293 public void handleStep(StepInterpolator interpolator, 294 boolean isLast) 295 throws DerivativeException { 296 297 ++nbSteps; 298 for (int a = 1; a < 100; ++a) { 299 300 double prev = interpolator.getPreviousTime(); 301 double curr = interpolator.getCurrentTime(); 302 double interp = ((100 - a) * prev + a * curr) / 100; 303 interpolator.setInterpolatedTime(interp); 304 305 double[] interpolatedY = interpolator.getInterpolatedState (); 306 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getInterpolatedTime()); 307 double dx = interpolatedY[0] - theoreticalY[0]; 308 double dy = interpolatedY[1] - theoreticalY[1]; 309 double error = dx * dx + dy * dy; 310 if (error > maxError) { 311 maxError = error; 312 } 313 } 314 if (isLast) { 315 assertTrue(maxError < 2.7e-6); 316 assertTrue(nbSteps < 80); 317 } 318 } 319 private int nbSteps; 320 private double maxError; 321 private TestProblem3 pb; 322 } 323 324 public static class VariableStepHandler implements StepHandler { 325 public VariableStepHandler() { 326 reset(); 327 } 328 public boolean requiresDenseOutput() { 329 return false; 330 } 331 public void reset() { 332 firstTime = true; 333 minStep = 0; 334 maxStep = 0; 335 } 336 public void handleStep(StepInterpolator interpolator, 337 boolean isLast) { 338 339 double step = Math.abs(interpolator.getCurrentTime() 340 - interpolator.getPreviousTime()); 341 if (firstTime) { 342 minStep = Math.abs(step); 343 maxStep = minStep; 344 firstTime = false; 345 } else { 346 if (step < minStep) { 347 minStep = step; 348 } 349 if (step > maxStep) { 350 maxStep = step; 351 } 352 } 353 354 if (isLast) { 355 assertTrue(minStep < 8.2e-3); 356 assertTrue(maxStep > 1.7); 357 } 358 } 359 private boolean firstTime; 360 private double minStep; 361 private double maxStep; 362 } 363 public static Test suite() { 364 return new TestSuite(GraggBulirschStoerIntegratorTest.class); 365 } 366 367 }