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

Use immutable step interpolators.

This commit is contained in:
Luc Maisonobe 2016-01-06 13:22:39 +01:00
parent 0ddec2917a
commit 771eb6a606
37 changed files with 984 additions and 571 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/main/java/org/apache/commons/math4/ode/AbstractFieldIntegrator.java
View File

@ -307,6 +307,7 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp
}
}
AbstractFieldStepInterpolator<T> restricted = interpolator;
while (!occurringEvents.isEmpty()) {
// handle the chronologically first event
@ -315,11 +316,10 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp
iterator.remove();
// get state at event time
final FieldODEStateAndDerivative<T> eventState = interpolator.getInterpolatedState(currentEvent.getEventTime());
final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
// restrict the interpolator to the first part of the step, up to the event
interpolator.setSoftPreviousState(previousState);
interpolator.setSoftCurrentState(eventState);
restricted = restricted.restrictStep(previousState, eventState);
// advance all event states to current time
for (final FieldEventState<T> state : eventsStates) {
@ -329,7 +329,7 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp
// handle the first part of the step, up to the event
for (final FieldStepHandler<T> handler : stepHandlers) {
handler.handleStep(interpolator, isLastStep);
handler.handleStep(restricted, isLastStep);
}
if (isLastStep) {
@ -351,11 +351,10 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp
// prepare handling of the remaining part of the step
previousState = eventState;
interpolator.setSoftPreviousState(eventState);
interpolator.setSoftCurrentState(currentState);
restricted = restricted.restrictStep(eventState, currentState);
// check if the same event occurs again in the remaining part of the step
if (currentEvent.evaluateStep(interpolator)) {
if (currentEvent.evaluateStep(restricted)) {
// the event occurs during the current step
occurringEvents.add(currentEvent);
}
@ -371,7 +370,7 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp
// handle the remaining part of the step, after all events if any
for (FieldStepHandler<T> handler : stepHandlers) {
handler.handleStep(interpolator, isLastStep);
handler.handleStep(restricted, isLastStep);
}
return currentState;

4
src/main/java/org/apache/commons/math4/ode/ContinuousOutputFieldModel.java
View File

@ -156,7 +156,7 @@ public class ContinuousOutputFieldModel<T extends RealFieldElement<T>>
}
for (FieldStepInterpolator<T> interpolator : model.steps) {
steps.add(interpolator.copy());
steps.add(interpolator);
}
index = steps.size() - 1;
@ -201,7 +201,7 @@ public class ContinuousOutputFieldModel<T extends RealFieldElement<T>>
forward = interpolator.isForward();
}
steps.add(interpolator.copy());
steps.add(interpolator);
if (isLast) {
finalTime = interpolator.getCurrentState().getTime();

11
src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
@ -100,8 +101,14 @@ public class ClassicalRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected ClassicalRungeKuttaFieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldEquationsMapper<T> mapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
}

42
src/main/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
View File

@ -62,26 +62,36 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
ClassicalRungeKuttaFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
ClassicalRungeKuttaFieldStepInterpolator(final ClassicalRungeKuttaFieldStepInterpolator<T> interpolator) {
super(interpolator);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
}
/** {@inheritDoc} */
@Override
protected ClassicalRungeKuttaFieldStepInterpolator<T> doCopy() {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(this);
protected ClassicalRungeKuttaFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@ -89,9 +99,9 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
final T one = getField().getOne();
final T one = time.getField().getOne();
final T oneMinusTheta = one.subtract(theta);
final T oneMinus2Theta = one.subtract(theta.multiply(2));
final T coeffDot1 = oneMinusTheta.multiply(oneMinus2Theta);
@ -102,7 +112,7 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T fourTheta2 = theta.multiply(theta).multiply(4);
final T s = theta.multiply(h).divide(6.0);
final T s = thetaH.divide(6.0);
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
final T coeff23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));

10
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
import org.apache.commons.math4.util.MathUtils;
@ -189,8 +190,13 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected DormandPrince54FieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new DormandPrince54FieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
return new DormandPrince54FieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
/** {@inheritDoc} */

64
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
View File

@ -75,11 +75,23 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
DormandPrince54FieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
final T one = field.getOne();
a70 = one.multiply( 35.0).divide( 384.0);
a72 = one.multiply( 500.0).divide(1113.0);
@ -94,43 +106,27 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
d6 = one.multiply( 69997945.0).divide( 29380423.0);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
DormandPrince54FieldStepInterpolator(final DormandPrince54FieldStepInterpolator<T> interpolator) {
super(interpolator);
a70 = interpolator.a70;
a72 = interpolator.a72;
a73 = interpolator.a73;
a74 = interpolator.a74;
a75 = interpolator.a75;
d0 = interpolator.d0;
d2 = interpolator.d2;
d3 = interpolator.d3;
d4 = interpolator.d4;
d5 = interpolator.d5;
d6 = interpolator.d6;
}
/** {@inheritDoc} */
@Override
protected DormandPrince54FieldStepInterpolator<T> doCopy() {
return new DormandPrince54FieldStepInterpolator<T>(this);
protected DormandPrince54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new DormandPrince54FieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
// interpolate
final T one = getField().getOne();
final T one = time.getField().getOne();
final T eta = one.subtract(theta);
final T twoTheta = theta.multiply(2);
final T dot2 = one.subtract(twoTheta);
@ -139,7 +135,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T f1 = h.multiply(theta);
final T f1 = thetaH;
final T f2 = f1.multiply(eta);
final T f3 = f2.multiply(theta);
final T f4 = f3.multiply(eta);
@ -147,7 +143,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
subtract(f2.multiply(a70.subtract(1))).
add(f3.multiply(a70.multiply(2).subtract(1))).
add(f4.multiply(d0));
final T coeff1 = getField().getZero();
final T coeff1 = time.getField().getZero();
final T coeff2 = f1.multiply(a72).
subtract(f2.multiply(a72)).
add(f3.multiply(a72.multiply(2))).
@ -169,7 +165,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
subtract(dot2.multiply(a70.subtract(1))).
add(dot3.multiply(a70.multiply(2).subtract(1))).
add(dot4.multiply(d0));
final T coeffDot1 = getField().getZero();
final T coeffDot1 = time.getField().getZero();
final T coeffDot2 = a72.
subtract(dot2.multiply(a72)).
add(dot3.multiply(a72.multiply(2))).
@ -200,7 +196,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
subtract(f2.multiply(a70.subtract(1))).
add(f3.multiply(a70.multiply(2).subtract(1))).
add(f4.multiply(d0));
final T coeff1 = getField().getZero();
final T coeff1 = time.getField().getZero();
final T coeff2 = f1.multiply(a72).
subtract(f2.multiply(a72)).
add(f3.multiply(a72.multiply(2))).
@ -222,7 +218,7 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
subtract(dot2.multiply(a70.subtract(1))).
add(dot3.multiply(a70.multiply(2).subtract(1))).
add(dot4.multiply(d0));
final T coeffDot1 = getField().getZero();
final T coeffDot1 = time.getField().getZero();
final T coeffDot2 = a72.
subtract(dot2.multiply(a72)).
add(dot3.multiply(a72.multiply(2))).

10
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
import org.apache.commons.math4.util.MathUtils;
@ -395,8 +396,13 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected DormandPrince853FieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new DormandPrince853FieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
return new DormandPrince853FieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
/** {@inheritDoc} */

227
src/main/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
View File

@ -45,32 +45,43 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
DormandPrince853FieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
// interpolation weights
d = MathArrays.buildArray(getField(), 7, 16);
d = MathArrays.buildArray(field, 7, 16);
// this row is the same as the b array
d[0][ 0] = fraction(104257, 1920240);
d[0][ 1] = getField().getZero();
d[0][ 2] = getField().getZero();
d[0][ 3] = getField().getZero();
d[0][ 4] = getField().getZero();
d[0][ 5] = fraction( 3399327.0, 763840.0);
d[0][ 6] = fraction( 66578432.0, 35198415.0);
d[0][ 7] = fraction( -1674902723.0, 288716400.0);
d[0][ 8] = fraction( 54980371265625.0, 176692375811392.0);
d[0][ 9] = fraction( -734375.0, 4826304.0);
d[0][10] = fraction( 171414593.0, 851261400.0);
d[0][11] = fraction( 137909.0, 3084480.0);
d[0][12] = getField().getZero();
d[0][13] = getField().getZero();
d[0][14] = getField().getZero();
d[0][15] = getField().getZero();
d[0][ 0] = fraction(field, 104257, 1920240);
d[0][ 1] = field.getZero();
d[0][ 2] = field.getZero();
d[0][ 3] = field.getZero();
d[0][ 4] = field.getZero();
d[0][ 5] = fraction(field, 3399327.0, 763840.0);
d[0][ 6] = fraction(field, 66578432.0, 35198415.0);
d[0][ 7] = fraction(field, -1674902723.0, 288716400.0);
d[0][ 8] = fraction(field, 54980371265625.0, 176692375811392.0);
d[0][ 9] = fraction(field, -734375.0, 4826304.0);
d[0][10] = fraction(field, 171414593.0, 851261400.0);
d[0][11] = fraction(field, 137909.0, 3084480.0);
d[0][12] = field.getZero();
d[0][13] = field.getZero();
d[0][14] = field.getZero();
d[0][15] = field.getZero();
d[1][ 0] = d[0][ 0].negate().add(1);
d[1][ 1] = d[0][ 1].negate();
@ -106,105 +117,97 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
d[2][14] = d[0][14].multiply(2); // really 0
d[2][15] = d[0][15].multiply(2); // really 0
d[3][ 0] = fraction( -17751989329.0, 2106076560.0);
d[3][ 1] = getField().getZero();
d[3][ 2] = getField().getZero();
d[3][ 3] = getField().getZero();
d[3][ 4] = getField().getZero();
d[3][ 5] = fraction( 4272954039.0, 7539864640.0);
d[3][ 6] = fraction( -118476319744.0, 38604839385.0);
d[3][ 7] = fraction( 755123450731.0, 316657731600.0);
d[3][ 8] = fraction( 3692384461234828125.0, 1744130441634250432.0);
d[3][ 9] = fraction( -4612609375.0, 5293382976.0);
d[3][10] = fraction( 2091772278379.0, 933644586600.0);
d[3][11] = fraction( 2136624137.0, 3382989120.0);
d[3][12] = fraction( -126493.0, 1421424.0);
d[3][13] = fraction( 98350000.0, 5419179.0);
d[3][14] = fraction( -18878125.0, 2053168.0);
d[3][15] = fraction( -1944542619.0, 438351368.0);
d[3][ 0] = fraction(field, -17751989329.0, 2106076560.0);
d[3][ 1] = field.getZero();
d[3][ 2] = field.getZero();
d[3][ 3] = field.getZero();
d[3][ 4] = field.getZero();
d[3][ 5] = fraction(field, 4272954039.0, 7539864640.0);
d[3][ 6] = fraction(field, -118476319744.0, 38604839385.0);
d[3][ 7] = fraction(field, 755123450731.0, 316657731600.0);
d[3][ 8] = fraction(field, 3692384461234828125.0, 1744130441634250432.0);
d[3][ 9] = fraction(field, -4612609375.0, 5293382976.0);
d[3][10] = fraction(field, 2091772278379.0, 933644586600.0);
d[3][11] = fraction(field, 2136624137.0, 3382989120.0);
d[3][12] = fraction(field, -126493.0, 1421424.0);
d[3][13] = fraction(field, 98350000.0, 5419179.0);
d[3][14] = fraction(field, -18878125.0, 2053168.0);
d[3][15] = fraction(field, -1944542619.0, 438351368.0);
d[4][ 0] = fraction( 32941697297.0, 3159114840.0);
d[4][ 1] = getField().getZero();
d[4][ 2] = getField().getZero();
d[4][ 3] = getField().getZero();
d[4][ 4] = getField().getZero();
d[4][ 5] = fraction( 456696183123.0, 1884966160.0);
d[4][ 6] = fraction( 19132610714624.0, 115814518155.0);
d[4][ 7] = fraction( -177904688592943.0, 474986597400.0);
d[4][ 8] = fraction(-4821139941836765625.0, 218016305204281304.0);
d[4][ 9] = fraction( 30702015625.0, 3970037232.0);
d[4][10] = fraction( -85916079474274.0, 2800933759800.0);
d[4][11] = fraction( -5919468007.0, 634310460.0);
d[4][12] = fraction( 2479159.0, 157936.0);
d[4][13] = fraction( -18750000.0, 602131.0);
d[4][14] = fraction( -19203125.0, 2053168.0);
d[4][15] = fraction( 15700361463.0, 438351368.0);
d[4][ 0] = fraction(field, 32941697297.0, 3159114840.0);
d[4][ 1] = field.getZero();
d[4][ 2] = field.getZero();
d[4][ 3] = field.getZero();
d[4][ 4] = field.getZero();
d[4][ 5] = fraction(field, 456696183123.0, 1884966160.0);
d[4][ 6] = fraction(field, 19132610714624.0, 115814518155.0);
d[4][ 7] = fraction(field, -177904688592943.0, 474986597400.0);
d[4][ 8] = fraction(field, -4821139941836765625.0, 218016305204281304.0);
d[4][ 9] = fraction(field, 30702015625.0, 3970037232.0);
d[4][10] = fraction(field, -85916079474274.0, 2800933759800.0);
d[4][11] = fraction(field, -5919468007.0, 634310460.0);
d[4][12] = fraction(field, 2479159.0, 157936.0);
d[4][13] = fraction(field, -18750000.0, 602131.0);
d[4][14] = fraction(field, -19203125.0, 2053168.0);
d[4][15] = fraction(field, 15700361463.0, 438351368.0);
d[5][ 0] = fraction( 12627015655.0, 631822968.0);
d[5][ 1] = getField().getZero();
d[5][ 2] = getField().getZero();
d[5][ 3] = getField().getZero();
d[5][ 4] = getField().getZero();
d[5][ 5] = fraction( -72955222965.0, 188496616.0);
d[5][ 6] = fraction( -13145744952320.0, 69488710893.0);
d[5][ 7] = fraction( 30084216194513.0, 56998391688.0);
d[5][ 8] = fraction( -296858761006640625.0, 25648977082856624.0);
d[5][ 9] = fraction( 569140625.0, 82709109.0);
d[5][10] = fraction( -18684190637.0, 18672891732.0);
d[5][11] = fraction( 69644045.0, 89549712.0);
d[5][12] = fraction( -11847025.0, 4264272.0);
d[5][13] = fraction( -978650000.0, 16257537.0);
d[5][14] = fraction( 519371875.0, 6159504.0);
d[5][15] = fraction( 5256837225.0, 438351368.0);
d[5][ 0] = fraction(field, 12627015655.0, 631822968.0);
d[5][ 1] = field.getZero();
d[5][ 2] = field.getZero();
d[5][ 3] = field.getZero();
d[5][ 4] = field.getZero();
d[5][ 5] = fraction(field, -72955222965.0, 188496616.0);
d[5][ 6] = fraction(field, -13145744952320.0, 69488710893.0);
d[5][ 7] = fraction(field, 30084216194513.0, 56998391688.0);
d[5][ 8] = fraction(field, -296858761006640625.0, 25648977082856624.0);
d[5][ 9] = fraction(field, 569140625.0, 82709109.0);
d[5][10] = fraction(field, -18684190637.0, 18672891732.0);
d[5][11] = fraction(field, 69644045.0, 89549712.0);
d[5][12] = fraction(field, -11847025.0, 4264272.0);
d[5][13] = fraction(field, -978650000.0, 16257537.0);
d[5][14] = fraction(field, 519371875.0, 6159504.0);
d[5][15] = fraction(field, 5256837225.0, 438351368.0);
d[6][ 0] = fraction( -450944925.0, 17550638.0);
d[6][ 1] = getField().getZero();
d[6][ 2] = getField().getZero();
d[6][ 3] = getField().getZero();
d[6][ 4] = getField().getZero();
d[6][ 5] = fraction( -14532122925.0, 94248308.0);
d[6][ 6] = fraction( -595876966400.0, 2573655959.0);
d[6][ 7] = fraction( 188748653015.0, 527762886.0);
d[6][ 8] = fraction( 2545485458115234375.0, 27252038150535163.0);
d[6][ 9] = fraction( -1376953125.0, 36759604.0);
d[6][10] = fraction( 53995596795.0, 518691437.0);
d[6][11] = fraction( 210311225.0, 7047894.0);
d[6][12] = fraction( -1718875.0, 39484.0);
d[6][13] = fraction( 58000000.0, 602131.0);
d[6][14] = fraction( -1546875.0, 39484.0);
d[6][15] = fraction( -1262172375.0, 8429834.0);
d[6][ 0] = fraction(field, -450944925.0, 17550638.0);
d[6][ 1] = field.getZero();
d[6][ 2] = field.getZero();
d[6][ 3] = field.getZero();
d[6][ 4] = field.getZero();
d[6][ 5] = fraction(field, -14532122925.0, 94248308.0);
d[6][ 6] = fraction(field, -595876966400.0, 2573655959.0);
d[6][ 7] = fraction(field, 188748653015.0, 527762886.0);
d[6][ 8] = fraction(field, 2545485458115234375.0, 27252038150535163.0);
d[6][ 9] = fraction(field, -1376953125.0, 36759604.0);
d[6][10] = fraction(field, 53995596795.0, 518691437.0);
d[6][11] = fraction(field, 210311225.0, 7047894.0);
d[6][12] = fraction(field, -1718875.0, 39484.0);
d[6][13] = fraction(field, 58000000.0, 602131.0);
d[6][14] = fraction(field, -1546875.0, 39484.0);
d[6][15] = fraction(field, -1262172375.0, 8429834.0);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
DormandPrince853FieldStepInterpolator(final DormandPrince853FieldStepInterpolator<T> interpolator) {
super(interpolator);
d = MathArrays.buildArray(getField(), 4, -1);
for (int i = 0; i < d.length; ++i) {
d[i] = interpolator.d[i].clone();
}
/** {@inheritDoc} */
protected DormandPrince853FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new DormandPrince853FieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** Create a fraction.
* @param field field to which the elements belong
* @param p numerator
* @param q denominator
* @return p/q computed in the instance field
*/
private T fraction(final double p, final double q) {
return getField().getOne().multiply(p).divide(q);
}
/** {@inheritDoc} */
@Override
protected DormandPrince853FieldStepInterpolator<T> doCopy() {
return new DormandPrince853FieldStepInterpolator<T>(this);
private T fraction(final Field<T> field, final double p, final double q) {
return field.getZero().add(p).divide(q);
}
/** {@inheritDoc} */
@ -212,10 +215,10 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH)
final T thetaH, final T oneMinusThetaH)
throws MaxCountExceededException {
final T one = getField().getOne();
final T one = time.getField().getOne();
final T eta = one.subtract(theta);
final T twoTheta = theta.multiply(2);
final T theta2 = theta.multiply(theta);
@ -230,15 +233,15 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T f0 = theta.multiply(h);
final T f0 = thetaH;
final T f1 = f0.multiply(eta);
final T f2 = f1.multiply(theta);
final T f3 = f2.multiply(eta);
final T f4 = f3.multiply(theta);
final T f5 = f4.multiply(eta);
final T f6 = f5.multiply(theta);
final T[] p = MathArrays.buildArray(getField(), 16);
final T[] q = MathArrays.buildArray(getField(), 16);
final T[] p = MathArrays.buildArray(time.getField(), 16);
final T[] q = MathArrays.buildArray(time.getField(), 16);
for (int i = 0; i < p.length; ++i) {
p[i] = f0.multiply(d[0][i]).
add(f1.multiply(d[1][i])).
@ -267,8 +270,8 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
final T f4 = f3.multiply(theta);
final T f5 = f4.multiply(eta);
final T f6 = f5.multiply(theta);
final T[] p = MathArrays.buildArray(getField(), 16);
final T[] q = MathArrays.buildArray(getField(), 16);
final T[] p = MathArrays.buildArray(time.getField(), 16);
final T[] q = MathArrays.buildArray(time.getField(), 16);
for (int i = 0; i < p.length; ++i) {
p[i] = f0.multiply(d[0][i]).
add(f1.multiply(d[1][i])).

22
src/main/java/org/apache/commons/math4/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
View File

@ -183,10 +183,15 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
/** Create an interpolator.
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param mapper equations mapper for the all equations
* @return external weights for the high order method from Butcher array
*/
protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward,
protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
FieldEquationsMapper<T> mapper);
/** Get the order of the method.
* @return order of the method
@ -226,11 +231,6 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
final T[] yTmp = MathArrays.buildArray(getField(), y0.length);
// set up an interpolator sharing the integrator arrays
final RungeKuttaFieldStepInterpolator<T> interpolator =
createInterpolator(forward, equations.getMapper());
interpolator.storeState(stepStart);
// set up integration control objects
T hNew = getField().getZero();
boolean firstTime = true;
@ -239,8 +239,6 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
isLastStep = false;
do {
interpolator.shift();
// iterate over step size, ensuring local normalized error is smaller than 1
T error = getField().getZero().add(10);
while (error.subtract(1.0).getReal() >= 0) {
@ -314,16 +312,12 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp);
// local error is small enough: accept the step, trigger events and step handlers
interpolator.setSlopes(yDotK);
interpolator.storeState(stateTmp);
System.arraycopy(yTmp, 0, y, 0, y0.length);
stepStart = acceptStep(interpolator, finalTime);
stepStart = acceptStep(createInterpolator(forward, yDotK, stepStart, stateTmp, equations.getMapper()),
finalTime);
if (!isLastStep) {
// prepare next step
interpolator.storeState(stepStart);
// stepsize control for next step
final T factor = MathUtils.min(maxGrowth,
MathUtils.max(minReduction, safety.multiply(error.pow(exp))));

11
src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
@ -85,8 +86,14 @@ public class EulerFieldIntegrator<T extends RealFieldElement<T>> extends RungeKu
/** {@inheritDoc} */
@Override
protected EulerFieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new EulerFieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldEquationsMapper<T> mapper) {
return new EulerFieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
}

44
src/main/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolator.java
View File

@ -52,26 +52,36 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
EulerFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
EulerFieldStepInterpolator(final EulerFieldStepInterpolator<T> interpolator) {
super(interpolator);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
}
/** {@inheritDoc} */
@Override
protected EulerFieldStepInterpolator<T> doCopy() {
return new EulerFieldStepInterpolator<T>(this);
protected EulerFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new EulerFieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@ -79,15 +89,15 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>>
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
interpolatedState = previousStateLinearCombination(theta.multiply(h));
interpolatedDerivatives = derivativeLinearCombination(getField().getOne());
interpolatedState = previousStateLinearCombination(thetaH);
interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne());
} else {
interpolatedState = currentStateLinearCombination(oneMinusThetaH.negate());
interpolatedDerivatives = derivativeLinearCombination(getField().getOne());
interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne());
}
return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);

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

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
@ -110,8 +111,14 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected GillFieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new GillFieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldEquationsMapper<T> mapper) {
return new GillFieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
}

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

@ -67,42 +67,49 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
GillFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
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.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
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} */
@Override
protected GillFieldStepInterpolator<T> doCopy() {
return new GillFieldStepInterpolator<T>(this);
protected GillFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new GillFieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
final T one = getField().getOne();
final T one = time.getField().getOne();
final T twoTheta = theta.multiply(2);
final T fourTheta2 = twoTheta.multiply(twoTheta);
final T coeffDot1 = theta.multiply(twoTheta.subtract(3)).add(1);
@ -114,7 +121,7 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h).divide(6.0);
final T s = thetaH.divide(6.0);
final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
final T coeff2 = c23.multiply(one_minus_inv_sqrt_2);

10
src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
import org.apache.commons.math4.util.MathUtils;
@ -164,8 +165,13 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected HighamHall54FieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new HighamHall54FieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState, final FieldEquationsMapper<T> mapper) {
return new HighamHall54FieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
/** {@inheritDoc} */

57
src/main/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolator.java
View File

@ -38,38 +38,47 @@ class HighamHall54FieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
HighamHall54FieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
HighamHall54FieldStepInterpolator(final HighamHall54FieldStepInterpolator<T> interpolator) {
super(interpolator);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
}
/** {@inheritDoc} */
@Override
protected HighamHall54FieldStepInterpolator<T> doCopy() {
return new HighamHall54FieldStepInterpolator<T>(this);
protected HighamHall54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new HighamHall54FieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
final T bDot0 = theta.multiply(theta.multiply(theta.multiply( -10.0 ).add( 16.0 )).add(-15.0 / 2.0)).add(1);
final T bDot1 = getField().getZero();
final T bDot1 = time.getField().getZero();
final T bDot2 = theta.multiply(theta.multiply(theta.multiply( 135.0 / 2.0).add(-729.0 / 8.0)).add(459.0 / 16.0));
final T bDot3 = theta.multiply(theta.multiply(theta.multiply(-120.0 ).add( 152.0 )).add(-44.0 ));
final T bDot4 = theta.multiply(theta.multiply(theta.multiply( 125.0 / 2.0).add(-625.0 / 8.0)).add(375.0 / 16.0));
@ -78,19 +87,19 @@ class HighamHall54FieldStepInterpolator<T extends RealFieldElement<T>>
final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T hTheta = h.multiply(theta);
final T b0 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply( -5.0 / 2.0).add( 16.0 / 3.0)).add(-15.0 / 4.0)).add(1));
final T b1 = getField().getZero();
final T b2 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(135.0 / 8.0).add(-243.0 / 8.0)).add(459.0 / 32.0)));
final T b3 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(-30.0 ).add( 152.0 / 3.0)).add(-22.0 )));
final T b4 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(125.0 / 8.0).add(-625.0 / 24.0)).add(375.0 / 32.0)));
final T b5 = hTheta.multiply(theta.multiply(theta.multiply( 5.0 / 12.0 ).add( -5.0 / 16.0)));
final T b0 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply( -5.0 / 2.0).add( 16.0 / 3.0)).add(-15.0 / 4.0)).add(1));
final T b1 = time.getField().getZero();
final T b2 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply(135.0 / 8.0).add(-243.0 / 8.0)).add(459.0 / 32.0)));
final T b3 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply(-30.0 ).add( 152.0 / 3.0)).add(-22.0 )));
final T b4 = thetaH.multiply(theta.multiply(theta.multiply(theta.multiply(125.0 / 8.0).add(-625.0 / 24.0)).add(375.0 / 32.0)));
final T b5 = thetaH.multiply(theta.multiply(theta.multiply( 5.0 / 12.0 ).add( -5.0 / 16.0)));
interpolatedState = previousStateLinearCombination(b0, b1, b2, b3, b4, b5);
interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5);
} else {
final T theta2 = theta.multiply(theta);
final T h = thetaH.divide(theta);
final T b0 = h.multiply( theta.multiply(theta.multiply(theta.multiply(theta.multiply(-5.0 / 2.0).add( 16.0 / 3.0)).add( -15.0 / 4.0)).add( 1.0 )).add( -1.0 / 12.0));
final T b1 = getField().getZero();
final T b1 = time.getField().getZero();
final T b2 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 135.0 / 8.0 ).add(-243.0 / 8.0)).add(459.0 / 32.0)).add( -27.0 / 32.0));
final T b3 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( -30.0 ).add( 152.0 / 3.0)).add(-22.0 )).add( 4.0 / 3.0));
final T b4 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 125.0 / 8.0 ).add(-625.0 / 24.0)).add(375.0 / 32.0)).add(-125.0 / 96.0));

11
src/main/java/org/apache/commons/math4/ode/nonstiff/LutherFieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
@ -135,8 +136,14 @@ public class LutherFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected LutherFieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new LutherFieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldEquationsMapper<T> mapper) {
return new LutherFieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
}

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

@ -83,11 +83,23 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
LutherFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
final T q = field.getZero().add(21).sqrt();
c5a = q.multiply( -49).add( -49);
c5b = q.multiply( 287).add( 392);
@ -103,45 +115,27 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
d6a = q.multiply( -49).add( 49);
d6b = q.multiply( 847).add(-1372);
d6c = q.multiply(-1029).add( 2254);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
LutherFieldStepInterpolator(final LutherFieldStepInterpolator<T> interpolator) {
super(interpolator);
c5a = interpolator.c5a;
c5b = interpolator.c5b;
c5c = interpolator.c5c;
c5d = interpolator.c5d;
c6a = interpolator.c6a;
c6b = interpolator.c6b;
c6c = interpolator.c6c;
c6d = interpolator.c6d;
d5a = interpolator.d5a;
d5b = interpolator.d5b;
d5c = interpolator.d5c;
d6a = interpolator.d6a;
d6b = interpolator.d6b;
d6c = interpolator.d6c;
}
/** {@inheritDoc} */
@Override
protected LutherFieldStepInterpolator<T> doCopy() {
return new LutherFieldStepInterpolator<T>(this);
protected LutherFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new LutherFieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
// the coefficients below have been computed by solving the
// order conditions from a theorem from Butcher (1963), using
@ -187,7 +181,7 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
// At the end, we get the b_i as polynomials in theta.
final T coeffDot1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 ).add( -47 )).add( 36 )).add( -54 / 5.0)).add(1);
final T coeffDot2 = getField().getZero();
final T coeffDot2 = time.getField().getZero();
final T coeffDot3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 ).add(-608 / 3.0)).add( 320 / 3.0 )).add(-208 / 15.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)));
@ -198,9 +192,9 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h);
final T s = thetaH;
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 = time.getField().getZero();
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(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(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))));
@ -212,7 +206,7 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
final T s = oneMinusThetaH;
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 = time.getField().getZero();
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(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.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));

11
src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
@ -85,8 +86,14 @@ public class MidpointFieldIntegrator<T extends RealFieldElement<T>> extends Rung
/** {@inheritDoc} */
@Override
protected MidpointFieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new MidpointFieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldEquationsMapper<T> mapper) {
return new MidpointFieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
}

45
src/main/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolator.java
View File

@ -54,44 +54,53 @@ class MidpointFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
MidpointFieldStepInterpolator(Field<T> field, final boolean forward,
MidpointFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
MidpointFieldStepInterpolator(final MidpointFieldStepInterpolator<T> interpolator) {
super(interpolator);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
}
/** {@inheritDoc} */
@Override
protected MidpointFieldStepInterpolator<T> doCopy() {
return new MidpointFieldStepInterpolator<T>(this);
protected MidpointFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new MidpointFieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
final T coeffDot2 = theta.multiply(2);
final T coeffDot1 = getField().getOne().subtract(coeffDot2);
final T coeffDot1 = time.getField().getOne().subtract(coeffDot2);
final T[] interpolatedState;
final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T coeff1 = theta.multiply(oneMinusThetaH);
final T coeff2 = theta.multiply(theta).multiply(h);
final T coeff2 = theta.multiply(thetaH);
interpolatedState = previousStateLinearCombination(coeff1, coeff2);
interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2);
} else {

22
src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldIntegrator.java
View File

@ -99,10 +99,15 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
/** Create an interpolator.
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param mapper equations mapper for the all equations
* @return external weights for the high order method from Butcher array
*/
protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward,
protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
FieldEquationsMapper<T> mapper);
/** {@inheritDoc} */
@ -124,11 +129,6 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
final T[] yTmp = MathArrays.buildArray(getField(), y0.length);
// set up an interpolator sharing the integrator arrays
final RungeKuttaFieldStepInterpolator<T> interpolator =
createInterpolator(forward, equations.getMapper());
interpolator.storeState(stepStart);
// set up integration control objects
if (forward) {
if (stepStart.getTime().add(step).subtract(finalTime).getReal() >= 0) {
@ -148,8 +148,6 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
isLastStep = false;
do {
interpolator.shift();
// first stage
y = equations.getMapper().mapState(stepStart);
yDotK[0] = equations.getMapper().mapDerivative(stepStart);
@ -182,16 +180,12 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp);
// discrete events handling
interpolator.setSlopes(yDotK);
interpolator.storeState(stateTmp);
System.arraycopy(yTmp, 0, y, 0, y0.length);
stepStart = acceptStep(interpolator, finalTime);
stepStart = acceptStep(createInterpolator(forward, yDotK, stepStart, stateTmp, equations.getMapper()),
finalTime);
if (!isLastStep) {
// prepare next step
interpolator.storeState(stepStart);
// stepsize control for next step
final T nextT = stepStart.getTime().add(stepSize);
final boolean nextIsLast = forward ?

83
src/main/java/org/apache/commons/math4/ode/nonstiff/RungeKuttaFieldStepInterpolator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.ode.sampling.AbstractFieldStepInterpolator;
import org.apache.commons.math4.util.MathArrays;
@ -40,57 +41,63 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
private final Field<T> field;
/** Slopes at the intermediate points. */
private T[][] yDotK;
private final T[][] yDotK;
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
protected RungeKuttaFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(forward, mapper);
super(forward, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
this.field = field;
this.yDotK = null;
this.yDotK = MathArrays.buildArray(field, yDotK.length, -1);
for (int i = 0; i < yDotK.length; ++i) {
this.yDotK[i] = yDotK[i].clone();
}
}
/** Copy constructor.
* <p>The copy is a deep copy: its arrays are separated from the
* original arrays of the instance.</p>
* @param interpolator interpolator to copy from.
/** {@inheritDoc} */
protected RungeKuttaFieldStepInterpolator<T> create(boolean newForward,
FieldODEStateAndDerivative<T> newGlobalPreviousState,
FieldODEStateAndDerivative<T> newGlobalCurrentState,
FieldODEStateAndDerivative<T> newSoftPreviousState,
FieldODEStateAndDerivative<T> newSoftCurrentState,
FieldEquationsMapper<T> newMapper) {
return create(field, newForward, yDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** Create a new instance.
* @param newField field to which the time and state vector elements belong
* @param newForward integration direction indicator
* @param newYDotK slopes at the intermediate points
* @param newGlobalPreviousState start of the global step
* @param newGlobalCurrentState end of the global step
* @param newSoftPreviousState start of the restricted step
* @param newSoftCurrentState end of the restricted step
* @param newMapper equations mapper for the all equations
* @return a new instance
*/
RungeKuttaFieldStepInterpolator(final RungeKuttaFieldStepInterpolator<T> interpolator) {
super(interpolator);
field = interpolator.field;
if (yDotK != null) {
yDotK = MathArrays.buildArray(field, interpolator.yDotK.length, -1);
for (int k = 0; k < yDotK.length; ++k) {
yDotK[k] = interpolator.yDotK[k].clone();
}
} else {
yDotK = null;
}
}
/** Get the field to which the time and state vector elements belong.
* @return to which the time and state vector elements belong
*/
protected Field<T> getField() {
return field;
}
/** Store the slopes at the intermediate points.
* @param slopes slopes at the intermediate points
*/
void setSlopes(final T[][] slopes) {
this.yDotK = slopes.clone();
}
protected abstract RungeKuttaFieldStepInterpolator<T> create(Field<T> newField, boolean newForward, T[][] newYDotK,
FieldODEStateAndDerivative<T> newGlobalPreviousState,
FieldODEStateAndDerivative<T> newGlobalCurrentState,
FieldODEStateAndDerivative<T> newSoftPreviousState,
FieldODEStateAndDerivative<T> newSoftCurrentState,
FieldEquationsMapper<T> newMapper);
/** Compute a state by linear combination added to previous state.
* @param coefficients coefficients to apply to the method staged derivatives

11
src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldIntegrator.java
View File

@ -20,6 +20,7 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.MathArrays;
/**
@ -99,8 +100,14 @@ public class ThreeEighthesFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected ThreeEighthesFieldStepInterpolator<T>
createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
return new ThreeEighthesFieldStepInterpolator<T>(getField(), forward, mapper);
createInterpolator(final boolean forward, T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldEquationsMapper<T> mapper) {
return new ThreeEighthesFieldStepInterpolator<T>(getField(), forward, yDotK,
globalPreviousState, globalCurrentState,
globalPreviousState, globalCurrentState,
mapper);
}
}

41
src/main/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java
View File

@ -64,35 +64,44 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param field field to which the time and state vector elements belong
* @param forward integration direction indicator
* @param yDotK slopes at the intermediate points
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param mapper equations mapper for the all equations
*/
ThreeEighthesFieldStepInterpolator(final Field<T> field, final boolean forward,
final T[][] yDotK,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(field, forward, mapper);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
ThreeEighthesFieldStepInterpolator(final ThreeEighthesFieldStepInterpolator<T> interpolator) {
super(interpolator);
super(field, forward, yDotK,
globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
mapper);
}
/** {@inheritDoc} */
@Override
protected ThreeEighthesFieldStepInterpolator<T> doCopy() {
return new ThreeEighthesFieldStepInterpolator<T>(this);
protected ThreeEighthesFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new ThreeEighthesFieldStepInterpolator<T>(newField, newForward, newYDotK,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
/** {@inheritDoc} */
@SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
final T thetaH, final T oneMinusThetaH) {
final T coeffDot3 = theta.multiply(0.75);
final T coeffDot1 = coeffDot3.multiply(theta.multiply(4).subtract(5)).add(1);
@ -102,7 +111,7 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
final T[] interpolatedDerivatives;
if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h).divide(8);
final T s = thetaH.divide(8);
final T fourTheta2 = theta.multiply(theta).multiply(4);
final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(15)).add(8));
final T coeff2 = s.multiply(theta.multiply(5).subtract(fourTheta2)).multiply(3);

141
src/main/java/org/apache/commons/math4/ode/sampling/AbstractFieldStepInterpolator.java
View File

@ -40,118 +40,76 @@ import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T>>
implements FieldStepInterpolator<T> {
/** Current time step. */
protected T h;
/** Global previous state. */
private FieldODEStateAndDerivative<T> globalPreviousState;
private final FieldODEStateAndDerivative<T> globalPreviousState;
/** Global current state. */
private FieldODEStateAndDerivative<T> globalCurrentState;
private final FieldODEStateAndDerivative<T> globalCurrentState;
/** Soft previous state. */
private FieldODEStateAndDerivative<T> softPreviousState;
private final FieldODEStateAndDerivative<T> softPreviousState;
/** Soft current state. */
private FieldODEStateAndDerivative<T> softCurrentState;
private final FieldODEStateAndDerivative<T> softCurrentState;
/** integration direction. */
private boolean forward;
private final boolean forward;
/** Mapper for ODE equations primary and secondary components. */
private FieldEquationsMapper<T> mapper;
/** Simple constructor.
* @param isForward integration direction indicator
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param softPreviousState start of the restricted step
* @param softCurrentState end of the restricted step
* @param equationsMapper mapper for ODE equations primary and secondary components
*/
protected AbstractFieldStepInterpolator(final boolean isForward,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> equationsMapper) {
globalPreviousState = null;
globalCurrentState = null;
softPreviousState = null;
softCurrentState = null;
h = null;
this.forward = isForward;
this.globalPreviousState = globalPreviousState;
this.globalCurrentState = globalCurrentState;
this.softPreviousState = softPreviousState;
this.softCurrentState = softCurrentState;
this.mapper = equationsMapper;
}
/** Copy constructor.
* <p>The copy is a deep copy: its arrays are separated from the
* original arrays of the instance.</p>
* @param interpolator interpolator to copy from.
*/
protected AbstractFieldStepInterpolator(final AbstractFieldStepInterpolator<T> interpolator) {
globalPreviousState = interpolator.globalPreviousState;
globalCurrentState = interpolator.globalCurrentState;
softPreviousState = interpolator.softPreviousState;
softCurrentState = interpolator.softCurrentState;
h = interpolator.h;
forward = interpolator.forward;
mapper = interpolator.mapper;
}
/** {@inheritDoc} */
public FieldStepInterpolator<T> copy() throws MaxCountExceededException {
// create the new independent instance
return doCopy();
}
/** Really copy the instance.
* @return a copy of the instance
*/
protected abstract FieldStepInterpolator<T> doCopy();
/** Shift one step forward.
* Copy the current time into the previous time, hence preparing the
* interpolator for future calls to {@link #storeTime storeTime}
*/
public void shift() {
globalPreviousState = globalCurrentState;
softPreviousState = globalPreviousState;
softCurrentState = globalCurrentState;
}
/** Store the current step state.
* @param state current state
*/
public void storeState(final FieldODEStateAndDerivative<T> state) {
globalCurrentState = state;
softCurrentState = globalCurrentState;
if (globalPreviousState != null) {
h = globalCurrentState.getTime().subtract(globalPreviousState.getTime());
}
}
/** Restrict step range to a limited part of the global step.
/** Create a new restricted version of the instance.
* <p>
* This method can be used to restrict a step and make it appear
* as if the original step was smaller. Calling this method
* <em>only</em> changes the value returned by {@link #getPreviousState()},
* it does not change any other property
* The instance is not changed at all.
* </p>
* @param softPreviousState start of the restricted step
* @param previousState start of the restricted step
* @param currentState end of the restricted step
* @return restricted version of the instance
* @see #getPreviousState()
* @see #getCurrentState()
*/
public void setSoftPreviousState(final FieldODEStateAndDerivative<T> softPreviousState) {
this.softPreviousState = softPreviousState;
public AbstractFieldStepInterpolator<T> restrictStep(final FieldODEStateAndDerivative<T> previousState,
final FieldODEStateAndDerivative<T> currentState) {
return create(forward, globalPreviousState, globalCurrentState, previousState, currentState, mapper);
}
/** Restrict step range to a limited part of the global step.
* <p>
* This method can be used to restrict a step and make it appear
* as if the original step was smaller. Calling this method
* <em>only</em> changes the value returned by {@link #getCurrentState()},
* it does not change any other property
* </p>
* @param softCurrentState end of the restricted step
/** Create a new instance.
* @param newForward integration direction indicator
* @param newGlobalPreviousState start of the global step
* @param newGlobalCurrentState end of the global step
* @param newSoftPreviousState start of the restricted step
* @param newSoftCurrentState end of the restricted step
* @param newMapper equations mapper for the all equations
* @return a new instance
*/
public void setSoftCurrentState(final FieldODEStateAndDerivative<T> softCurrentState) {
this.softCurrentState = softCurrentState;
}
protected abstract AbstractFieldStepInterpolator<T> create(boolean newForward,
FieldODEStateAndDerivative<T> newGlobalPreviousState,
FieldODEStateAndDerivative<T> newGlobalCurrentState,
FieldODEStateAndDerivative<T> newSoftPreviousState,
FieldODEStateAndDerivative<T> newSoftCurrentState,
FieldEquationsMapper<T> newMapper);
/**
* Get the previous global grid point state.
@ -169,25 +127,22 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
return globalCurrentState;
}
/** {@inheritDoc}
* @see #setSoftPreviousState(FieldODEStateAndDerivative)
*/
/** {@inheritDoc} */
public FieldODEStateAndDerivative<T> getPreviousState() {
return softPreviousState;
}
/** {@inheritDoc}
* @see #setSoftCurrentState(FieldODEStateAndDerivative)
*/
/** {@inheritDoc} */
public FieldODEStateAndDerivative<T> getCurrentState() {
return softCurrentState;
}
/** {@inheritDoc} */
public FieldODEStateAndDerivative<T> getInterpolatedState(final T time) {
final T thetaH = time.subtract(globalPreviousState.getTime());
final T oneMinusThetaH = globalCurrentState.getTime().subtract(time);
final T theta = (h.getReal() == 0) ? h.getField().getZero() : h.subtract(oneMinusThetaH).divide(h);
return computeInterpolatedStateAndDerivatives(mapper, time, theta, oneMinusThetaH);
final T theta = thetaH.divide(globalCurrentState.getTime().subtract(globalPreviousState.getTime()));
return computeInterpolatedStateAndDerivatives(mapper, time, theta, thetaH, oneMinusThetaH);
}
/** {@inheritDoc} */
@ -202,13 +157,15 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
* @param time interpolation time
* @param theta normalized interpolation abscissa within the step
* (theta is zero at the previous time step and one at the current time step)
* @param thetaH time gap between the previous time and the interpolated time
* @param oneMinusThetaH time gap between the interpolated time and
* the current time
* @return interpolated state and derivatives
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
*/
protected abstract FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(FieldEquationsMapper<T> equationsMapper,
T time, T theta, T oneMinusThetaH)
T time, T theta,
T thetaH, T oneMinusThetaH)
throws MaxCountExceededException;
}

11
src/main/java/org/apache/commons/math4/ode/sampling/FieldStepInterpolator.java
View File

@ -18,7 +18,6 @@
package org.apache.commons.math4.ode.sampling;
import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.exception.MaxCountExceededException;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
/** This interface represents an interpolator over the last step
@ -74,14 +73,4 @@ public interface FieldStepInterpolator<T extends RealFieldElement<T>> {
*/
boolean isForward();
/** Copy the instance.
* <p>The copied instance is guaranteed to be independent from the
* original one. Both can be used with different settings for
* interpolated time without any side effect.</p>
* @return a deep copy of the instance, which can be used independently.
* @exception MaxCountExceededException if the number of functions evaluations is exceeded
* during step finalization
*/
FieldStepInterpolator<T> copy() throws MaxCountExceededException;
}

219
src/test/java/org/apache/commons/math4/ode/ContinuousOutputFieldModelTest.java Normal file
View File

@ -0,0 +1,219 @@
/*
* 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;
import java.util.Random;
import org.apache.commons.math4.Field;
import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.nonstiff.DormandPrince54FieldIntegrator;
import org.apache.commons.math4.ode.nonstiff.DormandPrince853FieldIntegrator;
import org.apache.commons.math4.ode.sampling.DummyFieldStepInterpolator;
import org.apache.commons.math4.ode.sampling.FieldStepInterpolator;
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.MathUtils;
import org.junit.Assert;
import org.junit.Test;
public class ContinuousOutputFieldModelTest {
@Test
public void testBoundaries() {
doTestBoundaries(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestBoundaries(final Field<T> field) {
TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9));
double minStep = 0;
double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
FieldFirstOrderIntegrator<T> integ = new DormandPrince54FieldIntegrator<T>(field, minStep, maxStep, 1.0e-8, 1.0e-8);
integ.addStepHandler(new ContinuousOutputFieldModel<T>());
integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
ContinuousOutputFieldModel<T> cm = (ContinuousOutputFieldModel<T>) integ.getStepHandlers().iterator().next();
cm.getInterpolatedState(pb.getInitialState().getTime().multiply(2).subtract(pb.getFinalTime()));
cm.getInterpolatedState(pb.getFinalTime().multiply(2).subtract(pb.getInitialState().getTime()));
cm.getInterpolatedState(pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5));
}
@Test
public void testRandomAccess() {
doTestRandomAccess(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestRandomAccess(final Field<T> field) {
TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9));
double minStep = 0;
double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
FieldFirstOrderIntegrator<T> integ = new DormandPrince54FieldIntegrator<T>(field, minStep, maxStep, 1.0e-8, 1.0e-8);
ContinuousOutputFieldModel<T> cm = new ContinuousOutputFieldModel<T>();
integ.addStepHandler(cm);
integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
Random random = new Random(347588535632l);
T maxError = field.getZero();
T maxErrorDot = field.getZero();
for (int i = 0; i < 1000; ++i) {
double r = random.nextDouble();
T time = pb.getInitialState().getTime().multiply(r).add(pb.getFinalTime().multiply(1.0 - r));
FieldODEStateAndDerivative<T> interpolated = cm.getInterpolatedState(time);
T[] theoreticalY = pb.computeTheoreticalState(time);
T[] theoreticalYDot = pb.doComputeDerivatives(time, theoreticalY);
T dx = interpolated.getState()[0].subtract(theoreticalY[0]);
T dy = interpolated.getState()[1].subtract(theoreticalY[1]);
T error = dx.multiply(dx).add(dy.multiply(dy));
maxError = MathUtils.max(maxError, error);
T dxDot = interpolated.getDerivative()[0].subtract(theoreticalYDot[0]);
T dyDot = interpolated.getDerivative()[1].subtract(theoreticalYDot[1]);
T errorDot = dxDot.multiply(dxDot).add(dyDot.multiply(dyDot));
maxErrorDot = MathUtils.max(maxErrorDot, errorDot);
}
Assert.assertEquals(0.0, maxError.getReal(), 1.0e-9);
Assert.assertEquals(0.0, maxErrorDot.getReal(), 4.0e-7);
}
@Test
public void testModelsMerging() {
doTestModelsMerging(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestModelsMerging(final Field<T> field) {
// theoretical solution: y[0] = cos(t), y[1] = sin(t)
FieldFirstOrderDifferentialEquations<T> problem =
new FieldFirstOrderDifferentialEquations<T>() {
public T[] computeDerivatives(T t, T[] y) {
T[] yDot = MathArrays.buildArray(field, 2);
yDot[0] = y[1].negate();
yDot[1] = y[0];
return yDot;
}
public int getDimension() {
return 2;
}
public void init(T t0, T[] y0, T finalTime) {
}
};
// integrate backward from &pi; to 0;
ContinuousOutputFieldModel<T> cm1 = new ContinuousOutputFieldModel<T>();
FieldFirstOrderIntegrator<T> integ1 =
new DormandPrince853FieldIntegrator<T>(field, 0, 1.0, 1.0e-8, 1.0e-8);
integ1.addStepHandler(cm1);
T t0 = field.getZero().add(FastMath.PI);
T[] y0 = MathArrays.buildArray(field, 2);
y0[0] = field.getOne().negate();
y0[1] = field.getZero();
integ1.integrate(new FieldExpandableODE<T>(problem),
new FieldODEState<T>(t0, y0),
field.getZero());
// integrate backward from 2&pi; to &pi;
ContinuousOutputFieldModel<T> cm2 = new ContinuousOutputFieldModel<T>();
FieldFirstOrderIntegrator<T> integ2 =
new DormandPrince853FieldIntegrator<T>(field, 0, 0.1, 1.0e-12, 1.0e-12);
integ2.addStepHandler(cm2);
t0 = field.getZero().add(2.0 * FastMath.PI);
y0[0] = field.getOne();
y0[1] = field.getZero();
integ2.integrate(new FieldExpandableODE<T>(problem),
new FieldODEState<T>(t0, y0),
field.getZero().add(FastMath.PI));
// merge the two half circles
ContinuousOutputFieldModel<T> cm = new ContinuousOutputFieldModel<T>();
cm.append(cm2);
cm.append(new ContinuousOutputFieldModel<T>());
cm.append(cm1);
// check circle
Assert.assertEquals(2.0 * FastMath.PI, cm.getInitialTime().getReal(), 1.0e-12);
Assert.assertEquals(0, cm.getFinalTime().getReal(), 1.0e-12);
for (double t = 0; t < 2.0 * FastMath.PI; t += 0.1) {
FieldODEStateAndDerivative<T> interpolated = cm.getInterpolatedState(field.getZero().add(t));
Assert.assertEquals(FastMath.cos(t), interpolated.getState()[0].getReal(), 1.0e-7);
Assert.assertEquals(FastMath.sin(t), interpolated.getState()[1].getReal(), 1.0e-7);
}
}
@Test
public void testErrorConditions() {
doTestErrorConditions(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestErrorConditions(final Field<T> field) {
ContinuousOutputFieldModel<T> cm = new ContinuousOutputFieldModel<T>();
cm.handleStep(buildInterpolator(field, 0, 1, new double[] { 0.0, 1.0, -2.0 }), true);
// dimension mismatch
Assert.assertTrue(checkAppendError(field, cm, 1.0, 2.0, new double[] { 0.0, 1.0 }));
// hole between time ranges
Assert.assertTrue(checkAppendError(field, cm, 10.0, 20.0, new double[] { 0.0, 1.0, -2.0 }));
// propagation direction mismatch
Assert.assertTrue(checkAppendError(field, cm, 1.0, 0.0, new double[] { 0.0, 1.0, -2.0 }));
// no errors
Assert.assertFalse(checkAppendError(field, cm, 1.0, 2.0, new double[] { 0.0, 1.0, -2.0 }));
}
private <T extends RealFieldElement<T>> boolean checkAppendError(Field<T> field, ContinuousOutputFieldModel<T> cm,
double t0, double t1, double[] y) {
try {
ContinuousOutputFieldModel<T> otherCm = new ContinuousOutputFieldModel<T>();
otherCm.handleStep(buildInterpolator(field, t0, t1, y), true);
cm.append(otherCm);
} catch(IllegalArgumentException iae) {
return true; // there was an allowable error
}
return false; // no allowable error
}
private <T extends RealFieldElement<T>> FieldStepInterpolator<T> buildInterpolator(Field<T> field,
double t0, double t1, double[] y) {
T[] fieldY = MathArrays.buildArray(field, y.length);
for (int i = 0; i < y.length; ++i) {
fieldY[i] = field.getZero().add(y[i]);
}
final FieldODEStateAndDerivative<T> s0 = new FieldODEStateAndDerivative<T>(field.getZero().add(t0), fieldY, fieldY);
final FieldODEStateAndDerivative<T> s1 = new FieldODEStateAndDerivative<T>(field.getZero().add(t1), fieldY, fieldY);
final FieldEquationsMapper<T> mapper = new FieldExpandableODE<T>(new FieldFirstOrderDifferentialEquations<T>() {
public int getDimension() {
return s0.getStateDimension();
}
public void init(T t0, T[] y0, T finalTime) {
}
public T[] computeDerivatives(T t, T[] y) {
return y;
}
}).getMapper();
return new DummyFieldStepInterpolator<T>(t1 >= t0, s0, s1, s0, s1, mapper);
}
public void checkValue(double value, double reference) {
Assert.assertTrue(FastMath.abs(value - reference) < 1.0e-10);
}
}

71
src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractRungeKuttaFieldStepInterpolatorTest.java
View File

@ -18,8 +18,6 @@
package org.apache.commons.math4.ode.nonstiff;
import java.lang.reflect.InvocationTargetException;
import org.apache.commons.math4.Field;
import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.AbstractIntegrator;
@ -38,7 +36,15 @@ import org.junit.Test;
public abstract class AbstractRungeKuttaFieldStepInterpolatorTest {
protected abstract <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper);
protected abstract <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field);
@Test
public abstract void interpolationAtBounds();
@ -140,10 +146,8 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest {
final double t0, final double[] y0,
final double t1) {
RungeKuttaFieldStepInterpolator<T> interpolator = createInterpolator(field, t1 > t0,
new FieldExpandableODE<T>(eqn).getMapper());
// get the Butcher arrays from the field integrator
FieldButcherArrayProvider<T> provider = createButcherArrayProvider(field, interpolator);
FieldButcherArrayProvider<T> provider = createButcherArrayProvider(field);
T[][] a = provider.getA();
T[] b = provider.getB();
T[] c = provider.getC();
@ -156,8 +160,7 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest {
fieldY[i] = field.getZero().add(y0[i]);
}
fieldYDotK[0] = eqn.computeDerivatives(t, fieldY);
interpolator.storeState(new FieldODEStateAndDerivative<T>(t, fieldY, fieldYDotK[0]));
interpolator.shift();
FieldODEStateAndDerivative<T> s0 = new FieldODEStateAndDerivative<T>(t, fieldY, fieldYDotK[0]);
// perform one integration step, in order to get consistent derivatives
T h = field.getZero().add(t1 - t0);
@ -170,20 +173,20 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest {
}
fieldYDotK[k + 1] = eqn.computeDerivatives(h.multiply(c[k]).add(t0), fieldY);
}
interpolator.setSlopes(fieldYDotK);
// store state at step end
t = field.getZero().add(t1);
for (int i = 0; i < y0.length; ++i) {
fieldY[i] = field.getZero().add(y0[i]);
for (int s = 0; s < b.length; ++s) {
fieldY[i] = fieldY[i].add(h.multiply(b[s].multiply(fieldYDotK[s][i])));
}
}
interpolator.storeState(new FieldODEStateAndDerivative<T>(field.getZero().add(t1),
fieldY,
eqn.computeDerivatives(field.getZero().add(t1), fieldY)));
FieldODEStateAndDerivative<T> s1 = new FieldODEStateAndDerivative<T>(t, fieldY,
eqn.computeDerivatives(t, fieldY));
return interpolator;
return createInterpolator(field, t1 > t0, fieldYDotK, s0, s1, s0, s1,
new FieldExpandableODE<T>(eqn).getMapper());
}
@ -236,10 +239,10 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest {
}
@Override
public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
T fieldT = fieldInterpolator.getField().getZero().add(t);
T[] fieldY = MathArrays.buildArray(fieldInterpolator.getField(), y.length);
T fieldT = fieldInterpolator.getCurrentState().getTime().getField().getZero().add(t);
T[] fieldY = MathArrays.buildArray(fieldInterpolator.getCurrentState().getTime().getField(), y.length);
for (int i = 0; i < y.length; ++i) {
fieldY[i] = fieldInterpolator.getField().getZero().add(y[i]);
fieldY[i] = fieldInterpolator.getCurrentState().getTime().getField().getZero().add(y[i]);
}
T[] fieldYDot = eqn.computeDerivatives(fieldT, fieldY);
for (int i = 0; i < yDot.length; ++i) {
@ -281,42 +284,6 @@ public abstract class AbstractRungeKuttaFieldStepInterpolatorTest {
}
private <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field, final RungeKuttaFieldStepInterpolator<T> provider) {
FieldButcherArrayProvider<T> integrator = null;
try {
String interpolatorName = provider.getClass().getName();
String integratorName = interpolatorName.replaceAll("StepInterpolator", "Integrator");
@SuppressWarnings("unchecked")
Class<FieldButcherArrayProvider<T>> clz = (Class<FieldButcherArrayProvider<T>>) Class.forName(integratorName);
try {
integrator = clz.getConstructor(Field.class, RealFieldElement.class).
newInstance(field, field.getOne());
} catch (NoSuchMethodException nsme) {
try {
integrator = clz.getConstructor(Field.class,
Double.TYPE, Double.TYPE,
Double.TYPE, Double.TYPE).
newInstance(field, 0.001, 1.0, 1.0, 1.0);
} catch (NoSuchMethodException e) {
Assert.fail(e.getLocalizedMessage());
}
}
} catch (InvocationTargetException ite) {
Assert.fail(ite.getLocalizedMessage());
} catch (IllegalAccessException iae) {
Assert.fail(iae.getLocalizedMessage());
} catch (InstantiationException ie) {
Assert.fail(ie.getLocalizedMessage());
} catch (ClassNotFoundException cnfe) {
Assert.fail(cnfe.getLocalizedMessage());
}
return integrator;
}
private static class SinCos<T extends RealFieldElement<T>> implements FieldFirstOrderDifferentialEquations<T> {
private final Field<T> field;
protected SinCos(final Field<T> field) {

20
src/test/java/org/apache/commons/math4/ode/nonstiff/ClassicalRungKuttaFieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class ClassicalRungKuttaFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new ClassicalRungeKuttaFieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new ClassicalRungeKuttaFieldIntegrator<T>(field, field.getOne());
}
@Test
@ -43,7 +57,7 @@ public class ClassicalRungKuttaFieldStepInterpolatorTest extends AbstractRungeKu
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 1.2e-16, 2.3e-16, 1.0e-50);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 1.2e-16, 3.4e-16, 2.1e-17);
}
}

20
src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince54FieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class DormandPrince54FieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new DormandPrince54FieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new DormandPrince54FieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new DormandPrince54FieldIntegrator<T>(field, 0, 1, 1, 1);
}
@Test
@ -43,7 +57,7 @@ public class DormandPrince54FieldStepInterpolatorTest extends AbstractRungeKutta
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 2.3e-16, 4.5e-16, 5.6e-17);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 2.8e-17, 2.3e-16, 5.6e-16, 5.6e-17);
}
}

18
src/test/java/org/apache/commons/math4/ode/nonstiff/DormandPrince853FieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class DormandPrince853FieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new DormandPrince853FieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new DormandPrince853FieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new DormandPrince853FieldIntegrator<T>(field, 0, 1, 1, 1);
}
@Test

19
src/test/java/org/apache/commons/math4/ode/nonstiff/EulerFieldStepInterpolatorTest.java
View File

@ -21,14 +21,27 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class EulerFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new EulerFieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState, FieldEquationsMapper<T> mapper) {
return new EulerFieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new EulerFieldIntegrator<T>(field, field.getOne());
}
@Test
@ -43,7 +56,7 @@ public class EulerFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepI
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.0e-50, 1.0e-50, 1.0e-50, 1.0e-50);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 7.0e-18, 1.0e-50, 1.0e-50, 1.0e-50);
}
}

20
src/test/java/org/apache/commons/math4/ode/nonstiff/GillFieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class GillFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new GillFieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new GillFieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new GillFieldIntegrator<T>(field, field.getOne());
}
@Test
@ -43,7 +57,7 @@ public class GillFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepIn
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 1.0e-50, 1.0e-50);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 3.4e-16, 2.1e-17);
}
}

20
src/test/java/org/apache/commons/math4/ode/nonstiff/HighamHall54FieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class HighamHall54FieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new HighamHall54FieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new HighamHall54FieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new HighamHall54FieldIntegrator<T>(field, 0, 1, 1, 1);
}
@Test
@ -43,7 +57,7 @@ public class HighamHall54FieldStepInterpolatorTest extends AbstractRungeKuttaFie
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 1.0e-50, 1.0e-50);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.3e-16, 1.0e-50, 3.6e-15, 2.1e-16);
}
}

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

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class LutherFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new LutherFieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new LutherFieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new LutherFieldIntegrator<T>(field, field.getOne());
}
@Test

20
src/test/java/org/apache/commons/math4/ode/nonstiff/MidpointFieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class MidpointFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new MidpointFieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new MidpointFieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new MidpointFieldIntegrator<T>(field, field.getOne());
}
@Test
@ -43,7 +57,7 @@ public class MidpointFieldStepInterpolatorTest extends AbstractRungeKuttaFieldSt
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.0e-50, 1.0e-50, 1.0e-50, 1.0e-50);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.0e-50, 1.0e-50, 7.0e-18);
}
}

20
src/test/java/org/apache/commons/math4/ode/nonstiff/ThreeEighthesFieldStepInterpolatorTest.java
View File

@ -21,14 +21,28 @@ 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.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
import org.apache.commons.math4.util.Decimal64Field;
import org.junit.Test;
public class ThreeEighthesFieldStepInterpolatorTest extends AbstractRungeKuttaFieldStepInterpolatorTest {
protected <T extends RealFieldElement<T>> RungeKuttaFieldStepInterpolator<T>
createInterpolator(Field<T> field, boolean forward, FieldEquationsMapper<T> mapper) {
return new ThreeEighthesFieldStepInterpolator<T>(field, forward, mapper);
createInterpolator(Field<T> field, boolean forward, T[][] yDotK,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldODEStateAndDerivative<T> softPreviousState,
FieldODEStateAndDerivative<T> softCurrentState,
FieldEquationsMapper<T> mapper) {
return new ThreeEighthesFieldStepInterpolator<T>(field, forward, yDotK,
globalPreviousState, globalCurrentState,
softPreviousState, softCurrentState,
mapper);
}
protected <T extends RealFieldElement<T>> FieldButcherArrayProvider<T>
createButcherArrayProvider(final Field<T> field) {
return new ThreeEighthesFieldIntegrator<T>(field, field.getOne());
}
@Test
@ -43,7 +57,7 @@ public class ThreeEighthesFieldStepInterpolatorTest extends AbstractRungeKuttaFi
@Test
public void nonFieldInterpolatorConsistency() {
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.2e-16, 1.0e-50, 1.0e-50);
doNonFieldInterpolatorConsistency(Decimal64Field.getInstance(), 1.4e-17, 1.2e-16, 3.4e-16, 1.4e-17);
}
}

55
src/test/java/org/apache/commons/math4/ode/sampling/DummyFieldStepInterpolator.java Normal file
View File

@ -0,0 +1,55 @@
/*
* 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.sampling;
import org.apache.commons.math4.RealFieldElement;
import org.apache.commons.math4.ode.FieldEquationsMapper;
import org.apache.commons.math4.ode.FieldODEStateAndDerivative;
public class DummyFieldStepInterpolator<T extends RealFieldElement<T>>
extends AbstractFieldStepInterpolator<T> {
public DummyFieldStepInterpolator(final boolean forward,
final FieldODEStateAndDerivative<T> globalPreviousState,
final FieldODEStateAndDerivative<T> globalCurrentState,
final FieldODEStateAndDerivative<T> softPreviousState,
final FieldODEStateAndDerivative<T> softCurrentState,
final FieldEquationsMapper<T> mapper) {
super(forward, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
}
@Override
protected AbstractFieldStepInterpolator<T> create(final boolean newForward,
final FieldODEStateAndDerivative<T> newGlobalPreviousState,
final FieldODEStateAndDerivative<T> newGlobalCurrentState,
final FieldODEStateAndDerivative<T> newSoftPreviousState,
final FieldODEStateAndDerivative<T> newSoftCurrentState,
final FieldEquationsMapper<T> newMapper) {
return new DummyFieldStepInterpolator<T>(newForward,
newGlobalPreviousState, newGlobalCurrentState,
newSoftPreviousState, newSoftCurrentState,
newMapper);
}
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(FieldEquationsMapper<T> equationsMapper,
T time, T theta, T thetaH, T oneMinusThetaH) {
return getGlobalCurrentState();
}
}