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

Added test and fixed step interpolator for field version of Gill method.

This commit is contained in:
Luc Maisonobe 2016-01-06 12:40:54 +01:00
parent c4a093c1f2
commit 1f14ff6ce8
4 changed files with 410 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegrator.java
View File

@ -74,7 +74,7 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>>
@Override
protected T[][] getA() {
final T two = getField().getOne().multiply(2);
final T two = getField().getZero().add(2);
final T sqrtTwo = two.sqrt();
final T[][] a = MathArrays.buildArray(getField(), 3, -1);
@ -94,13 +94,13 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>>
@Override
protected T[] getB() {
final T two = getField().getOne().multiply(2);
final T two = getField().getZero().add(2);
final T sqrtTwo = two.sqrt();
final T[] b = MathArrays.buildArray(getField(), 4);
b[0] = fraction(1, 6);
b[1] = two.subtract(sqrtTwo).divide(6);
b[2] = two.add(sqrtTwo).divide(6);
b[1] = sqrtTwo.subtract(2).divide(-6);
b[2] = sqrtTwo.add(2).divide(6);
b[3] = b[0];
return b;

26
src/main/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolator.java
View File

@ -21,7 +21,6 @@ import org.apache.commons.math4.Field;
import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.FastMath;
/**
* This class implements a step interpolator for the Gill fourth
@ -60,10 +59,10 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
extends RungeKuttaFieldStepInterpolator<T> {
/** First Gill coefficient. */
private static final double ONE_MINUS_INV_SQRT_2 = 1 - FastMath.sqrt(0.5);
private final T one_minus_inv_sqrt_2;
/** Second Gill coefficient. */
private static final double ONE_PLUS_INV_SQRT_2 = 1 + FastMath.sqrt(0.5);
private final T one_plus_inv_sqrt_2;
/** Simple constructor.
* @param field field to which the time and state vector elements belong
@ -73,6 +72,9 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
GillFieldStepInterpolator(final Field<T> field, final boolean forward,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
final T sqrt = field.getZero().add(0.5).sqrt();
one_minus_inv_sqrt_2 = field.getOne().subtract(sqrt);
one_plus_inv_sqrt_2 = field.getOne().add(sqrt);
}
/** Copy constructor.
@ -82,6 +84,8 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
*/
GillFieldStepInterpolator(final GillFieldStepInterpolator<T> interpolator) {
super(interpolator);
one_minus_inv_sqrt_2 = interpolator.one_minus_inv_sqrt_2;
one_plus_inv_sqrt_2 = interpolator.one_plus_inv_sqrt_2;
}
/** {@inheritDoc} */
@ -103,8 +107,8 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
final T fourTheta2 = twoTheta.multiply(twoTheta);
final T coeffDot1 = theta.multiply(twoTheta.subtract(3)).add(1);
final T cDot23 = twoTheta.multiply(one.subtract(theta));
final T coeffDot2 = cDot23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeffDot3 = cDot23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeffDot2 = cDot23.multiply(one_minus_inv_sqrt_2);
final T coeffDot3 = cDot23.multiply(one_plus_inv_sqrt_2);
final T coeffDot4 = theta.multiply(twoTheta.subtract(1));
final T[] interpolatedState;
final T[] interpolatedDerivatives;
@ -112,9 +116,9 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h).divide(6.0);
final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6));
final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
final T coeff2 = c23.multiply(one_minus_inv_sqrt_2);
final T coeff3 = c23.multiply(one_plus_inv_sqrt_2);
final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
@ -122,10 +126,10 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
final T s = oneMinusThetaH.divide(-6.0);
final T c23 = s.multiply(twoTheta.add(2).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(5)).add(1));
final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeff2 = c23.multiply(one_minus_inv_sqrt_2);
final T coeff3 = c23.multiply(one_plus_inv_sqrt_2);
final T coeff4 = s.multiply(fourTheta2.add(theta).add(1));
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
}

313
src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldIntegratorTest.java Normal file
View File

@ -0,0 +1,313 @@
/*
* 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.FieldFirstOrderIntegrator;
import org.apache.commons.math4.ode.FieldODEState;
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.FastMath;
import org.apache.commons.math4.util.MathArrays;
import org.junit.Assert;
import org.junit.Test;
public class GillFieldIntegratorTest {
@Test(expected=DimensionMismatchException.class)
public void testDimensionCheck()
throws DimensionMismatchException, NumberIsTooSmallException,
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
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 = 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
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 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() < 2.0e-13);
Assert.assertTrue(handler.getMaximalValueError().getReal() < 4.0e-12);
Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
Assert.assertEquals("Gill", 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 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]);
}
}

78
src/test/java/org/apache/commons/math4/ode/nonstiff/StepFieldProblem.java Normal file
View File

@ -0,0 +1,78 @@
/*
* 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.ode.FieldFirstOrderDifferentialEquations;
import org.apache.commons.math4.ode.FieldODEState;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.ode.events.Action;
import org.apache.commons.math4.ode.events.FieldEventHandler;
import org.apache.commons.math4.util.MathArrays;
public class StepFieldProblem<T extends RealFieldElement<T>>
implements FieldFirstOrderDifferentialEquations<T>, FieldEventHandler<T> {
public StepFieldProblem(Field<T> field, T rateBefore, T rateAfter, T switchTime) {
this.field = field;
this.rateAfter = rateAfter;
this.switchTime = switchTime;
setRate(rateBefore);
}
public T[] computeDerivatives(T t, T[] y) {
T[] yDot = MathArrays.buildArray(field, 1);
yDot[0] = rate;
return yDot;
}
public int getDimension() {
return 1;
}
public void setRate(T rate) {
this.rate = rate;
}
public void init(T t0, T[] y0, T t) {
}
public void init(FieldODEStateAndDerivative<T> state0, T t) {
}
public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
setRate(rateAfter);
return Action.RESET_DERIVATIVES;
}
public T g(FieldODEStateAndDerivative<T> state) {
return state.getTime().subtract(switchTime);
}
public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
return state;
}
private Field<T> field;
private T rate;
private T rateAfter;
private T switchTime;
}