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;
019
020 import junit.framework.*;
021 import java.util.Random;
022
023 import org.apache.commons.math.ode.ContinuousOutputModel;
024 import org.apache.commons.math.ode.DerivativeException;
025 import org.apache.commons.math.ode.FirstOrderIntegrator;
026 import org.apache.commons.math.ode.IntegratorException;
027 import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator;
028 import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
029 import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
030 import org.apache.commons.math.ode.sampling.StepInterpolator;
031
032 public class ContinuousOutputModelTest
033 extends TestCase {
034
035 public ContinuousOutputModelTest(String name) {
036 super(name);
037 pb = null;
038 integ = null;
039 }
040
041 public void testBoundaries()
042 throws DerivativeException, IntegratorException {
043 integ.addStepHandler(new ContinuousOutputModel());
044 integ.integrate(pb,
045 pb.getInitialTime(), pb.getInitialState(),
046 pb.getFinalTime(), new double[pb.getDimension()]);
047 ContinuousOutputModel cm = (ContinuousOutputModel) integ.getStepHandlers().iterator().next();
048 cm.setInterpolatedTime(2.0 * pb.getInitialTime() - pb.getFinalTime());
049 cm.setInterpolatedTime(2.0 * pb.getFinalTime() - pb.getInitialTime());
050 cm.setInterpolatedTime(0.5 * (pb.getFinalTime() + pb.getInitialTime()));
051 }
052
053 public void testRandomAccess()
054 throws DerivativeException, IntegratorException {
055
056 ContinuousOutputModel cm = new ContinuousOutputModel();
057 integ.addStepHandler(cm);
058 integ.integrate(pb,
059 pb.getInitialTime(), pb.getInitialState(),
060 pb.getFinalTime(), new double[pb.getDimension()]);
061
062 Random random = new Random(347588535632l);
063 double maxError = 0.0;
064 for (int i = 0; i < 1000; ++i) {
065 double r = random.nextDouble();
066 double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
067 cm.setInterpolatedTime(time);
068 double[] interpolatedY = cm.getInterpolatedState ();
069 double[] theoreticalY = pb.computeTheoreticalState(time);
070 double dx = interpolatedY[0] - theoreticalY[0];
071 double dy = interpolatedY[1] - theoreticalY[1];
072 double error = dx * dx + dy * dy;
073 if (error > maxError) {
074 maxError = error;
075 }
076 }
077
078 assertTrue(maxError < 1.0e-9);
079
080 }
081
082 public void testModelsMerging()
083 throws DerivativeException, IntegratorException {
084
085 // theoretical solution: y[0] = cos(t), y[1] = sin(t)
086 FirstOrderDifferentialEquations problem =
087 new FirstOrderDifferentialEquations() {
088 private static final long serialVersionUID = 2472449657345878299L;
089 public void computeDerivatives(double t, double[] y, double[] dot)
090 throws DerivativeException {
091 dot[0] = -y[1];
092 dot[1] = y[0];
093 }
094 public int getDimension() {
095 return 2;
096 }
097 };
098
099 // integrate backward from π to 0;
100 ContinuousOutputModel cm1 = new ContinuousOutputModel();
101 FirstOrderIntegrator integ1 =
102 new DormandPrince853Integrator(0, 1.0, 1.0e-8, 1.0e-8);
103 integ1.addStepHandler(cm1);
104 integ1.integrate(problem, Math.PI, new double[] { -1.0, 0.0 },
105 0, new double[2]);
106
107 // integrate backward from 2π to π
108 ContinuousOutputModel cm2 = new ContinuousOutputModel();
109 FirstOrderIntegrator integ2 =
110 new DormandPrince853Integrator(0, 0.1, 1.0e-12, 1.0e-12);
111 integ2.addStepHandler(cm2);
112 integ2.integrate(problem, 2.0 * Math.PI, new double[] { 1.0, 0.0 },
113 Math.PI, new double[2]);
114
115 // merge the two half circles
116 ContinuousOutputModel cm = new ContinuousOutputModel();
117 cm.append(cm2);
118 cm.append(new ContinuousOutputModel());
119 cm.append(cm1);
120
121 // check circle
122 assertEquals(2.0 * Math.PI, cm.getInitialTime(), 1.0e-12);
123 assertEquals(0, cm.getFinalTime(), 1.0e-12);
124 assertEquals(cm.getFinalTime(), cm.getInterpolatedTime(), 1.0e-12);
125 for (double t = 0; t < 2.0 * Math.PI; t += 0.1) {
126 cm.setInterpolatedTime(t);
127 double[] y = cm.getInterpolatedState();
128 assertEquals(Math.cos(t), y[0], 1.0e-7);
129 assertEquals(Math.sin(t), y[1], 1.0e-7);
130 }
131
132 }
133
134 public void testErrorConditions()
135 throws DerivativeException {
136
137 ContinuousOutputModel cm = new ContinuousOutputModel();
138 cm.handleStep(buildInterpolator(0, new double[] { 0.0, 1.0, -2.0 }, 1), true);
139
140 // dimension mismatch
141 assertTrue(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0 }, 2.0));
142
143 // hole between time ranges
144 assertTrue(checkAppendError(cm, 10.0, new double[] { 0.0, 1.0, -2.0 }, 20.0));
145
146 // propagation direction mismatch
147 assertTrue(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0, -2.0 }, 0.0));
148
149 // no errors
150 assertFalse(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0, -2.0 }, 2.0));
151
152 }
153
154 private boolean checkAppendError(ContinuousOutputModel cm,
155 double t0, double[] y0, double t1)
156 throws DerivativeException {
157 try {
158 ContinuousOutputModel otherCm = new ContinuousOutputModel();
159 otherCm.handleStep(buildInterpolator(t0, y0, t1), true);
160 cm.append(otherCm);
161 } catch(IllegalArgumentException iae) {
162 //expected behavior
163 return true;
164 }
165 return false;
166 }
167
168 private StepInterpolator buildInterpolator(double t0, double[] y0, double t1) {
169 DummyStepInterpolator interpolator = new DummyStepInterpolator(y0, t1 >= t0);
170 interpolator.storeTime(t0);
171 interpolator.shift();
172 interpolator.storeTime(t1);
173 return interpolator;
174 }
175
176 public void checkValue(double value, double reference) {
177 assertTrue(Math.abs(value - reference) < 1.0e-10);
178 }
179
180 public static Test suite() {
181 return new TestSuite(ContinuousOutputModelTest.class);
182 }
183
184 @Override
185 public void setUp() {
186 pb = new TestProblem3(0.9);
187 double minStep = 0;
188 double maxStep = pb.getFinalTime() - pb.getInitialTime();
189 integ = new DormandPrince54Integrator(minStep, maxStep, 1.0e-8, 1.0e-8);
190 }
191
192 @Override
193 public void tearDown() {
194 pb = null;
195 integ = null;
196 }
197
198 TestProblem3 pb;
199 FirstOrderIntegrator integ;
200
201 }