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 }