commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject [17/50] [abbrv] [math] Set up test framework for field-based embedded Runge-Kutta integrators.
Date Wed, 06 Jan 2016 13:50:33 GMT
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/87edfd27
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/87edfd27
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/87edfd27

Branch: refs/heads/master
Commit: 87edfd2751ba8bf84f15aa1df288722c36ffcbf0
Parents: 7a5431e
Author: Luc Maisonobe <luc@apache.org>
Authored: Wed Jan 6 12:41:26 2016 +0100
Committer: Luc Maisonobe <luc@apache.org>
Committed: Wed Jan 6 12:41:26 2016 +0100

----------------------------------------------------------------------
 ...ctEmbeddedRungeKuttaFieldIntegratorTest.java | 445 +++++++++++++++++++
 1 file changed, 445 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/87edfd27/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
b/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractEmbeddedRungeKuttaFieldIntegratorTest.java
new file mode 100644
index 0000000..357d92a
--- /dev/null
+++ b/src/test/java/org/apache/commons/math4/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.math4.ode.nonstiff;
+
+
+import org.apache.commons.math4.Field;
+import org.apache.commons.math4.RealFieldElement;
+import org.apache.commons.math4.exception.DimensionMismatchException;
+import org.apache.commons.math4.exception.MaxCountExceededException;
+import org.apache.commons.math4.exception.NoBracketingException;
+import org.apache.commons.math4.exception.NumberIsTooSmallException;
+import org.apache.commons.math4.ode.FieldExpandableODE;
+import org.apache.commons.math4.ode.FieldFirstOrderDifferentialEquations;
+import org.apache.commons.math4.ode.FieldFirstOrderIntegrator;
+import org.apache.commons.math4.ode.FieldODEState;
+import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math4.ode.TestFieldProblem1;
+import org.apache.commons.math4.ode.TestFieldProblem3;
+import org.apache.commons.math4.ode.TestFieldProblem4;
+import org.apache.commons.math4.ode.TestFieldProblem5;
+import org.apache.commons.math4.ode.TestFieldProblemHandler;
+import org.apache.commons.math4.ode.events.Action;
+import org.apache.commons.math4.ode.events.FieldEventHandler;
+import org.apache.commons.math4.ode.sampling.FieldStepHandler;
+import org.apache.commons.math4.ode.sampling.FieldStepInterpolator;
+import org.apache.commons.math4.util.FastMath;
+import org.apache.commons.math4.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);
+            }
+        }
+    }
+
+}


Mime
View raw message