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

Fixed Luther step interpolator error. 

Now all tests related to the Luther integrator and step interpolator
pass.
This commit is contained in:
Luc Maisonobe 2016-01-06 12:41:21 +01:00
parent b5c1893a70
commit c23335bae4
3 changed files with 26 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolator.java
View File

@ -88,13 +88,13 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
LutherFieldStepInterpolator(final Field<T> field, final boolean forward, LutherFieldStepInterpolator(final Field<T> field, final boolean forward,
final FieldEquationsMapper<T> mapper) { final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper); super(field, forward, mapper);
final T q = field.getOne().multiply(21).sqrt(); final T q = field.getZero().add(21).sqrt();
c5a = q.multiply( -49).add( -49); c5a = q.multiply( -49).add( -49);
c5b = q.multiply( 287).add( 392); c5b = q.multiply( 287).add( 392);
c5c = q.multiply( -357).add( -637); c5c = q.multiply( -357).add( -637);
c5d = q.multiply( 343).add( 833); c5d = q.multiply( 343).add( 833);
c6a = q.multiply( 49).add( -49); c6a = q.multiply( 49).add( -49);
c6b = q.multiply( -287).add( -392); c6b = q.multiply( -287).add( 392);
c6c = q.multiply( 357).add( -637); c6c = q.multiply( 357).add( -637);
c6d = q.multiply( -343).add( 833); c6d = q.multiply( -343).add( 833);
d5a = q.multiply( 49).add( 49); d5a = q.multiply( 49).add( 49);
@ -192,32 +192,32 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
final T coeffDot4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567 / 5.0).add( 972 / 5.0)).add( -486 / 5.0 )).add( 324 / 25.0)); final T coeffDot4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567 / 5.0).add( 972 / 5.0)).add( -486 / 5.0 )).add( 324 / 25.0));
final T coeffDot5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150))); final T coeffDot5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150)));
final T coeffDot6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150))); final T coeffDot6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150)));
final T coeffDot7 = theta.multiply(theta.multiply(theta.multiply( 3 )).add( -3 )).add( 3 / 5.0); final T coeffDot7 = theta.multiply(theta.multiply(theta.multiply( 3.0 ).add( -3 )).add( 3 / 5.0));
final T[] interpolatedState; final T[] interpolatedState;
final T[] interpolatedDerivatives; final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) { if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(theta.multiply(h)); final T s = theta.multiply(h);
final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1); final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1));
final T coeff2 = getField().getZero(); final T coeff2 = getField().getZero();
final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0)); final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0)));
final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0)); final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0)));
final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))); final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))));
final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300))); final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300))));
final T coeff7 = s.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0); final T coeff7 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 3 / 4.0 ).add( -1 )).add( 3 / 10.0)));
interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7); interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7); interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
} else { } else {
final T s = oneMinusThetaH.multiply(theta); final T s = oneMinusThetaH;
final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0); final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add(- 1 / 20.0));
final T coeff2 = getField().getZero(); final T coeff2 = getField().getZero();
final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0); final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0));
final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0))); final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0))));
final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0); final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0));
final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0); final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0));
final T coeff7 = s.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0); final T coeff7 = s.multiply( theta.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0));
interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7); interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7); interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
} }

16
src/test/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegratorTest.java
View File

@ -64,40 +64,40 @@ public class LutherFieldIntegratorTest extends AbstractRungeKuttaFieldIntegrator
public void testSmallStep() public void testSmallStep()
throws DimensionMismatchException, NumberIsTooSmallException, throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException { MaxCountExceededException, NoBracketingException {
doTestSmallStep(Decimal64Field.getInstance(), 9.0e-17, 4.0e-15, 1.0e-12, "Luther"); doTestSmallStep(Decimal64Field.getInstance(), 8.7e-17, 3.6e-15, 1.0e-12, "Luther");
} }
@Test @Test
public void testBigStep() public void testBigStep()
throws DimensionMismatchException, NumberIsTooSmallException, throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException { MaxCountExceededException, NoBracketingException {
doTestBigStep(Decimal64Field.getInstance(), 0.00002, 0.001, 1.0e-12, "Luther"); doTestBigStep(Decimal64Field.getInstance(), 2.7e-5, 1.7e-3, 1.0e-12, "Luther");
} }
@Test @Test
public void testBackward() public void testBackward()
throws DimensionMismatchException, NumberIsTooSmallException, throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException { MaxCountExceededException, NoBracketingException {
doTestBackward(Decimal64Field.getInstance(), 3.0e-13, 5.0e-13, 1.0e-12, "Luther"); doTestBackward(Decimal64Field.getInstance(), 2.4e-13, 4.3e-13, 1.0e-12, "Luther");
} }
@Test @Test
public void testKepler() public void testKepler()
throws DimensionMismatchException, NumberIsTooSmallException, throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException { MaxCountExceededException, NoBracketingException {
doTestKepler(Decimal64Field.getInstance(), 2.2e-7, 1.0e-8); doTestKepler(Decimal64Field.getInstance(), 2.18e-7, 4.0e-10);
} }
@Test @Test
public void testStepSize() public void testStepSize()
throws DimensionMismatchException, NumberIsTooSmallException, throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException { MaxCountExceededException, NoBracketingException {
doTestStepSize(Decimal64Field.getInstance(), 1.0e-12); doTestStepSize(Decimal64Field.getInstance(), 1.0e-22);
} }
@Test @Test
public void testSingleStep() { public void testSingleStep() {
doTestSingleStep(Decimal64Field.getInstance(), 1.0e-11); doTestSingleStep(Decimal64Field.getInstance(), 6.0e-12);
} }
@Test @Test
@ -107,12 +107,12 @@ public class LutherFieldIntegratorTest extends AbstractRungeKuttaFieldIntegrator
@Test @Test
public void testUnstableDerivative() { public void testUnstableDerivative() {
doTestUnstableDerivative(Decimal64Field.getInstance(), 1.0e-12); doTestUnstableDerivative(Decimal64Field.getInstance(), 4.0e-15);
} }
@Test @Test
public void testDerivativesConsistency() { public void testDerivativesConsistency() {
doTestDerivativesConsistency(Decimal64Field.getInstance(), 1.0e-10); doTestDerivativesConsistency(Decimal64Field.getInstance(), 1.0e-20);
} }
} }

2
src/test/java/org/apache/commons/math4/ode/nonstiff/LutherFieldStepInterpolatorTest.java
View File

@ -38,7 +38,7 @@ public class LutherFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStep
@Test @Test
public void interpolationInside() { public void interpolationInside() {
doInterpolationInside(Decimal64Field.getInstance(), 3.3e-14, 7.9e-13); doInterpolationInside(Decimal64Field.getInstance(), 1.1e-7, 9.6e-9);
} }
} }