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