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.GillIntegrator; 034 import org.apache.commons.math.ode.sampling.StepHandler; 035 import org.apache.commons.math.ode.sampling.StepInterpolator; 036 037 public class GillIntegratorTest 038 extends TestCase { 039 040 public GillIntegratorTest(String name) { 041 super(name); 042 } 043 044 public void testDimensionCheck() { 045 try { 046 TestProblem1 pb = new TestProblem1(); 047 new GillIntegrator(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 = 5; 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 GillIntegrator(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 > 5) { 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 GillIntegrator(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("Gill", 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 GillIntegrator(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 GillIntegrator(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("Gill", 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 GillIntegrator(step); 159 integ.addStepHandler(new KeplerStepHandler(pb)); 160 integ.integrate(pb, 161 pb.getInitialTime(), pb.getInitialState(), 162 pb.getFinalTime(), new double[pb.getDimension()]); 163 } 164 165 public void testUnstableDerivative() 166 throws DerivativeException, IntegratorException { 167 final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0); 168 FirstOrderIntegrator integ = new GillIntegrator(0.3); 169 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000); 170 double[] y = { Double.NaN }; 171 integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y); 172 assertEquals(8.0, y[0], 1.0e-12); 173 } 174 175 private static class KeplerStepHandler implements StepHandler { 176 public KeplerStepHandler(TestProblem3 pb) { 177 this.pb = pb; 178 reset(); 179 } 180 public boolean requiresDenseOutput() { 181 return false; 182 } 183 public void reset() { 184 maxError = 0; 185 } 186 public void handleStep(StepInterpolator interpolator, 187 boolean isLast) throws DerivativeException { 188 189 double[] interpolatedY = interpolator.getInterpolatedState(); 190 double[] theoreticalY = pb.computeTheoreticalState(interpolator.getCurrentTime()); 191 double dx = interpolatedY[0] - theoreticalY[0]; 192 double dy = interpolatedY[1] - theoreticalY[1]; 193 double error = dx * dx + dy * dy; 194 if (error > maxError) { 195 maxError = error; 196 } 197 if (isLast) { 198 // even with more than 1000 evaluations per period, 199 // RK4 is not able to integrate such an eccentric 200 // orbit with a good accuracy 201 assertTrue(maxError > 0.001); 202 } 203 } 204 private double maxError; 205 private TestProblem3 pb; 206 } 207 208 public void testStepSize() 209 throws DerivativeException, IntegratorException { 210 final double step = 1.23456; 211 FirstOrderIntegrator integ = new GillIntegrator(step); 212 integ.addStepHandler(new StepHandler() { 213 public void handleStep(StepInterpolator interpolator, boolean isLast) { 214 if (! isLast) { 215 assertEquals(step, 216 interpolator.getCurrentTime() - interpolator.getPreviousTime(), 217 1.0e-12); 218 } 219 } 220 public boolean requiresDenseOutput() { 221 return false; 222 } 223 public void reset() { 224 } 225 }); 226 integ.integrate(new FirstOrderDifferentialEquations() { 227 private static final long serialVersionUID = 0L; 228 public void computeDerivatives(double t, double[] y, double[] dot) { 229 dot[0] = 1.0; 230 } 231 public int getDimension() { 232 return 1; 233 } 234 }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); 235 } 236 237 public static Test suite() { 238 return new TestSuite(GillIntegratorTest.class); 239 } 240 241 }