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.*; 021 022 import org.apache.commons.math.ode.DerivativeException; 023 import org.apache.commons.math.ode.FirstOrderDifferentialEquations; 024 import org.apache.commons.math.ode.FirstOrderIntegrator; 025 import org.apache.commons.math.ode.IntegratorException; 026 import org.apache.commons.math.ode.TestProblem1; 027 import org.apache.commons.math.ode.TestProblem3; 028 import org.apache.commons.math.ode.TestProblem5; 029 import org.apache.commons.math.ode.TestProblemAbstract; 030 import org.apache.commons.math.ode.TestProblemFactory; 031 import org.apache.commons.math.ode.TestProblemHandler; 032 import org.apache.commons.math.ode.events.EventHandler; 033 import org.apache.commons.math.ode.nonstiff.ClassicalRungeKuttaIntegrator; 034 import org.apache.commons.math.ode.sampling.StepHandler; 035 import org.apache.commons.math.ode.sampling.StepInterpolator; 036 037 public class ClassicalRungeKuttaIntegratorTest 038 extends TestCase { 039 040 public ClassicalRungeKuttaIntegratorTest(String name) { 041 super(name); 042 } 043 044 public void testSanityChecks() { 045 try { 046 TestProblem1 pb = new TestProblem1(); 047 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb, 048 0.0, new double[pb.getDimension()+10], 049 1.0, new double[pb.getDimension()]); 050 fail("an exception should have been thrown"); 051 } catch(DerivativeException de) { 052 fail("wrong exception caught"); 053 } catch(IntegratorException ie) { 054 } 055 try { 056 TestProblem1 pb = new TestProblem1(); 057 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb, 058 0.0, new double[pb.getDimension()], 059 1.0, new double[pb.getDimension()+10]); 060 fail("an exception should have been thrown"); 061 } catch(DerivativeException de) { 062 fail("wrong exception caught"); 063 } catch(IntegratorException ie) { 064 } 065 try { 066 TestProblem1 pb = new TestProblem1(); 067 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb, 068 0.0, new double[pb.getDimension()], 069 0.0, new double[pb.getDimension()]); 070 fail("an exception should have been thrown"); 071 } catch(DerivativeException de) { 072 fail("wrong exception caught"); 073 } catch(IntegratorException ie) { 074 } 075 } 076 077 public void testDecreasingSteps() 078 throws DerivativeException, IntegratorException { 079 080 TestProblemAbstract[] problems = TestProblemFactory.getProblems(); 081 for (int k = 0; k < problems.length; ++k) { 082 083 double previousError = Double.NaN; 084 for (int i = 4; i < 10; ++i) { 085 086 TestProblemAbstract pb = problems[k].copy(); 087 double step = (pb.getFinalTime() - pb.getInitialTime()) * Math.pow(2.0, -i); 088 089 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 090 TestProblemHandler handler = new TestProblemHandler(pb, integ); 091 integ.addStepHandler(handler); 092 EventHandler[] functions = pb.getEventsHandlers(); 093 for (int l = 0; l < functions.length; ++l) { 094 integ.addEventHandler(functions[l], 095 Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000); 096 } 097 assertEquals(functions.length, integ.getEventHandlers().size()); 098 double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 099 pb.getFinalTime(), new double[pb.getDimension()]); 100 if (functions.length == 0) { 101 assertEquals(pb.getFinalTime(), stopTime, 1.0e-10); 102 } 103 104 double error = handler.getMaximalValueError(); 105 if (i > 4) { 106 assertTrue(error < Math.abs(previousError)); 107 } 108 previousError = error; 109 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 110 integ.clearEventHandlers(); 111 assertEquals(0, integ.getEventHandlers().size()); 112 } 113 114 } 115 116 } 117 118 public void testSmallStep() 119 throws DerivativeException, IntegratorException { 120 121 TestProblem1 pb = new TestProblem1(); 122 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001; 123 124 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 125 TestProblemHandler handler = new TestProblemHandler(pb, integ); 126 integ.addStepHandler(handler); 127 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 128 pb.getFinalTime(), new double[pb.getDimension()]); 129 130 assertTrue(handler.getLastError() < 2.0e-13); 131 assertTrue(handler.getMaximalValueError() < 4.0e-12); 132 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 133 assertEquals("classical Runge-Kutta", integ.getName()); 134 } 135 136 public void testBigStep() 137 throws DerivativeException, IntegratorException { 138 139 TestProblem1 pb = new TestProblem1(); 140 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2; 141 142 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 143 TestProblemHandler handler = new TestProblemHandler(pb, integ); 144 integ.addStepHandler(handler); 145 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 146 pb.getFinalTime(), new double[pb.getDimension()]); 147 148 assertTrue(handler.getLastError() > 0.0004); 149 assertTrue(handler.getMaximalValueError() > 0.005); 150 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 151 152 } 153 154 public void testBackward() 155 throws DerivativeException, IntegratorException { 156 157 TestProblem5 pb = new TestProblem5(); 158 double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001; 159 160 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 161 TestProblemHandler handler = new TestProblemHandler(pb, integ); 162 integ.addStepHandler(handler); 163 integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(), 164 pb.getFinalTime(), new double[pb.getDimension()]); 165 166 assertTrue(handler.getLastError() < 5.0e-10); 167 assertTrue(handler.getMaximalValueError() < 7.0e-10); 168 assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); 169 assertEquals("classical Runge-Kutta", integ.getName()); 170 } 171 172 public void testKepler() 173 throws DerivativeException, IntegratorException { 174 175 final TestProblem3 pb = new TestProblem3(0.9); 176 double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.0003; 177 178 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 179 integ.addStepHandler(new KeplerHandler(pb)); 180 integ.integrate(pb, 181 pb.getInitialTime(), pb.getInitialState(), 182 pb.getFinalTime(), new double[pb.getDimension()]); 183 } 184 185 private static class KeplerHandler implements StepHandler { 186 public KeplerHandler(TestProblem3 pb) { 187 this.pb = pb; 188 reset(); 189 } 190 public boolean requiresDenseOutput() { 191 return false; 192 } 193 public void reset() { 194 maxError = 0; 195 } 196 public void handleStep(StepInterpolator interpolator, 197 boolean isLast) throws DerivativeException { 198 199 double[] interpolatedY = interpolator.getInterpolatedState (); 200 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime()); 201 double dx = interpolatedY[0] - theoreticalY[0]; 202 double dy = interpolatedY[1] - theoreticalY[1]; 203 double error = dx * dx + dy * dy; 204 if (error > maxError) { 205 maxError = error; 206 } 207 if (isLast) { 208 // even with more than 1000 evaluations per period, 209 // RK4 is not able to integrate such an eccentric 210 // orbit with a good accuracy 211 assertTrue(maxError > 0.005); 212 } 213 } 214 private double maxError = 0; 215 private TestProblem3 pb; 216 } 217 218 public void testStepSize() 219 throws DerivativeException, IntegratorException { 220 final double step = 1.23456; 221 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); 222 integ.addStepHandler(new StepHandler() { 223 public void handleStep(StepInterpolator interpolator, boolean isLast) { 224 if (! isLast) { 225 assertEquals(step, 226 interpolator.getCurrentTime() - interpolator.getPreviousTime(), 227 1.0e-12); 228 } 229 } 230 public boolean requiresDenseOutput() { 231 return false; 232 } 233 public void reset() { 234 } 235 }); 236 integ.integrate(new FirstOrderDifferentialEquations() { 237 private static final long serialVersionUID = 0L; 238 public void computeDerivatives(double t, double[] y, double[] dot) { 239 dot[0] = 1.0; 240 } 241 public int getDimension() { 242 return 1; 243 } 244 }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); 245 } 246 247 public static Test suite() { 248 return new TestSuite(ClassicalRungeKuttaIntegratorTest.class); 249 } 250 251 }