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 static org.junit.Assert.assertTrue;
021    
022    import java.io.ByteArrayInputStream;
023    import java.io.ByteArrayOutputStream;
024    import java.io.IOException;
025    import java.io.ObjectInputStream;
026    import java.io.ObjectOutputStream;
027    import java.util.Random;
028    
029    import org.apache.commons.math.ode.ContinuousOutputModel;
030    import org.apache.commons.math.ode.DerivativeException;
031    import org.apache.commons.math.ode.IntegratorException;
032    import org.apache.commons.math.ode.TestProblem1;
033    import org.apache.commons.math.ode.TestProblem3;
034    import org.apache.commons.math.ode.sampling.StepHandler;
035    import org.apache.commons.math.ode.sampling.StepInterpolatorTestUtils;
036    import org.junit.Test;
037    
038    public class EulerStepInterpolatorTest {
039    
040      @Test
041      public void noReset() throws DerivativeException {
042    
043        double[]   y    =   { 0.0, 1.0, -2.0 };
044        double[][] yDot = { { 1.0, 2.0, -2.0 } };
045        EulerStepInterpolator interpolator = new EulerStepInterpolator();
046        interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
047        interpolator.storeTime(0);
048        interpolator.shift();
049        interpolator.storeTime(1);
050    
051        double[] result = interpolator.getInterpolatedState();
052        for (int i = 0; i < result.length; ++i) {
053          assertTrue(Math.abs(result[i] - y[i]) < 1.0e-10);
054        }
055    
056      }
057    
058      @Test
059      public void interpolationAtBounds()
060        throws DerivativeException {
061    
062        double   t0 = 0;
063        double[] y0 = {0.0, 1.0, -2.0};
064    
065        double[] y = y0.clone();
066        double[][] yDot = { new double[y0.length] };
067        EulerStepInterpolator interpolator = new EulerStepInterpolator();
068        interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
069        interpolator.storeTime(t0);
070    
071        double dt = 1.0;
072        y[0] =  1.0;
073        y[1] =  3.0;
074        y[2] = -4.0;
075        yDot[0][0] = (y[0] - y0[0]) / dt;
076        yDot[0][1] = (y[1] - y0[1]) / dt;
077        yDot[0][2] = (y[2] - y0[2]) / dt;
078        interpolator.shift();
079        interpolator.storeTime(t0 + dt);
080    
081        interpolator.setInterpolatedTime(interpolator.getPreviousTime());
082        double[] result = interpolator.getInterpolatedState();
083        for (int i = 0; i < result.length; ++i) {
084          assertTrue(Math.abs(result[i] - y0[i]) < 1.0e-10);
085        }
086    
087        interpolator.setInterpolatedTime(interpolator.getCurrentTime());
088        result = interpolator.getInterpolatedState();
089        for (int i = 0; i < result.length; ++i) {
090          assertTrue(Math.abs(result[i] - y[i]) < 1.0e-10);
091        }
092    
093      }
094    
095      @Test
096      public void interpolationInside()
097      throws DerivativeException {
098    
099        double[]   y    =   { 1.0, 3.0, -4.0 };
100        double[][] yDot = { { 1.0, 2.0, -2.0 } };
101        EulerStepInterpolator interpolator = new EulerStepInterpolator();
102        interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
103        interpolator.storeTime(0);
104        interpolator.shift();
105        interpolator.storeTime(1);
106    
107        interpolator.setInterpolatedTime(0.1);
108        double[] result = interpolator.getInterpolatedState();
109        assertTrue(Math.abs(result[0] - 0.1) < 1.0e-10);
110        assertTrue(Math.abs(result[1] - 1.2) < 1.0e-10);
111        assertTrue(Math.abs(result[2] + 2.2) < 1.0e-10);
112    
113        interpolator.setInterpolatedTime(0.5);
114        result = interpolator.getInterpolatedState();
115        assertTrue(Math.abs(result[0] - 0.5) < 1.0e-10);
116        assertTrue(Math.abs(result[1] - 2.0) < 1.0e-10);
117        assertTrue(Math.abs(result[2] + 3.0) < 1.0e-10);
118    
119      }
120    
121      @Test
122      public void derivativesConsistency()
123      throws DerivativeException, IntegratorException {
124        TestProblem3 pb = new TestProblem3();
125        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
126        EulerIntegrator integ = new EulerIntegrator(step);
127        StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
128      }
129    
130      @Test
131      public void serialization()
132        throws DerivativeException, IntegratorException,
133               IOException, ClassNotFoundException {
134    
135        TestProblem1 pb = new TestProblem1();
136        double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
137        EulerIntegrator integ = new EulerIntegrator(step);
138        integ.addStepHandler(new ContinuousOutputModel());
139        integ.integrate(pb,
140                        pb.getInitialTime(), pb.getInitialState(),
141                        pb.getFinalTime(), new double[pb.getDimension()]);
142    
143        ByteArrayOutputStream bos = new ByteArrayOutputStream();
144        ObjectOutputStream    oos = new ObjectOutputStream(bos);
145        for (StepHandler handler : integ.getStepHandlers()) {
146            oos.writeObject(handler);
147        }
148        
149        ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
150        ObjectInputStream     ois = new ObjectInputStream(bis);
151        ContinuousOutputModel cm  = (ContinuousOutputModel) ois.readObject();
152    
153        Random random = new Random(347588535632l);
154        double maxError = 0.0;
155        for (int i = 0; i < 1000; ++i) {
156          double r = random.nextDouble();
157          double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
158          cm.setInterpolatedTime(time);
159          double[] interpolatedY = cm.getInterpolatedState ();
160          double[] theoreticalY  = pb.computeTheoreticalState(time);
161          double dx = interpolatedY[0] - theoreticalY[0];
162          double dy = interpolatedY[1] - theoreticalY[1];
163          double error = dx * dx + dy * dy;
164          if (error > maxError) {
165            maxError = error;
166          }
167        }
168        assertTrue(maxError < 0.001);
169    
170      }
171    
172      private static class DummyIntegrator extends RungeKuttaIntegrator {
173    
174          
175          protected DummyIntegrator(RungeKuttaStepInterpolator prototype) {
176              super("dummy", new double[0], new double[0][0], new double[0], prototype, Double.NaN);
177          }
178    
179      }
180    
181    }