Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Factor test code.

This commit is contained in:
Luc Maisonobe 2016-01-06 12:40:58 +01:00
parent d53a6f834b
commit 85b941ea1f
4 changed files with 657 additions and 847 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

510
src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractRungeKuttaFieldIntegratorTest.java Normal file
View File

@ -0,0 +1,510 @@
/*
* 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 java.lang.reflect.Array;
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.FieldODEState;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.ode.TestFieldProblem1;
import org.apache.commons.math4.ode.TestFieldProblem2;
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.TestFieldProblem6;
import org.apache.commons.math4.ode.TestFieldProblemAbstract;
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 AbstractRungeKuttaFieldIntegratorTest {
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);
FieldFirstOrderDifferentialEquations<T> ode = new FieldFirstOrderDifferentialEquations<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<>(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<>(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 espilonLast,
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(), espilonLast);
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 espilonLast,
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<>(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 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<>(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 FieldFirstOrderDifferentialEquations<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);
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) {
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<>(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<>(stepProblem),
new FieldODEState<>(field.getZero(), MathArrays.buildArray(field, 1)),
field.getZero().add(10.0));
Assert.assertEquals(8.0, result.getState()[0].getReal(), epsilon);
}
}

440
src/test/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegratorTest.java
View File

@ -18,411 +18,77 @@
package org.apache.commons.math4.ode.nonstiff; package org.apache.commons.math4.ode.nonstiff;
import java.lang.reflect.Array;
import org.apache.commons.math4.Field; import org.apache.commons.math4.Field;
import org.apache.commons.math4.RealFieldElement; 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.TestFieldProblem2;
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.TestFieldProblem6;
import org.apache.commons.math4.ode.TestFieldProblemAbstract;
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.Decimal64Field; import org.apache.commons.math4.util.Decimal64Field;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathArrays;
import org.junit.Assert;
import org.junit.Test; import org.junit.Test;
public class ClassicalRungeKuttaFieldIntegratorTest { public class ClassicalRungeKuttaFieldIntegratorTest extends AbstractRungeKuttaFieldIntegratorTest {
@Test protected <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T>
public void testMissedEndEvent() createIntegrator(Field<T> field, T step) {
throws DimensionMismatchException, NumberIsTooSmallException, return new ClassicalRungeKuttaFieldIntegrator<T>(field, step);
MaxCountExceededException, NoBracketingException { }
doTestMissedEndEvent(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestMissedEndEvent(final Field<T> field) @Test
throws DimensionMismatchException, NumberIsTooSmallException, public void testNonFieldIntegratorConsistency() {
MaxCountExceededException, NoBracketingException { doTestNonFieldIntegratorConsistency(Decimal64Field.getInstance());
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);
FieldFirstOrderDifferentialEquations<T> ode = new FieldFirstOrderDifferentialEquations<T>() {
public int getDimension() { @Test
return k.length; public void testMissedEndEvent() {
} doTestMissedEndEvent(Decimal64Field.getInstance(), 5.0e-6, 1.0e-9);
}
public void init(T t0, T[] y0, T t) { @Test
} public void testSanityChecks() {
doTestSanityChecks(Decimal64Field.getInstance());
}
public T[] computeDerivatives(T t, T[] y) { @Test
T[] yDot = MathArrays.buildArray(field, k.length); public void testDecreasingSteps() {
for (int i = 0; i < y.length; ++i) { doTestDecreasingSteps(Decimal64Field.getInstance(), 1.0, 1.0, 1.0e-10);
yDot[i] = k[i].multiply(y[i]); }
}
return yDot;
}
};
ClassicalRungeKuttaFieldIntegrator<T> integrator = @Test
new ClassicalRungeKuttaFieldIntegrator<T>(field, field.getZero().add(60.0)); public void testSmallStep() {
doTestSmallStep(Decimal64Field.getInstance(), 2.0e-13, 4.0e-12, 1.0e-12, "classical Runge-Kutta");
}
T[] y0 = MathArrays.buildArray(field, k.length); @Test
for (int i = 0; i < y0.length; ++i) { public void testBigStep() {
y0[i] = field.getOne().add(i); doTestBigStep(Decimal64Field.getInstance(), 0.0004, 0.005, 1.0e-12, "classical Runge-Kutta");
}
FieldODEStateAndDerivative<T> result = integrator.integrate(new FieldExpandableODE<T>(ode), }
new FieldODEState<T>(t0, y0),
tEvent);
Assert.assertEquals(tEvent.getReal(), result.getTime().getReal(), 5.0e-6);
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(),
1.0e-9);
}
integrator.addEventHandler(new FieldEventHandler<T>() { @Test
public void testBackward() {
doTestBackward(Decimal64Field.getInstance(), 5.0e-10, 7.0e-10, 1.0e-12, "classical Runge-Kutta");
}
public void init(FieldODEStateAndDerivative<T> state0, T t) { @Test
} public void testKepler() {
doTestKepler(Decimal64Field.getInstance(), 5.82e-3, 1.0e-5);
}
public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) { @Test
return state; public void testStepSize() {
} doTestStepSize(Decimal64Field.getInstance(), 1.0e-12);
}
public T g(FieldODEStateAndDerivative<T> state) { @Test
return state.getTime().subtract(tEvent); public void testSingleStep() {
} doTestSingleStep(Decimal64Field.getInstance(), 9.3e-9);
}
public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) { @Test
Assert.assertEquals(tEvent.getReal(), state.getTime().getReal(), 5.0e-6); public void testTooLargeFirstStep() {
return Action.CONTINUE; doTestTooLargeFirstStep(Decimal64Field.getInstance());
} }
}, 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(), 5.0e-6);
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(),
1.0e-9);
}
} @Test
public void testUnstableDerivative() {
@Test doTestUnstableDerivative(Decimal64Field.getInstance(), 1.0e-12);
public void testSanityChecks() }
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestSanityChecks(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestSanityChecks(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
ClassicalRungeKuttaFieldIntegrator<T> integrator =
new ClassicalRungeKuttaFieldIntegrator<T>(field, field.getZero().add(0.01));
try {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
integrator.integrate(new FieldExpandableODE<>(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<>(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 void testDecreasingSteps()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestDecreasingSteps(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestDecreasingSteps(Field<T> field)
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));
FieldFirstOrderIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(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(), 1.0e-10);
}
T error = handler.getMaximalValueError();
if (i > 4) {
Assert.assertTrue(error.subtract(previousValueError.abs().multiply(1.01)).getReal() < 0);
}
previousValueError = error;
T timeError = handler.getMaximalTimeError();
if (i > 4) {
Assert.assertTrue(timeError.subtract(previousTimeError.abs()).getReal() <= 0);
}
previousTimeError = timeError;
integ.clearEventHandlers();
Assert.assertEquals(0, integ.getEventHandlers().size());
}
}
}
@Test
public void testSmallStep()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestSmallStep(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestSmallStep(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
FieldFirstOrderIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(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() < 2.0e-13);
Assert.assertTrue(handler.getMaximalValueError().getReal() < 4.0e-12);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
Assert.assertEquals("classical Runge-Kutta", integ.getName());
}
@Test
public void testBigStep()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestBigStep(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestBigStep(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2);
FieldFirstOrderIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(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() > 0.0004);
Assert.assertTrue(handler.getMaximalValueError().getReal() > 0.005);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
}
@Test
public void testBackward()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestBackward(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestBackward(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs();
FieldFirstOrderIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(field, step);
TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
integ.addStepHandler(handler);
integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
Assert.assertTrue(handler.getLastError().getReal() < 5.0e-10);
Assert.assertTrue(handler.getMaximalValueError().getReal() < 7.0e-10);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
Assert.assertEquals("classical Runge-Kutta", integ.getName());
}
@Test
public void testKepler()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestKepler(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestKepler(Field<T> field)
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);
FieldFirstOrderIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(field, step);
integ.addStepHandler(new KeplerHandler<T>(pb));
integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
}
private static class KeplerHandler<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
public KeplerHandler(TestFieldProblem3<T> pb) {
this.pb = pb;
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) {
// even with more than 1000 evaluations per period,
// RK4 is not able to integrate such an eccentric
// orbit with a good accuracy
Assert.assertTrue(maxError.getReal() > 0.005);
}
}
private T maxError;
private TestFieldProblem3<T> pb;
}
@Test
public void testStepSize()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestStepSize(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestStepSize(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
final T step = field.getZero().add(1.23456);
FieldFirstOrderIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(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(),
1.0e-12);
}
}
public void init(FieldODEStateAndDerivative<T> s0, T t) {
}
});
integ.integrate(new FieldExpandableODE<T>(new FieldFirstOrderDifferentialEquations<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 void testTooLargeFirstStep() {
doTestTooLargeFirstStep(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestTooLargeFirstStep(final Field<T> field) {
RungeKuttaFieldIntegrator<T> integ = new ClassicalRungeKuttaFieldIntegrator<T>(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);
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) {
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<>(equations), new FieldODEState<T>(t0, y0), t);
}
} }

236
src/test/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegratorTest.java
View File

@ -18,230 +18,78 @@
package org.apache.commons.math4.ode.nonstiff; package org.apache.commons.math4.ode.nonstiff;
import java.lang.reflect.Array;
import org.apache.commons.math4.Field; import org.apache.commons.math4.Field;
import org.apache.commons.math4.RealFieldElement; 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.TestFieldProblem2;
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.TestFieldProblem6;
import org.apache.commons.math4.ode.TestFieldProblemAbstract;
import org.apache.commons.math4.ode.TestFieldProblemHandler;
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.Decimal64Field; import org.apache.commons.math4.util.Decimal64Field;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathArrays;
import org.junit.Assert;
import org.junit.Test; import org.junit.Test;
public class EulerFieldIntegratorTest { public class EulerFieldIntegratorTest extends AbstractRungeKuttaFieldIntegratorTest {
@Test(expected=DimensionMismatchException.class) protected <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T>
public void testDimensionCheck() createIntegrator(Field<T> field, T step) {
throws DimensionMismatchException, NumberIsTooSmallException, return new EulerFieldIntegrator<T>(field, step);
MaxCountExceededException, NoBracketingException {
doTestDimensionCheck(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestDimensionCheck(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblemAbstract<T> pb = new TestFieldProblem1<T>(field);
FieldExpandableODE<T> ode = new FieldExpandableODE<T>(pb);
FieldFirstOrderIntegrator<T> integ =
new EulerFieldIntegrator<T>(field, field.getZero().add(0.01));
FieldODEState<T> start =
new FieldODEState<T>(field.getZero().add(0.0),
MathArrays.buildArray(field, pb.getDimension() + 10));
integ.integrate(ode, start, field.getZero().add(1.0));
} }
@Test @Test
public void testDecreasingSteps() public void testNonFieldIntegratorConsistency() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestNonFieldIntegratorConsistency(Decimal64Field.getInstance());
MaxCountExceededException, NoBracketingException {
dotTestDecreasingSteps(Decimal64Field.getInstance());
} }
private <T extends RealFieldElement<T>> void dotTestDecreasingSteps(Field<T> field) @Test
throws DimensionMismatchException, NumberIsTooSmallException, public void testMissedEndEvent() {
MaxCountExceededException, NoBracketingException { doTestMissedEndEvent(Decimal64Field.getInstance(), 1.0e-15, 6.0e-5);
}
@SuppressWarnings("unchecked") @Test
TestFieldProblemAbstract<T>[] allProblems = public void testSanityChecks() {
(TestFieldProblemAbstract<T>[]) Array.newInstance(TestFieldProblemAbstract.class, 6); doTestSanityChecks(Decimal64Field.getInstance());
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; @Test
T previousTimeError = null; public void testDecreasingSteps() {
for (int i = 4; i < 8; ++i) { doTestDecreasingSteps(Decimal64Field.getInstance(), 1.0, 1.5, 1.0e-10);
}
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(FastMath.pow(2.0, -i)); @Test
public void testSmallStep() {
doTestSmallStep(Decimal64Field.getInstance(), 2.0e-4, 1.0e-3, 1.0e-12, "Euler");
}
FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step); @Test
TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ); public void testBigStep() {
integ.addStepHandler(handler); doTestBigStep(Decimal64Field.getInstance(), 0.01, 0.2, 1.0e-12, "Euler");
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);
}
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(), 1.0e-10);
}
T valueError = handler.getMaximalValueError();
if (i > 4) {
Assert.assertTrue(valueError.subtract(previousValueError.abs()).getReal() < 0);
}
previousValueError = valueError;
T timeError = handler.getMaximalTimeError();
if (i > 4) {
Assert.assertTrue(timeError.subtract(previousTimeError.abs()).getReal() <= 0);
}
previousTimeError = timeError;
}
}
} }
@Test @Test
public void testSmallStep() public void testBackward() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestBackward(Decimal64Field.getInstance(),0.45, 0.45, 1.0e-12, "Euler");
MaxCountExceededException, NoBracketingException {
doTestSmallStep(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestSmallStep(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(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() < 2.0e-4);
Assert.assertTrue(handler.getMaximalValueError().getReal() < 1.0e-3);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
Assert.assertEquals("Euler", integ.getName());
} }
@Test @Test
public void testBigStep() public void testKepler() {
throws DimensionMismatchException, NumberIsTooSmallException, // Euler integrator is clearly not able to solve this problem
MaxCountExceededException, NoBracketingException { doTestKepler(Decimal64Field.getInstance(), 881.176, 0.001);
doTestBigStep(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestBigStep(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2);
FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(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() > 0.01);
Assert.assertTrue(handler.getMaximalValueError().getReal() > 0.2);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
} }
@Test @Test
public void testBackward() public void testStepSize() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestStepSize(Decimal64Field.getInstance(), 1.0e-12);
MaxCountExceededException, NoBracketingException {
doTestBackward(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestBackward(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs();
FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
integ.addStepHandler(handler);
integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
Assert.assertTrue(handler.getLastError().getReal() < 0.45);
Assert.assertTrue(handler.getMaximalValueError().getReal() < 0.45);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
Assert.assertEquals("Euler", integ.getName());
} }
@Test @Test
public void testStepSize() public void testSingleStep() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestSingleStep(Decimal64Field.getInstance(), 0.21);
MaxCountExceededException, NoBracketingException {
doTestStepSize(Decimal64Field.getInstance());
} }
private <T extends RealFieldElement<T>> void doTestStepSize(Field<T> field) @Test
throws DimensionMismatchException, NumberIsTooSmallException, public void testTooLargeFirstStep() {
MaxCountExceededException, NoBracketingException { doTestTooLargeFirstStep(Decimal64Field.getInstance());
final T step = field.getZero().add(1.23456); }
FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
integ.addStepHandler(new FieldStepHandler<T>() { @Test
public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) { public void testUnstableDerivative() {
if (! isLast) { doTestUnstableDerivative(Decimal64Field.getInstance(), 1.0e-12);
Assert.assertEquals(step.getReal(),
interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).getReal(),
1.0e-12);
}
}
public void init(FieldODEStateAndDerivative<T> s0, T t) {
}
});
integ.integrate(new FieldExpandableODE<T>(new FieldFirstOrderDifferentialEquations<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));
} }
} }

318
src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegratorTest.java
View File

@ -18,296 +18,82 @@
package org.apache.commons.math4.ode.nonstiff; package org.apache.commons.math4.ode.nonstiff;
import java.lang.reflect.Array;
import org.apache.commons.math4.Field; import org.apache.commons.math4.Field;
import org.apache.commons.math4.RealFieldElement; 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.FieldExpandableODE;
import org.apache.commons.math4.ode.FieldFirstOrderIntegrator;
import org.apache.commons.math4.ode.FieldODEState; import org.apache.commons.math4.ode.FieldODEState;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative; import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math4.ode.FirstOrderIntegrator;
import org.apache.commons.math4.ode.TestFieldProblem1;
import org.apache.commons.math4.ode.TestFieldProblem2;
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.TestFieldProblem6;
import org.apache.commons.math4.ode.TestFieldProblemAbstract;
import org.apache.commons.math4.ode.TestFieldProblemHandler;
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.ode.sampling.StepHandler;
import org.apache.commons.math4.ode.sampling.StepInterpolator;
import org.apache.commons.math4.util.Decimal64Field; import org.apache.commons.math4.util.Decimal64Field;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathArrays; import org.apache.commons.math4.util.MathArrays;
import org.junit.Assert; import org.junit.Assert;
import org.junit.Test; import org.junit.Test;
public class GillFieldIntegratorTest { public class GillFieldIntegratorTest extends AbstractRungeKuttaFieldIntegratorTest {
@Test(expected=DimensionMismatchException.class) protected <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T>
public void testDimensionCheck() createIntegrator(Field<T> field, T step) {
throws DimensionMismatchException, NumberIsTooSmallException, return new GillFieldIntegrator<T>(field, step);
MaxCountExceededException, NoBracketingException {
doTestDimensionCheck(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestDimensionCheck(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
new GillFieldIntegrator<T>(field, field.getZero().add(0.01)).integrate(new FieldExpandableODE<>(pb),
new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension() + 10)),
field.getOne());
Assert.fail("an exception should have been thrown");
} }
@Test @Test
public void testDecreasingSteps() public void testNonFieldIntegratorConsistency() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestNonFieldIntegratorConsistency(Decimal64Field.getInstance());
MaxCountExceededException, NoBracketingException {
doTestDecreasingSteps(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestDecreasingSteps(Field<T> field)
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 = 5; i < 10; ++i) {
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(FastMath.pow(2.0, -i));
FieldFirstOrderIntegrator<T> integ = new GillFieldIntegrator<T>(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(), 1.0e-10);
}
T error = handler.getMaximalValueError();
if (i > 5) {
Assert.assertTrue(error.subtract(previousValueError.abs().multiply(1.01)).getReal() < 0);
}
previousValueError = error;
T timeError = handler.getMaximalTimeError();
if (i > 5) {
Assert.assertTrue(timeError.subtract(previousTimeError.abs()).getReal() <= 0);
}
previousTimeError = timeError;
integ.clearEventHandlers();
Assert.assertEquals(0, integ.getEventHandlers().size());
}
}
} }
@Test @Test
public void testSmallStep() public void testMissedEndEvent() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestMissedEndEvent(Decimal64Field.getInstance(), 1.0e-15, 6.0e-5);
MaxCountExceededException, NoBracketingException {
doTestSmallStep(Decimal64Field.getInstance());
} }
private <T extends RealFieldElement<T>>void doTestSmallStep(Field<T> field) @Test
throws DimensionMismatchException, NumberIsTooSmallException, public void testSanityChecks() {
MaxCountExceededException, NoBracketingException { doTestSanityChecks(Decimal64Field.getInstance());
}
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field); @Test
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001); public void testDecreasingSteps() {
doTestDecreasingSteps(Decimal64Field.getInstance(), 1.0, 1.0, 1.0e-10);
}
FieldFirstOrderIntegrator<T> integ = new GillFieldIntegrator<T>(field, step); @Test
TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ); public void testSmallStep() {
integ.addStepHandler(handler); doTestSmallStep(Decimal64Field.getInstance(), 2.0e-13, 4.0e-12, 1.0e-12, "Gill");
integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime()); }
Assert.assertTrue(handler.getLastError().getReal() < 2.0e-13); @Test
Assert.assertTrue(handler.getMaximalValueError().getReal() < 4.0e-12); public void testBigStep() {
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12); doTestBigStep(Decimal64Field.getInstance(), 0.0004, 0.005, 1.0e-12, "Gill");
Assert.assertEquals("Gill", integ.getName());
} }
@Test
public void testBackward() {
doTestBackward(Decimal64Field.getInstance(), 5.0e-10, 7.0e-10, 1.0e-12, "Gill");
}
@Test
public void testKepler() {
doTestKepler(Decimal64Field.getInstance(), 1.72e-3, 1.0e-5);
}
@Test
public void testStepSize() {
doTestStepSize(Decimal64Field.getInstance(), 1.0e-12);
}
@Test
public void testSingleStep() {
doTestSingleStep(Decimal64Field.getInstance(), 0.21);
}
@Test
public void testTooLargeFirstStep() {
doTestTooLargeFirstStep(Decimal64Field.getInstance());
}
@Test @Test
public void testBigStep() public void testUnstableDerivative() {
throws DimensionMismatchException, NumberIsTooSmallException, doTestUnstableDerivative(Decimal64Field.getInstance(), 1.0e-12);
MaxCountExceededException, NoBracketingException {
doTestBigStep(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestBigStep(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2);
FieldFirstOrderIntegrator<T> integ = new GillFieldIntegrator<T>(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() > 0.0004);
Assert.assertTrue(handler.getMaximalValueError().getReal() > 0.005);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
}
@Test
public void testBackward()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestBackward(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestBackward(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs();
FieldFirstOrderIntegrator<T> integ = new GillFieldIntegrator<T>(field, step);
TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
integ.addStepHandler(handler);
integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
Assert.assertTrue(handler.getLastError().getReal() < 5.0e-10);
Assert.assertTrue(handler.getMaximalValueError().getReal() < 7.0e-10);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
Assert.assertEquals("Gill", integ.getName());
}
@Test
public void testKepler()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestKepler(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestKepler(Field<T> field)
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);
FieldFirstOrderIntegrator<T> integ = new GillFieldIntegrator<T>(field, step);
integ.addStepHandler(new KeplerHandler<T>(pb));
integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
}
@Test
public void testUnstableDerivative()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestUnstableDerivative(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestUnstableDerivative(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
final StepFieldProblem<T> stepProblem = new StepFieldProblem<T>(field,
field.getZero().add(0.0),
field.getZero().add(1.0),
field.getZero().add(2.0));
FieldFirstOrderIntegrator<T> integ = new GillFieldIntegrator<T>(field, field.getZero().add(0.3));
integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
FieldODEStateAndDerivative<T> result = integ.integrate(new FieldExpandableODE<>(stepProblem),
new FieldODEState<>(field.getZero(), MathArrays.buildArray(field, 1)),
field.getZero().add(10.0));
Assert.assertEquals(8.0, result.getState()[0].getReal(), 1.0e-12);
}
private static class KeplerHandler<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
public KeplerHandler(TestFieldProblem3<T> pb) {
this.pb = pb;
maxError = pb.getField().getZero();
}
public void init(FieldODEStateAndDerivative<T> state0, T t) {
maxError = pb.getField().getZero();
}
public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) {
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) {
// even with more than 1000 evaluations per period,
// Gill is not able to integrate such an eccentric
// orbit with a good accuracy
Assert.assertTrue(maxError.getReal() > 0.001);
}
}
private T maxError;
private TestFieldProblem3<T> pb;
}
@Test
public void testStepSize()
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
doTestStepSize(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>>void doTestStepSize(Field<T> field)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
final double step = 1.23456;
FirstOrderIntegrator integ = new GillIntegrator(step);
integ.addStepHandler(new StepHandler() {
public void handleStep(StepInterpolator interpolator, boolean isLast) {
if (! isLast) {
Assert.assertEquals(step,
interpolator.getCurrentTime() - interpolator.getPreviousTime(),
1.0e-12);
}
}
public void init(double t0, double[] y0, double t) {
}
});
integ.integrate(new FirstOrderDifferentialEquations() {
public void computeDerivatives(double t, double[] y, double[] dot) {
dot[0] = 1.0;
}
public int getDimension() {
return 1;
}
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
} }
} }