View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math4.legacy.ode;
19  
20  import java.util.Random;
21  
22  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
23  import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
24  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25  import org.apache.commons.math4.legacy.exception.NoBracketingException;
26  import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
27  import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince54Integrator;
28  import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853Integrator;
29  import org.apache.commons.math4.legacy.ode.sampling.DummyStepInterpolator;
30  import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
31  import org.apache.commons.math4.core.jdkmath.JdkMath;
32  import org.junit.After;
33  import org.junit.Assert;
34  import org.junit.Before;
35  import org.junit.Test;
36  
37  public class ContinuousOutputModelTest {
38  
39    public ContinuousOutputModelTest() {
40      pb    = null;
41      integ = null;
42    }
43  
44    @Test
45    public void testBoundaries() throws DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException {
46      integ.addStepHandler(new ContinuousOutputModel());
47      integ.integrate(pb,
48                      pb.getInitialTime(), pb.getInitialState(),
49                      pb.getFinalTime(), new double[pb.getDimension()]);
50      ContinuousOutputModel cm = (ContinuousOutputModel) integ.getStepHandlers().iterator().next();
51      cm.setInterpolatedTime(2.0 * pb.getInitialTime() - pb.getFinalTime());
52      cm.setInterpolatedTime(2.0 * pb.getFinalTime() - pb.getInitialTime());
53      cm.setInterpolatedTime(0.5 * (pb.getFinalTime() + pb.getInitialTime()));
54    }
55  
56    @Test
57    public void testRandomAccess() throws DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException {
58  
59      ContinuousOutputModel cm = new ContinuousOutputModel();
60      integ.addStepHandler(cm);
61      integ.integrate(pb,
62                      pb.getInitialTime(), pb.getInitialState(),
63                      pb.getFinalTime(), new double[pb.getDimension()]);
64  
65      Random random = new Random(347588535632L);
66      double maxError    = 0.0;
67      double maxErrorDot = 0.0;
68      for (int i = 0; i < 1000; ++i) {
69        double r = random.nextDouble();
70        double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
71        cm.setInterpolatedTime(time);
72        double[] interpolatedY    = cm.getInterpolatedState();
73        double[] interpolatedYDot = cm.getInterpolatedDerivatives();
74        double[] theoreticalY     = pb.computeTheoreticalState(time);
75        double[] theoreticalYDot  = new double[pb.getDimension()];
76        pb.doComputeDerivatives(time, theoreticalY, theoreticalYDot);
77        double dx = interpolatedY[0] - theoreticalY[0];
78        double dy = interpolatedY[1] - theoreticalY[1];
79        double error = dx * dx + dy * dy;
80        maxError = JdkMath.max(maxError, error);
81        double dxDot = interpolatedYDot[0] - theoreticalYDot[0];
82        double dyDot = interpolatedYDot[1] - theoreticalYDot[1];
83        double errorDot = dxDot * dxDot + dyDot * dyDot;
84        maxErrorDot = JdkMath.max(maxErrorDot, errorDot);
85      }
86  
87      Assert.assertEquals(0.0, maxError,    1.0e-9);
88      Assert.assertEquals(0.0, maxErrorDot, 4.0e-7);
89    }
90  
91    @Test
92    public void testModelsMerging() throws MaxCountExceededException, MathIllegalArgumentException {
93  
94        // theoretical solution: y[0] = cos(t), y[1] = sin(t)
95        FirstOrderDifferentialEquations problem =
96            new FirstOrderDifferentialEquations() {
97                @Override
98              public void computeDerivatives(double t, double[] y, double[] dot) {
99                    dot[0] = -y[1];
100                   dot[1] =  y[0];
101               }
102               @Override
103             public int getDimension() {
104                   return 2;
105               }
106           };
107 
108       // integrate backward from &pi; to 0;
109       ContinuousOutputModel cm1 = new ContinuousOutputModel();
110       FirstOrderIntegrator integ1 =
111           new DormandPrince853Integrator(0, 1.0, 1.0e-8, 1.0e-8);
112       integ1.addStepHandler(cm1);
113       integ1.integrate(problem, JdkMath.PI, new double[] { -1.0, 0.0 },
114                        0, new double[2]);
115 
116       // integrate backward from 2&pi; to &pi;
117       ContinuousOutputModel cm2 = new ContinuousOutputModel();
118       FirstOrderIntegrator integ2 =
119           new DormandPrince853Integrator(0, 0.1, 1.0e-12, 1.0e-12);
120       integ2.addStepHandler(cm2);
121       integ2.integrate(problem, 2.0 * JdkMath.PI, new double[] { 1.0, 0.0 },
122                        JdkMath.PI, new double[2]);
123 
124       // merge the two half circles
125       ContinuousOutputModel cm = new ContinuousOutputModel();
126       cm.append(cm2);
127       cm.append(new ContinuousOutputModel());
128       cm.append(cm1);
129 
130       // check circle
131       Assert.assertEquals(2.0 * JdkMath.PI, cm.getInitialTime(), 1.0e-12);
132       Assert.assertEquals(0, cm.getFinalTime(), 1.0e-12);
133       Assert.assertEquals(cm.getFinalTime(), cm.getInterpolatedTime(), 1.0e-12);
134       for (double t = 0; t < 2.0 * JdkMath.PI; t += 0.1) {
135           cm.setInterpolatedTime(t);
136           double[] y = cm.getInterpolatedState();
137           Assert.assertEquals(JdkMath.cos(t), y[0], 1.0e-7);
138           Assert.assertEquals(JdkMath.sin(t), y[1], 1.0e-7);
139       }
140   }
141 
142   @Test
143   public void testErrorConditions() throws MaxCountExceededException, MathIllegalArgumentException {
144 
145       ContinuousOutputModel cm = new ContinuousOutputModel();
146       cm.handleStep(buildInterpolator(0, new double[] { 0.0, 1.0, -2.0 }, 1), true);
147 
148       // dimension mismatch
149       Assert.assertTrue(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0 }, 2.0));
150 
151       // hole between time ranges
152       Assert.assertTrue(checkAppendError(cm, 10.0, new double[] { 0.0, 1.0, -2.0 }, 20.0));
153 
154       // propagation direction mismatch
155       Assert.assertTrue(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0, -2.0 }, 0.0));
156 
157       // no errors
158       Assert.assertFalse(checkAppendError(cm, 1.0, new double[] { 0.0, 1.0, -2.0 }, 2.0));
159   }
160 
161   private boolean checkAppendError(ContinuousOutputModel cm,
162                                    double t0, double[] y0, double t1)
163       throws MaxCountExceededException, MathIllegalArgumentException {
164       try {
165           ContinuousOutputModel otherCm = new ContinuousOutputModel();
166           otherCm.handleStep(buildInterpolator(t0, y0, t1), true);
167           cm.append(otherCm);
168       } catch(MathIllegalArgumentException iae) {
169           return true; // there was an allowable error
170       }
171       return false; // no allowable error
172   }
173 
174     private StepInterpolator buildInterpolator(double t0, double[] y0, double t1) {
175         DummyStepInterpolator interpolator  = new DummyStepInterpolator(y0, new double[y0.length], t1 >= t0);
176         interpolator.storeTime(t0);
177         interpolator.shift();
178         interpolator.storeTime(t1);
179         return interpolator;
180     }
181 
182     public void checkValue(double value, double reference) {
183         Assert.assertTrue(JdkMath.abs(value - reference) < 1.0e-10);
184     }
185 
186     @Before
187     public void setUp() {
188         pb = new TestProblem3(0.9);
189         double minStep = 0;
190         double maxStep = pb.getFinalTime() - pb.getInitialTime();
191         integ = new DormandPrince54Integrator(minStep, maxStep, 1.0e-8, 1.0e-8);
192     }
193 
194     @After
195     public void tearDown() {
196         pb = null;
197         integ = null;
198     }
199 
200     private TestProblem3 pb;
201     private FirstOrderIntegrator integ;
202 }