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 2015/12/08 20:34:08 UTC
[3/8] [math] Set up test framework for field-based embedded
Runge-Kutta integrators.
Set up test framework for field-based embedded Runge-Kutta integrators.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/a39a8527
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/a39a8527
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/a39a8527
Branch: refs/heads/field-ode
Commit: a39a8527ddc71711357cc9f1ff84771fe24d0109
Parents: d671ede
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 8 18:09:41 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 8 18:09:41 2015 +0100
----------------------------------------------------------------------
...ctEmbeddedRungeKuttaFieldIntegratorTest.java | 445 +++++++++++++++++++
1 file changed, 445 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a39a8527/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
new file mode 100644
index 0000000..06a05c9
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
@@ -0,0 +1,445 @@
+/*
+ * 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.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.FieldFirstOrderDifferentialEquations;
+import org.apache.commons.math3.ode.FieldFirstOrderIntegrator;
+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.TestFieldProblem3;
+import org.apache.commons.math3.ode.TestFieldProblem4;
+import org.apache.commons.math3.ode.TestFieldProblem5;
+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.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
+import org.junit.Assert;
+import org.junit.Test;
+
+public abstract class AbstractEmbeddedRungeKuttaFieldIntegratorTest {
+
+ protected abstract <T extends RealFieldElement<T>> EmbeddedRungeKuttaFieldIntegrator<T>
+ createIntegrator(Field<T> field, final double minStep, final double maxStep,
+ final double scalAbsoluteTolerance, final double scalRelativeTolerance) ;
+
+ protected abstract <T extends RealFieldElement<T>> EmbeddedRungeKuttaFieldIntegrator<T>
+ createIntegrator(Field<T> field, final double minStep, final double maxStep,
+ final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
+
+ @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
+ EmbeddedRungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, 0.001, 1.0, 1.0, 1.0);
+ 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 testForwardBackwardExceptions();
+
+ protected <T extends RealFieldElement<T>> void doTestForwardBackwardExceptions(final Field<T> field) {
+ FieldFirstOrderDifferentialEquations<T> equations = new FieldFirstOrderDifferentialEquations<T>() {
+
+ public int getDimension() {
+ return 1;
+ }
+
+ public void init(T t0, T[] y0, T t) {
+ }
+
+ public T[] computeDerivatives(T t, T[] y) {
+ if (t.getReal() < -0.5) {
+ throw new LocalException();
+ } else {
+ throw new RuntimeException("oops");
+ }
+ }
+ };
+
+ EmbeddedRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, 0.0, 1.0, 1.0e-10, 1.0e-10);
+
+ try {
+ integrator.integrate(new FieldExpandableODE<>(equations),
+ new FieldODEState<T>(field.getOne().negate(),
+ MathArrays.buildArray(field, 1)),
+ field.getZero());
+ Assert.fail("an exception should have been thrown");
+ } catch(LocalException de) {
+ // expected behavior
+ }
+
+ try {
+ integrator.integrate(new FieldExpandableODE<>(equations),
+ new FieldODEState<T>(field.getZero(),
+ MathArrays.buildArray(field, 1)),
+ field.getOne());
+ Assert.fail("an exception should have been thrown");
+ } catch(RuntimeException de) {
+ // expected behavior
+ }
+ }
+
+ protected static class LocalException extends RuntimeException {
+ private static final long serialVersionUID = 20151208L;
+ }
+
+ @Test(expected=NumberIsTooSmallException.class)
+ public abstract void testMinStep();
+
+ protected <T extends RealFieldElement<T>> void doTestMinStep(final Field<T> field)
+ throws NumberIsTooSmallException {
+
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
+ double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
+
+ FieldFirstOrderIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ vecAbsoluteTolerance, vecRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
+ Assert.fail("an exception should have been thrown");
+
+ }
+
+ @Test
+ public abstract void testIncreasingTolerance();
+
+ protected <T extends RealFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
+ double factor,
+ double epsilon) {
+
+ int previousCalls = Integer.MAX_VALUE;
+ for (int i = -12; i < -2; ++i) {
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double scalAbsoluteTolerance = FastMath.pow(10.0, i);
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ FieldFirstOrderIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ scalAbsoluteTolerance, scalRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertTrue(handler.getMaximalValueError().getReal() < (factor * scalAbsoluteTolerance));
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilon);
+
+ int calls = pb.getCalls();
+ Assert.assertEquals(integ.getEvaluations(), calls);
+ Assert.assertTrue(calls <= previousCalls);
+ previousCalls = calls;
+
+ }
+
+ }
+
+ @Test
+ public abstract void testEvents();
+
+ protected <T extends RealFieldElement<T>> void doTestEvents(final Field<T> field,
+ final double epsilonMaxValue,
+ final String name) {
+
+ TestFieldProblem4<T> pb = new TestFieldProblem4<T>(field);
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double scalAbsoluteTolerance = 1.0e-8;
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ FieldFirstOrderIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ scalAbsoluteTolerance, scalRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ FieldEventHandler<T>[] functions = pb.getEventsHandlers();
+ double convergence = 1.0e-8 * maxStep;
+ for (int l = 0; l < functions.length; ++l) {
+ integ.addEventHandler(functions[l], Double.POSITIVE_INFINITY, convergence, 1000);
+ }
+ Assert.assertEquals(functions.length, integ.getEventHandlers().size());
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), convergence);
+ Assert.assertEquals(12.0, handler.getLastTime().getReal(), convergence);
+ Assert.assertEquals(name, integ.getName());
+ integ.clearEventHandlers();
+ Assert.assertEquals(0, integ.getEventHandlers().size());
+
+ }
+
+ @Test(expected=LocalException.class)
+ public abstract void testEventsErrors();
+
+ protected <T extends RealFieldElement<T>> void doTestEventsErrors(final Field<T> field)
+ throws LocalException {
+ final TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double scalAbsoluteTolerance = 1.0e-8;
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ FieldFirstOrderIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ scalAbsoluteTolerance, scalRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+
+ integ.addEventHandler(new FieldEventHandler<T>() {
+ public void init(FieldODEStateAndDerivative<T> state0, T t) {
+ }
+ public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
+ return Action.CONTINUE;
+ }
+ public T g(FieldODEStateAndDerivative<T> state) {
+ T middle = pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5);
+ T offset = state.getTime().subtract(middle);
+ if (offset.getReal() > 0) {
+ throw new LocalException();
+ }
+ return offset;
+ }
+ public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
+ return state;
+ }
+ }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
+
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ }
+
+ @Test
+ public abstract void testEventsNoConvergence();
+
+ protected <T extends RealFieldElement<T>> void doTestEventsNoConvergence(final Field<T> field){
+
+ final TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double scalAbsoluteTolerance = 1.0e-8;
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ FieldFirstOrderIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ scalAbsoluteTolerance, scalRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+
+ integ.addEventHandler(new FieldEventHandler<T>() {
+ public void init(FieldODEStateAndDerivative<T> state0, T t) {
+ }
+ public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
+ return Action.CONTINUE;
+ }
+ public T g(FieldODEStateAndDerivative<T> state) {
+ T middle = pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5);
+ T offset = state.getTime().subtract(middle);
+ return (offset.getReal() > 0) ? offset.add(0.5) : offset.subtract(0.5);
+ }
+ public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
+ return state;
+ }
+ }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
+
+ try {
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+ Assert.fail("an exception should have been thrown");
+ } catch (MaxCountExceededException mcee) {
+ // Expected.
+ }
+
+ }
+
+ @Test
+ public abstract void testSanityChecks();
+
+ protected <T extends RealFieldElement<T>> void doTestSanityChecks(Field<T> field) {
+ TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field);
+ try {
+ EmbeddedRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, 0,
+ pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
+ new double[4], new double[4]);
+ integrator.integrate(new FieldExpandableODE<>(pb),
+ new FieldODEState<T>(pb.getInitialState().getTime(),
+ MathArrays.buildArray(field, 6)),
+ pb.getFinalTime());
+ Assert.fail("an exception should have been thrown");
+ } catch(DimensionMismatchException ie) {
+ }
+ try {
+ EmbeddedRungeKuttaFieldIntegrator<T> integrator =
+ createIntegrator(field, 0,
+ pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
+ new double[2], new double[4]);
+ integrator.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
+ Assert.fail("an exception should have been thrown");
+ } catch(DimensionMismatchException ie) {
+ }
+ try {
+ EmbeddedRungeKuttaFieldIntegrator<T> integrator =
+ createIntegrator(field, 0,
+ pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
+ new double[4], new double[4]);
+ integrator.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getInitialState().getTime());
+ Assert.fail("an exception should have been thrown");
+ } catch(NumberIsTooSmallException ie) {
+ }
+ }
+
+ @Test
+ public abstract void testBackward();
+
+ protected <T extends RealFieldElement<T>> void doTestBackward(Field<T> field,
+ final double espilonLast,
+ final double epsilonMaxValue,
+ final double epsilonMaxTime,
+ final String name)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).abs().getReal();
+ double scalAbsoluteTolerance = 1.0e-8;
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ EmbeddedRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ scalAbsoluteTolerance,
+ scalRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertEquals(0, handler.getLastError().getReal(), espilonLast);
+ 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 epsilon) {
+
+ final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9));
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
+ double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
+
+ FieldFirstOrderIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
+ vecAbsoluteTolerance, vecRelativeTolerance);
+ integ.addStepHandler(new KeplerHandler<T>(pb, epsilon));
+ integ.integrate(new FieldExpandableODE<>(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 epsilon;
+ public KeplerHandler(TestFieldProblem3<T> pb, double epsilon) {
+ this.pb = pb;
+ 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(0.0, maxError.getReal(), epsilon);
+ }
+ }
+ }
+
+}