You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2016/01/02 20:07:05 UTC
[05/14] [math] Renamed abstract test classes to match build
environment filters.
http://git-wip-us.apache.org/repos/asf/commons-math/blob/29649959/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegratorAbstractTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegratorAbstractTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegratorAbstractTest.java
new file mode 100644
index 0000000..518690a
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegratorAbstractTest.java
@@ -0,0 +1,662 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+
+import java.lang.reflect.Array;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.ode.FieldExpandableODE;
+import org.apache.commons.math3.ode.FirstOrderFieldDifferentialEquations;
+import org.apache.commons.math3.ode.FieldODEState;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.TestFieldProblem1;
+import org.apache.commons.math3.ode.TestFieldProblem2;
+import org.apache.commons.math3.ode.TestFieldProblem3;
+import org.apache.commons.math3.ode.TestFieldProblem4;
+import org.apache.commons.math3.ode.TestFieldProblem5;
+import org.apache.commons.math3.ode.TestFieldProblem6;
+import org.apache.commons.math3.ode.TestFieldProblemAbstract;
+import org.apache.commons.math3.ode.TestFieldProblemHandler;
+import org.apache.commons.math3.ode.events.Action;
+import org.apache.commons.math3.ode.events.FieldEventHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
+import org.apache.commons.math3.ode.sampling.StepInterpolatorTestUtils;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
+import org.junit.Assert;
+import org.junit.Test;
+
+public abstract class RungeKuttaFieldIntegratorAbstractTest {
+
+ protected abstract <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T>
+ createIntegrator(Field<T> field, T step);
+
+ @Test
+ public abstract void testNonFieldIntegratorConsistency();
+
+ protected <T extends RealFieldElement<T>> void doTestNonFieldIntegratorConsistency(final Field<T> field) {
+ try {
+
+ // get the Butcher arrays from the field integrator
+ RungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, field.getZero().add(1));
+ T[][] fieldA = fieldIntegrator.getA();
+ T[] fieldB = fieldIntegrator.getB();
+ T[] fieldC = fieldIntegrator.getC();
+
+ String fieldName = fieldIntegrator.getClass().getName();
+ String regularName = fieldName.replaceAll("Field", "");
+
+ // get the Butcher arrays from the regular integrator
+ @SuppressWarnings("unchecked")
+ Class<RungeKuttaIntegrator> c = (Class<RungeKuttaIntegrator>) Class.forName(regularName);
+ java.lang.reflect.Field jlrFieldA = c.getDeclaredField("STATIC_A");
+ jlrFieldA.setAccessible(true);
+ double[][] regularA = (double[][]) jlrFieldA.get(null);
+ java.lang.reflect.Field jlrFieldB = c.getDeclaredField("STATIC_B");
+ jlrFieldB.setAccessible(true);
+ double[] regularB = (double[]) jlrFieldB.get(null);
+ java.lang.reflect.Field jlrFieldC = c.getDeclaredField("STATIC_C");
+ jlrFieldC.setAccessible(true);
+ double[] regularC = (double[]) jlrFieldC.get(null);
+
+ Assert.assertEquals(regularA.length, fieldA.length);
+ for (int i = 0; i < regularA.length; ++i) {
+ checkArray(regularA[i], fieldA[i]);
+ }
+ checkArray(regularB, fieldB);
+ checkArray(regularC, fieldC);
+
+ } catch (ClassNotFoundException cnfe) {
+ Assert.fail(cnfe.getLocalizedMessage());
+ } catch (IllegalAccessException iae) {
+ Assert.fail(iae.getLocalizedMessage());
+ } catch (IllegalArgumentException iae) {
+ Assert.fail(iae.getLocalizedMessage());
+ } catch (SecurityException se) {
+ Assert.fail(se.getLocalizedMessage());
+ } catch (NoSuchFieldException nsfe) {
+ Assert.fail(nsfe.getLocalizedMessage());
+ }
+ }
+
+ private <T extends RealFieldElement<T>> void checkArray(double[] regularArray, T[] fieldArray) {
+ Assert.assertEquals(regularArray.length, fieldArray.length);
+ for (int i = 0; i < regularArray.length; ++i) {
+ if (regularArray[i] == 0) {
+ Assert.assertTrue(0.0 == fieldArray[i].getReal());
+ } else {
+ Assert.assertEquals(regularArray[i], fieldArray[i].getReal(), FastMath.ulp(regularArray[i]));
+ }
+ }
+ }
+
+ @Test
+ public abstract void testMissedEndEvent();
+
+ protected <T extends RealFieldElement<T>> void doTestMissedEndEvent(final Field<T> field,
+ final double epsilonT, final double epsilonY)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ final T t0 = field.getZero().add(1878250320.0000029);
+ final T tEvent = field.getZero().add(1878250379.9999986);
+ final T[] k = MathArrays.buildArray(field, 3);
+ k[0] = field.getZero().add(1.0e-4);
+ k[1] = field.getZero().add(1.0e-5);
+ k[2] = field.getZero().add(1.0e-6);
+ FirstOrderFieldDifferentialEquations<T> ode = new FirstOrderFieldDifferentialEquations<T>() {
+
+ public int getDimension() {
+ return k.length;
+ }
+
+ public void init(T t0, T[] y0, T t) {
+ }
+
+ public T[] computeDerivatives(T t, T[] y) {
+ T[] yDot = MathArrays.buildArray(field, k.length);
+ for (int i = 0; i < y.length; ++i) {
+ yDot[i] = k[i].multiply(y[i]);
+ }
+ return yDot;
+ }
+ };
+
+ RungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(60.0));
+
+ T[] y0 = MathArrays.buildArray(field, k.length);
+ for (int i = 0; i < y0.length; ++i) {
+ y0[i] = field.getOne().add(i);
+ }
+
+ FieldODEStateAndDerivative<T> result = integrator.integrate(new FieldExpandableODE<T>(ode),
+ new FieldODEState<T>(t0, y0),
+ tEvent);
+ Assert.assertEquals(tEvent.getReal(), result.getTime().getReal(), epsilonT);
+ T[] y = result.getState();
+ for (int i = 0; i < y.length; ++i) {
+ Assert.assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(),
+ y[i].getReal(),
+ epsilonY);
+ }
+
+ integrator.addEventHandler(new FieldEventHandler<T>() {
+
+ public void init(FieldODEStateAndDerivative<T> state0, T t) {
+ }
+
+ public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
+ return state;
+ }
+
+ public T g(FieldODEStateAndDerivative<T> state) {
+ return state.getTime().subtract(tEvent);
+ }
+
+ public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
+ Assert.assertEquals(tEvent.getReal(), state.getTime().getReal(), epsilonT);
+ return Action.CONTINUE;
+ }
+ }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
+ result = integrator.integrate(new FieldExpandableODE<T>(ode),
+ new FieldODEState<T>(t0, y0),
+ tEvent.add(120));
+ Assert.assertEquals(tEvent.add(120).getReal(), result.getTime().getReal(), epsilonT);
+ y = result.getState();
+ for (int i = 0; i < y.length; ++i) {
+ Assert.assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(),
+ y[i].getReal(),
+ epsilonY);
+ }
+
+ }
+
+ @Test
+ public abstract void testSanityChecks();
+
+ protected <T extends RealFieldElement<T>> void doTestSanityChecks(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ RungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(0.01));
+ try {
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ integrator.integrate(new FieldExpandableODE<T>(pb),
+ new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension() + 10)),
+ field.getOne());
+ Assert.fail("an exception should have been thrown");
+ } catch(DimensionMismatchException ie) {
+ }
+ try {
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ integrator.integrate(new FieldExpandableODE<T>(pb),
+ new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension())),
+ field.getZero());
+ Assert.fail("an exception should have been thrown");
+ } catch(NumberIsTooSmallException ie) {
+ }
+ }
+
+ @Test
+ public abstract void testDecreasingSteps();
+
+ protected <T extends RealFieldElement<T>> void doTestDecreasingSteps(Field<T> field,
+ final double safetyValueFactor,
+ final double safetyTimeFactor,
+ final double epsilonT)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ @SuppressWarnings("unchecked")
+ TestFieldProblemAbstract<T>[] allProblems =
+ (TestFieldProblemAbstract<T>[]) Array.newInstance(TestFieldProblemAbstract.class, 6);
+ allProblems[0] = new TestFieldProblem1<T>(field);
+ allProblems[1] = new TestFieldProblem2<T>(field);
+ allProblems[2] = new TestFieldProblem3<T>(field);
+ allProblems[3] = new TestFieldProblem4<T>(field);
+ allProblems[4] = new TestFieldProblem5<T>(field);
+ allProblems[5] = new TestFieldProblem6<T>(field);
+ for (TestFieldProblemAbstract<T> pb : allProblems) {
+
+ T previousValueError = null;
+ T previousTimeError = null;
+ for (int i = 4; i < 10; ++i) {
+
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(FastMath.pow(2.0, -i));
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ FieldEventHandler<T>[] functions = pb.getEventsHandlers();
+ for (int l = 0; l < functions.length; ++l) {
+ integ.addEventHandler(functions[l],
+ Double.POSITIVE_INFINITY, 1.0e-6 * step.getReal(), 1000);
+ }
+ Assert.assertEquals(functions.length, integ.getEventHandlers().size());
+ FieldODEStateAndDerivative<T> stop = integ.integrate(new FieldExpandableODE<T>(pb),
+ pb.getInitialState(),
+ pb.getFinalTime());
+ if (functions.length == 0) {
+ Assert.assertEquals(pb.getFinalTime().getReal(), stop.getTime().getReal(), epsilonT);
+ }
+
+ T error = handler.getMaximalValueError();
+ if (i > 4) {
+ Assert.assertTrue(error.subtract(previousValueError.abs().multiply(safetyValueFactor)).getReal() < 0);
+ }
+ previousValueError = error;
+
+ T timeError = handler.getMaximalTimeError();
+ if (i > 4) {
+ Assert.assertTrue(timeError.subtract(previousTimeError.abs().multiply(safetyTimeFactor)).getReal() <= 0);
+ }
+ previousTimeError = timeError;
+
+ integ.clearEventHandlers();
+ Assert.assertEquals(0, integ.getEventHandlers().size());
+ }
+
+ }
+
+ }
+
+ @Test
+ public abstract void testSmallStep();
+
+ protected <T extends RealFieldElement<T>> void doTestSmallStep(Field<T> field,
+ final double epsilonLast,
+ final double epsilonMaxValue,
+ final double epsilonMaxTime,
+ final String name)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertEquals(0, handler.getLastError().getReal(), epsilonLast);
+ Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
+ Assert.assertEquals(name, integ.getName());
+
+ }
+
+ @Test
+ public abstract void testBigStep();
+
+ protected <T extends RealFieldElement<T>> void doTestBigStep(Field<T> field,
+ final double belowLast,
+ final double belowMaxValue,
+ final double epsilonMaxTime,
+ final String name)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2);
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertTrue(handler.getLastError().getReal() > belowLast);
+ Assert.assertTrue(handler.getMaximalValueError().getReal() > belowMaxValue);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
+ Assert.assertEquals(name, integ.getName());
+
+ }
+
+ @Test
+ public abstract void testBackward();
+
+ protected <T extends RealFieldElement<T>> void doTestBackward(Field<T> field,
+ final double epsilonLast,
+ final double epsilonMaxValue,
+ final double epsilonMaxTime,
+ final String name)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs();
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertEquals(0, handler.getLastError().getReal(), epsilonLast);
+ Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
+ Assert.assertEquals(name, integ.getName());
+
+ }
+
+ @Test
+ public abstract void testKepler();
+
+ protected <T extends RealFieldElement<T>> void doTestKepler(Field<T> field, double expectedMaxError, double epsilon)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9));
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003);
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ integ.addStepHandler(new KeplerHandler<T>(pb, expectedMaxError, epsilon));
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+ }
+
+ private static class KeplerHandler<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
+ private T maxError;
+ private final TestFieldProblem3<T> pb;
+ private final double expectedMaxError;
+ private final double epsilon;
+ public KeplerHandler(TestFieldProblem3<T> pb, double expectedMaxError, double epsilon) {
+ this.pb = pb;
+ this.expectedMaxError = expectedMaxError;
+ this.epsilon = epsilon;
+ maxError = pb.getField().getZero();
+ }
+ public void init(FieldODEStateAndDerivative<T> state0, T t) {
+ maxError = pb.getField().getZero();
+ }
+ public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
+ throws MaxCountExceededException {
+
+ FieldODEStateAndDerivative<T> current = interpolator.getCurrentState();
+ T[] theoreticalY = pb.computeTheoreticalState(current.getTime());
+ T dx = current.getState()[0].subtract(theoreticalY[0]);
+ T dy = current.getState()[1].subtract(theoreticalY[1]);
+ T error = dx.multiply(dx).add(dy.multiply(dy));
+ if (error.subtract(maxError).getReal() > 0) {
+ maxError = error;
+ }
+ if (isLast) {
+ Assert.assertEquals(expectedMaxError, maxError.getReal(), epsilon);
+ }
+ }
+ }
+
+ @Test
+ public abstract void testStepSize();
+
+ protected <T extends RealFieldElement<T>> void doTestStepSize(final Field<T> field, final double epsilon)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ final T step = field.getZero().add(1.23456);
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ integ.addStepHandler(new FieldStepHandler<T>() {
+ public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) {
+ if (! isLast) {
+ Assert.assertEquals(step.getReal(),
+ interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).getReal(),
+ epsilon);
+ }
+ }
+ public void init(FieldODEStateAndDerivative<T> s0, T t) {
+ }
+ });
+ integ.integrate(new FieldExpandableODE<T>(new FirstOrderFieldDifferentialEquations<T>() {
+ public void init(T t0, T[] y0, T t) {
+ }
+ public T[] computeDerivatives(T t, T[] y) {
+ T[] dot = MathArrays.buildArray(t.getField(), 1);
+ dot[0] = t.getField().getOne();
+ return dot;
+ }
+ public int getDimension() {
+ return 1;
+ }
+ }), new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)), field.getZero().add(5.0));
+ }
+
+ @Test
+ public abstract void testSingleStep();
+
+ protected <T extends RealFieldElement<T>> void doTestSingleStep(final Field<T> field, final double epsilon) {
+
+ final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9));
+ T h = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003);
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(Double.NaN));
+ T t = pb.getInitialState().getTime();
+ T[] y = pb.getInitialState().getState();
+ for (int i = 0; i < 100; ++i) {
+ y = integ.singleStep(pb, t, y, t.add(h));
+ t = t.add(h);
+ }
+ T[] yth = pb.computeTheoreticalState(t);
+ T dx = y[0].subtract(yth[0]);
+ T dy = y[1].subtract(yth[1]);
+ T error = dx.multiply(dx).add(dy.multiply(dy));
+ Assert.assertEquals(0.0, error.getReal(), epsilon);
+ }
+
+ @Test
+ public abstract void testTooLargeFirstStep();
+
+ protected <T extends RealFieldElement<T>> void doTestTooLargeFirstStep(final Field<T> field) {
+
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.5));
+ final T t0 = field.getZero();
+ final T[] y0 = MathArrays.buildArray(field, 1);
+ y0[0] = field.getOne();
+ final T t = field.getZero().add(0.001);
+ FirstOrderFieldDifferentialEquations<T> equations = new FirstOrderFieldDifferentialEquations<T>() {
+
+ public int getDimension() {
+ return 1;
+ }
+
+ public void init(T t0, T[] y0, T t) {
+ }
+
+ public T[] computeDerivatives(T t, T[] y) {
+ Assert.assertTrue(t.getReal() >= FastMath.nextAfter(t0.getReal(), Double.NEGATIVE_INFINITY));
+ Assert.assertTrue(t.getReal() <= FastMath.nextAfter(t.getReal(), Double.POSITIVE_INFINITY));
+ T[] yDot = MathArrays.buildArray(field, 1);
+ yDot[0] = y[0].multiply(-100.0);
+ return yDot;
+ }
+
+ };
+
+ integ.integrate(new FieldExpandableODE<T>(equations), new FieldODEState<T>(t0, y0), t);
+
+ }
+
+ @Test
+ public abstract void testUnstableDerivative();
+
+ protected <T extends RealFieldElement<T>> void doTestUnstableDerivative(Field<T> field, double epsilon) {
+ final StepFieldProblem<T> stepProblem = new StepFieldProblem<T>(field,
+ field.getZero().add(0.0),
+ field.getZero().add(1.0),
+ field.getZero().add(2.0));
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.3));
+ integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
+ FieldODEStateAndDerivative<T> result = integ.integrate(new FieldExpandableODE<T>(stepProblem),
+ new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)),
+ field.getZero().add(10.0));
+ Assert.assertEquals(8.0, result.getState()[0].getReal(), epsilon);
+ }
+
+ @Test
+ public abstract void testDerivativesConsistency();
+
+ protected <T extends RealFieldElement<T>> void doTestDerivativesConsistency(final Field<T> field, double epsilon) {
+ TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
+ RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ }
+
+ @Test
+ public abstract void testPartialDerivatives();
+
+ protected void doTestPartialDerivatives(final double epsilonY,
+ final double[] epsilonPartials) {
+
+ // parameters indices
+ final int parameters = 5;
+ final int order = 1;
+ final int parOmega = 0;
+ final int parTO = 1;
+ final int parY00 = 2;
+ final int parY01 = 3;
+ final int parT = 4;
+
+ DerivativeStructure omega = new DerivativeStructure(parameters, order, parOmega, 1.3);
+ DerivativeStructure t0 = new DerivativeStructure(parameters, order, parTO, 1.3);
+ DerivativeStructure[] y0 = new DerivativeStructure[] {
+ new DerivativeStructure(parameters, order, parY00, 3.0),
+ new DerivativeStructure(parameters, order, parY01, 4.0)
+ };
+ DerivativeStructure t = new DerivativeStructure(parameters, order, parT, 6.0);
+ SinCos sinCos = new SinCos(omega);
+
+ RungeKuttaFieldIntegrator<DerivativeStructure> integrator =
+ createIntegrator(omega.getField(), t.subtract(t0).multiply(0.001));
+ FieldODEStateAndDerivative<DerivativeStructure> result =
+ integrator.integrate(new FieldExpandableODE<DerivativeStructure>(sinCos),
+ new FieldODEState<DerivativeStructure>(t0, y0),
+ t);
+
+ // check values
+ for (int i = 0; i < sinCos.getDimension(); ++i) {
+ Assert.assertEquals(sinCos.theoreticalY(t.getReal())[i], result.getState()[i].getValue(), epsilonY);
+ }
+
+ // check derivatives
+ final double[][] derivatives = sinCos.getDerivatives(t.getReal());
+ for (int i = 0; i < sinCos.getDimension(); ++i) {
+ for (int parameter = 0; parameter < parameters; ++parameter) {
+ Assert.assertEquals(derivatives[i][parameter],
+ dYdP(result.getState()[i], parameter),
+ epsilonPartials[parameter]);
+ }
+ }
+
+ }
+
+ private double dYdP(final DerivativeStructure y, final int parameter) {
+ int[] orders = new int[y.getFreeParameters()];
+ orders[parameter] = 1;
+ return y.getPartialDerivative(orders);
+ }
+
+ private static class SinCos implements FirstOrderFieldDifferentialEquations<DerivativeStructure> {
+
+ private final DerivativeStructure omega;
+ private DerivativeStructure r;
+ private DerivativeStructure alpha;
+
+ private double dRdY00;
+ private double dRdY01;
+ private double dAlphadOmega;
+ private double dAlphadT0;
+ private double dAlphadY00;
+ private double dAlphadY01;
+
+ protected SinCos(final DerivativeStructure omega) {
+ this.omega = omega;
+ }
+
+ public int getDimension() {
+ return 2;
+ }
+
+ public void init(final DerivativeStructure t0, final DerivativeStructure[] y0,
+ final DerivativeStructure finalTime) {
+
+ // theoretical solution is y(t) = { r * sin(omega * t + alpha), r * cos(omega * t + alpha) }
+ // so we retrieve alpha by identification from the initial state
+ final DerivativeStructure r2 = y0[0].multiply(y0[0]).add(y0[1].multiply(y0[1]));
+
+ this.r = r2.sqrt();
+ this.dRdY00 = y0[0].divide(r).getReal();
+ this.dRdY01 = y0[1].divide(r).getReal();
+
+ this.alpha = y0[0].atan2(y0[1]).subtract(t0.multiply(omega));
+ this.dAlphadOmega = -t0.getReal();
+ this.dAlphadT0 = -omega.getReal();
+ this.dAlphadY00 = y0[1].divide(r2).getReal();
+ this.dAlphadY01 = y0[0].negate().divide(r2).getReal();
+
+ }
+
+ public DerivativeStructure[] computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y) {
+ return new DerivativeStructure[] {
+ omega.multiply(y[1]),
+ omega.multiply(y[0]).negate()
+ };
+ }
+
+ public double[] theoreticalY(final double t) {
+ final double theta = omega.getReal() * t + alpha.getReal();
+ return new double[] {
+ r.getReal() * FastMath.sin(theta), r.getReal() * FastMath.cos(theta)
+ };
+ }
+
+ public double[][] getDerivatives(final double t) {
+
+ // intermediate angle and state
+ final double theta = omega.getReal() * t + alpha.getReal();
+ final double sin = FastMath.sin(theta);
+ final double cos = FastMath.cos(theta);
+ final double y0 = r.getReal() * sin;
+ final double y1 = r.getReal() * cos;
+
+ // partial derivatives of the state first component
+ final double dY0dOmega = y1 * (t + dAlphadOmega);
+ final double dY0dT0 = y1 * dAlphadT0;
+ final double dY0dY00 = dRdY00 * sin + y1 * dAlphadY00;
+ final double dY0dY01 = dRdY01 * sin + y1 * dAlphadY01;
+ final double dY0dT = y1 * omega.getReal();
+
+ // partial derivatives of the state second component
+ final double dY1dOmega = - y0 * (t + dAlphadOmega);
+ final double dY1dT0 = - y0 * dAlphadT0;
+ final double dY1dY00 = dRdY00 * cos - y0 * dAlphadY00;
+ final double dY1dY01 = dRdY01 * cos - y0 * dAlphadY01;
+ final double dY1dT = - y0 * omega.getReal();
+
+ return new double[][] {
+ { dY0dOmega, dY0dT0, dY0dY00, dY0dY01, dY0dT },
+ { dY1dOmega, dY1dT0, dY1dY00, dY1dY01, dY1dT }
+ };
+
+ }
+
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/29649959/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolatorAbstractTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolatorAbstractTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolatorAbstractTest.java
new file mode 100644
index 0000000..601e7b0
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolatorAbstractTest.java
@@ -0,0 +1,305 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.ode.nonstiff;
+
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.ode.AbstractIntegrator;
+import org.apache.commons.math3.ode.EquationsMapper;
+import org.apache.commons.math3.ode.ExpandableStatefulODE;
+import org.apache.commons.math3.ode.FieldEquationsMapper;
+import org.apache.commons.math3.ode.FieldExpandableODE;
+import org.apache.commons.math3.ode.FirstOrderFieldDifferentialEquations;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
+import org.junit.Assert;
+import org.junit.Test;
+
+public abstract class RungeKuttaFieldStepInterpolatorAbstractTest {
+
+ protected abstract <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
+ createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
+ FieldODEStateAndDerivative<T> globalPreviousState,
+ FieldODEStateAndDerivative<T> globalCurrentState,
+ FieldODEStateAndDerivative<T> softPreviousState,
+ FieldODEStateAndDerivative<T> softCurrentState,
+ FieldEquationsMapper<T> mapper);
+
+ protected abstract <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
+ createButcherArrayProvider(final Field<T> field);
+
+ @Test
+ public abstract void interpolationAtBounds();
+
+ protected <T extends RealFieldElement<T>> void doInterpolationAtBounds(final Field<T> field, double epsilon) {
+
+ RungeKuttaFieldStepInterpolator<T> interpolator = setUpInterpolator(field,
+ new SinCos<T>(field),
+ 0.0, new double[] { 0.0, 1.0 }, 0.125);
+
+ Assert.assertEquals(0.0, interpolator.getPreviousState().getTime().getReal(), 1.0e-15);
+ for (int i = 0; i < 2; ++i) {
+ Assert.assertEquals(interpolator.getPreviousState().getState()[i].getReal(),
+ interpolator.getInterpolatedState(interpolator.getPreviousState().getTime()).getState()[i].getReal(),
+ epsilon);
+ }
+ Assert.assertEquals(0.125, interpolator.getCurrentState().getTime().getReal(), 1.0e-15);
+ for (int i = 0; i < 2; ++i) {
+ Assert.assertEquals(interpolator.getCurrentState().getState()[i].getReal(),
+ interpolator.getInterpolatedState(interpolator.getCurrentState().getTime()).getState()[i].getReal(),
+ epsilon);
+ }
+
+ }
+
+ @Test
+ public abstract void interpolationInside();
+
+ protected <T extends RealFieldElement<T>> void doInterpolationInside(final Field<T> field,
+ double epsilonSin, double epsilonCos) {
+
+ RungeKuttaFieldStepInterpolator<T> interpolator = setUpInterpolator(field,
+ new SinCos<T>(field),
+ 0.0, new double[] { 0.0, 1.0 }, 0.0125);
+
+ int n = 100;
+ double maxErrorSin = 0;
+ double maxErrorCos = 0;
+ for (int i = 0; i <= n; ++i) {
+ T t = interpolator.getPreviousState().getTime().multiply(n - i).
+ add(interpolator.getCurrentState().getTime().multiply(i)).
+ divide(n);
+ FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(t);
+ maxErrorSin = FastMath.max(maxErrorSin, state.getState()[0].subtract(t.sin()).abs().getReal());
+ maxErrorCos = FastMath.max(maxErrorCos, state.getState()[1].subtract(t.cos()).abs().getReal());
+ }
+ Assert.assertEquals(0.0, maxErrorSin, epsilonSin);
+ Assert.assertEquals(0.0, maxErrorCos, epsilonCos);
+
+ }
+
+ @Test
+ public abstract void nonFieldInterpolatorConsistency();
+
+ protected <T extends RealFieldElement<T>> void doNonFieldInterpolatorConsistency(final Field<T> field,
+ double epsilonSin, double epsilonCos,
+ double epsilonSinDot, double epsilonCosDot) {
+
+ FirstOrderFieldDifferentialEquations<T> eqn = new SinCos<T>(field);
+ RungeKuttaFieldStepInterpolator<T> fieldInterpolator =
+ setUpInterpolator(field, eqn, 0.0, new double[] { 0.0, 1.0 }, 0.125);
+ RungeKuttaStepInterpolator regularInterpolator = convertInterpolator(fieldInterpolator, eqn);
+
+ int n = 100;
+ double maxErrorSin = 0;
+ double maxErrorCos = 0;
+ double maxErrorSinDot = 0;
+ double maxErrorCosDot = 0;
+ for (int i = 0; i <= n; ++i) {
+
+ T t = fieldInterpolator.getPreviousState().getTime().multiply(n - i).
+ add(fieldInterpolator.getCurrentState().getTime().multiply(i)).
+ divide(n);
+
+ FieldODEStateAndDerivative<T> state = fieldInterpolator.getInterpolatedState(t);
+ T[] fieldY = state.getState();
+ T[] fieldYDot = state.getDerivative();
+
+ regularInterpolator.setInterpolatedTime(t.getReal());
+ double[] regularY = regularInterpolator.getInterpolatedState();
+ double[] regularYDot = regularInterpolator.getInterpolatedDerivatives();
+
+ maxErrorSin = FastMath.max(maxErrorSin, fieldY[0].subtract(regularY[0]).abs().getReal());
+ maxErrorCos = FastMath.max(maxErrorCos, fieldY[1].subtract(regularY[1]).abs().getReal());
+ maxErrorSinDot = FastMath.max(maxErrorSinDot, fieldYDot[0].subtract(regularYDot[0]).abs().getReal());
+ maxErrorCosDot = FastMath.max(maxErrorCosDot, fieldYDot[1].subtract(regularYDot[1]).abs().getReal());
+
+ }
+ Assert.assertEquals(0.0, maxErrorSin, epsilonSin);
+ Assert.assertEquals(0.0, maxErrorCos, epsilonCos);
+ Assert.assertEquals(0.0, maxErrorSinDot, epsilonSinDot);
+ Assert.assertEquals(0.0, maxErrorCosDot, epsilonCosDot);
+
+ }
+
+ private <T extends RealFieldElement<T>>
+ RungeKuttaFieldStepInterpolator<T> setUpInterpolator(final Field<T> field,
+ final FirstOrderFieldDifferentialEquations<T> eqn,
+ final double t0, final double[] y0,
+ final double t1) {
+
+ // get the Butcher arrays from the field integrator
+ FieldButcherArrayProvider<T> provider = createButcherArrayProvider(field);
+ T[][] a = provider.getA();
+ T[] b = provider.getB();
+ T[] c = provider.getC();
+
+ // store initial state
+ T t = field.getZero().add(t0);
+ T[] fieldY = MathArrays.buildArray(field, eqn.getDimension());
+ T[][] fieldYDotK = MathArrays.buildArray(field, b.length, -1);
+ for (int i = 0; i < y0.length; ++i) {
+ fieldY[i] = field.getZero().add(y0[i]);
+ }
+ fieldYDotK[0] = eqn.computeDerivatives(t, fieldY);
+ FieldODEStateAndDerivative<T> s0 = new FieldODEStateAndDerivative<T>(t, fieldY, fieldYDotK[0]);
+
+ // perform one integration step, in order to get consistent derivatives
+ T h = field.getZero().add(t1 - t0);
+ for (int k = 0; k < a.length; ++k) {
+ for (int i = 0; i < y0.length; ++i) {
+ fieldY[i] = field.getZero().add(y0[i]);
+ for (int s = 0; s <= k; ++s) {
+ fieldY[i] = fieldY[i].add(h.multiply(a[k][s].multiply(fieldYDotK[s][i])));
+ }
+ }
+ fieldYDotK[k + 1] = eqn.computeDerivatives(h.multiply(c[k]).add(t0), fieldY);
+ }
+
+ // store state at step end
+ t = field.getZero().add(t1);
+ for (int i = 0; i < y0.length; ++i) {
+ fieldY[i] = field.getZero().add(y0[i]);
+ for (int s = 0; s < b.length; ++s) {
+ fieldY[i] = fieldY[i].add(h.multiply(b[s].multiply(fieldYDotK[s][i])));
+ }
+ }
+ FieldODEStateAndDerivative<T> s1 = new FieldODEStateAndDerivative<T>(t, fieldY,
+ eqn.computeDerivatives(t, fieldY));
+
+ return createInterpolator(field, t1 > t0, fieldYDotK, s0, s1, s0, s1,
+ new FieldExpandableODE<T>(eqn).getMapper());
+
+ }
+
+ private <T extends RealFieldElement<T>>
+ RungeKuttaStepInterpolator convertInterpolator(final RungeKuttaFieldStepInterpolator<T> fieldInterpolator,
+ final FirstOrderFieldDifferentialEquations<T> eqn) {
+
+ RungeKuttaStepInterpolator regularInterpolator = null;
+ try {
+
+ String interpolatorName = fieldInterpolator.getClass().getName();
+ String integratorName = interpolatorName.replaceAll("Field", "");
+ @SuppressWarnings("unchecked")
+ Class<RungeKuttaStepInterpolator> clz = (Class<RungeKuttaStepInterpolator>) Class.forName(integratorName);
+ regularInterpolator = clz.newInstance();
+
+ double[][] yDotArray = null;
+ java.lang.reflect.Field fYD = RungeKuttaFieldStepInterpolator.class.getDeclaredField("yDotK");
+ fYD.setAccessible(true);
+ @SuppressWarnings("unchecked")
+ T[][] fieldYDotk = (T[][]) fYD.get(fieldInterpolator);
+ yDotArray = new double[fieldYDotk.length][];
+ for (int i = 0; i < yDotArray.length; ++i) {
+ yDotArray[i] = new double[fieldYDotk[i].length];
+ for (int j = 0; j < yDotArray[i].length; ++j) {
+ yDotArray[i][j] = fieldYDotk[i][j].getReal();
+ }
+ }
+ double[] y = new double[yDotArray[0].length];
+
+ EquationsMapper primaryMapper = null;
+ EquationsMapper[] secondaryMappers = null;
+ java.lang.reflect.Field fMapper = AbstractFieldStepInterpolator.class.getDeclaredField("mapper");
+ fMapper.setAccessible(true);
+ @SuppressWarnings("unchecked")
+ FieldEquationsMapper<T> mapper = (FieldEquationsMapper<T>) fMapper.get(fieldInterpolator);
+ java.lang.reflect.Field fStart = FieldEquationsMapper.class.getDeclaredField("start");
+ fStart.setAccessible(true);
+ int[] start = (int[]) fStart.get(mapper);
+ primaryMapper = new EquationsMapper(start[0], start[1]);
+ secondaryMappers = new EquationsMapper[mapper.getNumberOfEquations() - 1];
+ for (int i = 0; i < secondaryMappers.length; ++i) {
+ secondaryMappers[i] = new EquationsMapper(start[i + 1], start[i + 2]);
+ }
+
+ AbstractIntegrator dummyIntegrator = new AbstractIntegrator("dummy") {
+ @Override
+ public void integrate(ExpandableStatefulODE equations, double t) {
+ Assert.fail("this method should not be called");
+ }
+ @Override
+ public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
+ T fieldT = fieldInterpolator.getCurrentState().getTime().getField().getZero().add(t);
+ T[] fieldY = MathArrays.buildArray(fieldInterpolator.getCurrentState().getTime().getField(), y.length);
+ for (int i = 0; i < y.length; ++i) {
+ fieldY[i] = fieldInterpolator.getCurrentState().getTime().getField().getZero().add(y[i]);
+ }
+ T[] fieldYDot = eqn.computeDerivatives(fieldT, fieldY);
+ for (int i = 0; i < yDot.length; ++i) {
+ yDot[i] = fieldYDot[i].getReal();
+ }
+ }
+ };
+ regularInterpolator.reinitialize(dummyIntegrator, y, yDotArray,
+ fieldInterpolator.isForward(),
+ primaryMapper, secondaryMappers);
+
+ T[] fieldPreviousY = fieldInterpolator.getPreviousState().getState();
+ for (int i = 0; i < y.length; ++i) {
+ y[i] = fieldPreviousY[i].getReal();
+ }
+ regularInterpolator.storeTime(fieldInterpolator.getPreviousState().getTime().getReal());
+
+ regularInterpolator.shift();
+
+ T[] fieldCurrentY = fieldInterpolator.getCurrentState().getState();
+ for (int i = 0; i < y.length; ++i) {
+ y[i] = fieldCurrentY[i].getReal();
+ }
+ regularInterpolator.storeTime(fieldInterpolator.getCurrentState().getTime().getReal());
+
+ } catch (ClassNotFoundException cnfe) {
+ Assert.fail(cnfe.getLocalizedMessage());
+ } catch (InstantiationException ie) {
+ Assert.fail(ie.getLocalizedMessage());
+ } catch (IllegalAccessException iae) {
+ Assert.fail(iae.getLocalizedMessage());
+ } catch (NoSuchFieldException nsfe) {
+ Assert.fail(nsfe.getLocalizedMessage());
+ } catch (IllegalArgumentException iae) {
+ Assert.fail(iae.getLocalizedMessage());
+ }
+
+ return regularInterpolator;
+
+ }
+
+ private static class SinCos<T extends RealFieldElement<T>> implements FirstOrderFieldDifferentialEquations<T> {
+ private final Field<T> field;
+ protected SinCos(final Field<T> field) {
+ this.field = field;
+ }
+ public int getDimension() {
+ return 2;
+ }
+ public void init(final T t0, final T[] y0, final T finalTime) {
+ }
+ public T[] computeDerivatives(final T t, final T[] y) {
+ T[] yDot = MathArrays.buildArray(field, 2);
+ yDot[0] = y[1];
+ yDot[1] = y[0].negate();
+ return yDot;
+ }
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/29649959/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegratorTest.java
index 700661d..971dac8 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegratorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegratorTest.java
@@ -22,7 +22,7 @@ import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.util.Decimal64Field;
-public class ThreeEighthesFieldIntegratorTest extends AbstractRungeKuttaFieldIntegratorTest {
+public class ThreeEighthesFieldIntegratorTest extends RungeKuttaFieldIntegratorAbstractTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T>
createIntegrator(Field<T> field, T step) {
http://git-wip-us.apache.org/repos/asf/commons-math/blob/29649959/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java
index f1a2d28..f326a2b 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java
@@ -25,7 +25,7 @@ import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
import org.apache.commons.math3.util.Decimal64Field;
import org.junit.Test;
-public class ThreeEighthesFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
+public class ThreeEighthesFieldStepInterpolatorTest extends RungeKuttaFieldStepInterpolatorAbstractTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,