1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.nonstiff;
19
20 import junit.framework.*;
21
22 import org.apache.commons.math.ode.DerivativeException;
23 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
24 import org.apache.commons.math.ode.FirstOrderIntegrator;
25 import org.apache.commons.math.ode.IntegratorException;
26 import org.apache.commons.math.ode.TestProblem1;
27 import org.apache.commons.math.ode.TestProblem3;
28 import org.apache.commons.math.ode.TestProblem5;
29 import org.apache.commons.math.ode.TestProblemAbstract;
30 import org.apache.commons.math.ode.TestProblemFactory;
31 import org.apache.commons.math.ode.TestProblemHandler;
32 import org.apache.commons.math.ode.events.EventHandler;
33 import org.apache.commons.math.ode.nonstiff.ClassicalRungeKuttaIntegrator;
34 import org.apache.commons.math.ode.sampling.StepHandler;
35 import org.apache.commons.math.ode.sampling.StepInterpolator;
36
37 public class ClassicalRungeKuttaIntegratorTest
38 extends TestCase {
39
40 public ClassicalRungeKuttaIntegratorTest(String name) {
41 super(name);
42 }
43
44 public void testSanityChecks() {
45 try {
46 TestProblem1 pb = new TestProblem1();
47 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
48 0.0, new double[pb.getDimension()+10],
49 1.0, new double[pb.getDimension()]);
50 fail("an exception should have been thrown");
51 } catch(DerivativeException de) {
52 fail("wrong exception caught");
53 } catch(IntegratorException ie) {
54 }
55 try {
56 TestProblem1 pb = new TestProblem1();
57 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
58 0.0, new double[pb.getDimension()],
59 1.0, new double[pb.getDimension()+10]);
60 fail("an exception should have been thrown");
61 } catch(DerivativeException de) {
62 fail("wrong exception caught");
63 } catch(IntegratorException ie) {
64 }
65 try {
66 TestProblem1 pb = new TestProblem1();
67 new ClassicalRungeKuttaIntegrator(0.01).integrate(pb,
68 0.0, new double[pb.getDimension()],
69 0.0, new double[pb.getDimension()]);
70 fail("an exception should have been thrown");
71 } catch(DerivativeException de) {
72 fail("wrong exception caught");
73 } catch(IntegratorException ie) {
74 }
75 }
76
77 public void testDecreasingSteps()
78 throws DerivativeException, IntegratorException {
79
80 TestProblemAbstract[] problems = TestProblemFactory.getProblems();
81 for (int k = 0; k < problems.length; ++k) {
82
83 double previousError = Double.NaN;
84 for (int i = 4; i < 10; ++i) {
85
86 TestProblemAbstract pb = problems[k].copy();
87 double step = (pb.getFinalTime() - pb.getInitialTime()) * Math.pow(2.0, -i);
88
89 FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
90 TestProblemHandler handler = new TestProblemHandler(pb, integ);
91 integ.addStepHandler(handler);
92 EventHandler[] functions = pb.getEventsHandlers();
93 for (int l = 0; l < functions.length; ++l) {
94 integ.addEventHandler(functions[l],
95 Double.POSITIVE_INFINITY, 1.0e-6 * step, 1000);
96 }
97 assertEquals(functions.length, integ.getEventHandlers().size());
98 double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
99 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
209
210
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 }