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 &pi; 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&pi; to &pi;
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    }