added the two methods getCurrentStepStart and getCurrentStepsize to interface FirstOrderIntegrator allowing ODE problems to retrieve the current step start and size during integration for each step trial (i.e. even before the step is accepted)

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@580753 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2007-09-30 16:58:50 +00:00
parent e92a93d217
commit 67b86d4ccf
5 changed files with 138 additions and 65 deletions

View File

@ -75,6 +75,8 @@ public abstract class AdaptiveStepsizeIntegrator
switchesHandler = new SwitchingFunctionsHandler();
resetInternalState();
}
/** Build an integrator with the given stepsize bounds.
@ -104,6 +106,8 @@ public abstract class AdaptiveStepsizeIntegrator
switchesHandler = new SwitchingFunctionsHandler();
resetInternalState();
}
/** Set the initial step size.
@ -324,6 +328,20 @@ public abstract class AdaptiveStepsizeIntegrator
double t, double[] y)
throws DerivativeException, IntegratorException;
public double getCurrentStepStart() {
return stepStart;
}
public double getCurrentStepsize() {
return stepSize;
}
/** Reset internal state to dummy values. */
protected void resetInternalState() {
stepStart = Double.NaN;
stepSize = Math.sqrt(minStep * maxStep);
}
/** Get the minimal step.
* @return minimal step
*/
@ -365,4 +383,10 @@ public abstract class AdaptiveStepsizeIntegrator
/** Switching functions handler. */
protected SwitchingFunctionsHandler switchesHandler;
/** Current step start time. */
protected double stepStart;
/** Current stepsize. */
protected double stepSize;
}

View File

@ -63,12 +63,16 @@ public interface FirstOrderIntegrator {
double maxCheckInterval,
double convergence);
/** Integrate the differential equations up to the given time
/** Integrate the differential equations up to the given time.
* <p>This method solves an Initial Value Problem (IVP).</p>
* <p>Since this method stores some internal state variables made
* available in its public interface during integration ({@link
* #getCurrentStepsize()}), it is <em>not</em> thread-safe.</p>
* @param equations differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param t target time for the integration
* (can be set to a value smaller thant <code>t0</code> for backward integration)
* (can be set to a value smaller than <code>t0</code> for backward integration)
* @param y placeholder where to put the state vector at each successful
* step (and hence at the end of integration), can be the same object as y0
* @throws IntegratorException if the integrator cannot perform integration
@ -80,4 +84,26 @@ public interface FirstOrderIntegrator {
double t, double[] y)
throws DerivativeException, IntegratorException;
/** Get the current value of the step start time t<sub>i</sub>.
* <p>This method can be called during integration (typically by
* the object implementing the {@link FirstOrderDifferentialEquations
* differential equations} problem) if the value of the current step that
* is attempted is needed.</p>
* <p>The result is undefined if the method is called outside of
* calls to {@link #integrate}</p>
* @return current value of the step start time t<sub>i</sub>
*/
public double getCurrentStepStart();
/** Get the current value of the integration stepsize.
* <p>This method can be called during integration (typically by
* the object implementing the {@link FirstOrderDifferentialEquations
* differential equations} problem) if the value of the current stepsize
* that is tried is needed.</p>
* <p>The result is undefined if the method is called outside of
* calls to {@link #integrate}</p>
* @return current value of the stepsize
*/
public double getCurrentStepsize();
}

View File

@ -579,7 +579,7 @@ public class GraggBulirschStoerIntegrator
}
interpolator.storeTime(t0);
double currentT = t0;
stepStart = t0;
double hNew = 0;
double maxError = Double.MAX_VALUE;
boolean previousRejected = false;
@ -591,7 +591,6 @@ public class GraggBulirschStoerIntegrator
costPerTimeUnit[0] = 0;
while (! lastStep) {
double h;
double error;
boolean reject = false;
@ -601,14 +600,14 @@ public class GraggBulirschStoerIntegrator
// first evaluation, at the beginning of the step
if (! firstStepAlreadyComputed) {
equations.computeDerivatives(currentT, y, yDot0);
equations.computeDerivatives(stepStart, y, yDot0);
}
if (firstTime) {
hNew = initializeStep(equations, forward,
2 * targetIter + 1, scale,
currentT, y, yDot0, yTmp, yTmpDot);
stepStart, y, yDot0, yTmp, yTmpDot);
if (! forward) {
hNew = -hNew;
@ -620,14 +619,14 @@ public class GraggBulirschStoerIntegrator
}
h = hNew;
stepSize = hNew;
// step adjustment near bounds
if ((forward && (currentT + h > t))
|| ((! forward) && (currentT + h < t))) {
h = t - currentT;
if ((forward && (stepStart + stepSize > t))
|| ((! forward) && (stepStart + stepSize < t))) {
stepSize = t - stepStart;
}
double nextT = currentT + h;
double nextT = stepStart + stepSize;
lastStep = forward ? (nextT >= t) : (nextT <= t);
// iterate over several substep sizes
@ -637,13 +636,13 @@ public class GraggBulirschStoerIntegrator
++k;
// modified midpoint integration with the current substep
if ( ! tryStep(equations, currentT, y, h, k, scale, fk[k],
if ( ! tryStep(equations, stepStart, y, stepSize, k, scale, fk[k],
(k == 0) ? yMidDots[0] : diagonal[k-1],
(k == 0) ? y1 : y1Diag[k-1],
yTmp)) {
// the stability check failed, we reduce the global step
hNew = Math.abs(filterStep(h * stabilityReduction, false));
hNew = Math.abs(filterStep(stepSize * stabilityReduction, false));
reject = true;
loop = false;
@ -667,7 +666,7 @@ public class GraggBulirschStoerIntegrator
if ((error > 1.0e15) || ((k > 1) && (error > maxError))) {
// error is too big, we reduce the global step
hNew = Math.abs(filterStep(h * stabilityReduction, false));
hNew = Math.abs(filterStep(stepSize * stabilityReduction, false));
reject = true;
loop = false;
} else {
@ -679,7 +678,7 @@ public class GraggBulirschStoerIntegrator
double fac = stepControl2 / Math.pow(error / stepControl1, exp);
double pow = Math.pow(stepControl3, exp);
fac = Math.max(pow / stepControl4, Math.min(1 / pow, fac));
optimalStep[k] = Math.abs(filterStep(h * fac, true));
optimalStep[k] = Math.abs(filterStep(stepSize * fac, true));
costPerTimeUnit[k] = costPerStep[k] / optimalStep[k];
// check convergence
@ -775,7 +774,7 @@ public class GraggBulirschStoerIntegrator
}
// derivative at end of step
equations.computeDerivatives(currentT + h, y1, yDot1);
equations.computeDerivatives(stepStart + stepSize, y1, yDot1);
int mu = 2 * k - mudif + 3;
@ -797,7 +796,7 @@ public class GraggBulirschStoerIntegrator
extrapolate(l2, j, diagonal, yMidDots[l+1]);
}
for (int i = 0; i < y0.length; ++i) {
yMidDots[l+1][i] *= h;
yMidDots[l+1][i] *= stepSize;
}
// compute centered differences to evaluate next derivatives
@ -816,13 +815,13 @@ public class GraggBulirschStoerIntegrator
// estimate the dense output coefficients
GraggBulirschStoerStepInterpolator gbsInterpolator
= (GraggBulirschStoerStepInterpolator) interpolator;
gbsInterpolator.computeCoefficients(mu, h);
gbsInterpolator.computeCoefficients(mu, stepSize);
if (useInterpolationError) {
// use the interpolation error to limit stepsize
double interpError = gbsInterpolator.estimateError(scale);
hInt = Math.abs(h / Math.max(Math.pow(interpError, 1.0 / (mu+4)),
0.01));
hInt = Math.abs(stepSize / Math.max(Math.pow(interpError, 1.0 / (mu+4)),
0.01));
if (interpError > 10.0) {
hNew = hInt;
reject = true;
@ -831,10 +830,10 @@ public class GraggBulirschStoerIntegrator
// Switching functions handling
if (!reject) {
interpolator.storeTime(currentT + h);
interpolator.storeTime(stepStart + stepSize);
if (switchesHandler.evaluateStep(interpolator)) {
reject = true;
hNew = Math.abs(switchesHandler.getEventTime() - currentT);
hNew = Math.abs(switchesHandler.getEventTime() - stepStart);
}
}
@ -851,19 +850,19 @@ public class GraggBulirschStoerIntegrator
if (! reject) {
// store end of step state
currentT += h;
stepStart += stepSize;
System.arraycopy(y1, 0, y, 0, y0.length);
switchesHandler.stepAccepted(currentT, y);
switchesHandler.stepAccepted(stepStart, y);
if (switchesHandler.stop()) {
lastStep = true;
}
// provide the step data to the step handler
interpolator.storeTime(currentT);
interpolator.storeTime(stepStart);
handler.handleStep(interpolator, lastStep);
if (switchesHandler.reset(currentT, y) && ! lastStep) {
if (switchesHandler.reset(stepStart, y) && ! lastStep) {
// some switching function has triggered changes that
// invalidate the derivatives, we need to recompute them
firstStepAlreadyComputed = false;
@ -897,7 +896,7 @@ public class GraggBulirschStoerIntegrator
// after a rejected step neither order nor stepsize
// should increase
targetIter = Math.min(optimalIter, k);
hNew = Math.min(Math.abs(h), optimalStep[targetIter]);
hNew = Math.min(Math.abs(stepSize), optimalStep[targetIter]);
} else {
// stepsize control
if (optimalIter <= k) {

View File

@ -187,7 +187,7 @@ public abstract class RungeKuttaFehlbergIntegrator
}
interpolator.storeTime(t0);
double currentT = t0;
stepStart = t0;
double hNew = 0;
boolean firstTime = true;
boolean lastStep;
@ -196,13 +196,12 @@ public abstract class RungeKuttaFehlbergIntegrator
interpolator.shift();
double h = 0;
double error = 0;
for (boolean loop = true; loop;) {
if (firstTime || !fsal) {
// first stage
equations.computeDerivatives(currentT, y, yDotK[0]);
equations.computeDerivatives(stepStart, y, yDotK[0]);
}
if (firstTime) {
@ -216,16 +215,16 @@ public abstract class RungeKuttaFehlbergIntegrator
}
}
hNew = initializeStep(equations, forward, getOrder(), scale,
currentT, y, yDotK[0], yTmp, yDotK[1]);
stepStart, y, yDotK[0], yTmp, yDotK[1]);
firstTime = false;
}
h = hNew;
stepSize = hNew;
// step adjustment near bounds
if ((forward && (currentT + h > t))
|| ((! forward) && (currentT + h < t))) {
h = t - currentT;
if ((forward && (stepStart + stepSize > t))
|| ((! forward) && (stepStart + stepSize < t))) {
stepSize = t - stepStart;
}
// next stages
@ -236,10 +235,10 @@ public abstract class RungeKuttaFehlbergIntegrator
for (int l = 1; l < k; ++l) {
sum += a[k-1][l] * yDotK[l][j];
}
yTmp[j] = y[j] + h * sum;
yTmp[j] = y[j] + stepSize * sum;
}
equations.computeDerivatives(currentT + c[k-1] * h, yTmp, yDotK[k]);
equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
}
@ -249,18 +248,18 @@ public abstract class RungeKuttaFehlbergIntegrator
for (int l = 1; l < stages; ++l) {
sum += b[l] * yDotK[l][j];
}
yTmp[j] = y[j] + h * sum;
yTmp[j] = y[j] + stepSize * sum;
}
// estimate the error at the end of the step
error = estimateError(yDotK, y, yTmp, h);
error = estimateError(yDotK, y, yTmp, stepSize);
if (error <= 1.0) {
// Switching functions handling
interpolator.storeTime(currentT + h);
interpolator.storeTime(stepStart + stepSize);
if (switchesHandler.evaluateStep(interpolator)) {
// reject the step to match exactly the next switch time
hNew = switchesHandler.getEventTime() - currentT;
hNew = switchesHandler.getEventTime() - stepStart;
} else {
// accept the step
loop = false;
@ -271,23 +270,23 @@ public abstract class RungeKuttaFehlbergIntegrator
double factor = Math.min(maxGrowth,
Math.max(minReduction,
safety * Math.pow(error, exp)));
hNew = filterStep(h * factor, false);
hNew = filterStep(stepSize * factor, false);
}
}
// the step has been accepted
currentT += h;
stepStart += stepSize;
System.arraycopy(yTmp, 0, y, 0, y0.length);
switchesHandler.stepAccepted(currentT, y);
switchesHandler.stepAccepted(stepStart, y);
if (switchesHandler.stop()) {
lastStep = true;
} else {
lastStep = forward ? (currentT >= t) : (currentT <= t);
lastStep = forward ? (stepStart >= t) : (stepStart <= t);
}
// provide the step data to the step handler
interpolator.storeTime(currentT);
interpolator.storeTime(stepStart);
handler.handleStep(interpolator, lastStep);
if (fsal) {
@ -295,10 +294,10 @@ public abstract class RungeKuttaFehlbergIntegrator
System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
}
if (switchesHandler.reset(currentT, y) && ! lastStep) {
if (switchesHandler.reset(stepStart, y) && ! lastStep) {
// some switching function has triggered changes that
// invalidate the derivatives, we need to recompute them
equations.computeDerivatives(currentT, y, yDotK[0]);
equations.computeDerivatives(stepStart, y, yDotK[0]);
}
if (! lastStep) {
@ -306,14 +305,16 @@ public abstract class RungeKuttaFehlbergIntegrator
double factor = Math.min(maxGrowth,
Math.max(minReduction,
safety * Math.pow(error, exp)));
double scaledH = h * factor;
double nextT = currentT + scaledH;
double scaledH = stepSize * factor;
double nextT = stepStart + scaledH;
boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
hNew = filterStep(scaledH, nextIsLast);
}
} while (! lastStep);
resetInternalState();
}
/** Get the minimal reduction factor for stepsize control.

View File

@ -78,6 +78,7 @@ public abstract class RungeKuttaIntegrator
this.step = step;
handler = DummyStepHandler.getInstance();
switchesHandler = new SwitchingFunctionsHandler();
resetInternalState();
}
/** Get the name of the method.
@ -180,11 +181,11 @@ public abstract class RungeKuttaIntegrator
interpolator.storeTime(t0);
// recompute the step
double currentT = t0;
long nbStep = Math.max(1l, Math.abs(Math.round((t - t0) / step)));
double h = (t - t0) / nbStep;
boolean firstTime = true;
boolean lastStep = false;
stepStart = t0;
stepSize = (t - t0) / nbStep;
handler.reset();
for (long i = 0; ! lastStep; ++i) {
@ -195,7 +196,7 @@ public abstract class RungeKuttaIntegrator
if (firstTime || !fsal) {
// first stage
equations.computeDerivatives(currentT, y, yDotK[0]);
equations.computeDerivatives(stepStart, y, yDotK[0]);
firstTime = false;
}
@ -207,10 +208,10 @@ public abstract class RungeKuttaIntegrator
for (int l = 1; l < k; ++l) {
sum += a[k-1][l] * yDotK[l][j];
}
yTmp[j] = y[j] + h * sum;
yTmp[j] = y[j] + stepSize * sum;
}
equations.computeDerivatives(currentT + c[k-1] * h, yTmp, yDotK[k]);
equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
}
@ -220,14 +221,14 @@ public abstract class RungeKuttaIntegrator
for (int l = 1; l < stages; ++l) {
sum += b[l] * yDotK[l][j];
}
yTmp[j] = y[j] + h * sum;
yTmp[j] = y[j] + stepSize * sum;
}
// Switching functions handling
interpolator.storeTime(currentT + h);
interpolator.storeTime(stepStart + stepSize);
if (switchesHandler.evaluateStep(interpolator)) {
needUpdate = true;
h = switchesHandler.getEventTime() - currentT;
stepSize = switchesHandler.getEventTime() - stepStart;
} else {
loop = false;
}
@ -235,9 +236,9 @@ public abstract class RungeKuttaIntegrator
}
// the step has been accepted
currentT += h;
stepStart += stepSize;
System.arraycopy(yTmp, 0, y, 0, y0.length);
switchesHandler.stepAccepted(currentT, y);
switchesHandler.stepAccepted(stepStart, y);
if (switchesHandler.stop()) {
lastStep = true;
} else {
@ -245,7 +246,7 @@ public abstract class RungeKuttaIntegrator
}
// provide the step data to the step handler
interpolator.storeTime(currentT);
interpolator.storeTime(stepStart);
handler.handleStep(interpolator, lastStep);
if (fsal) {
@ -253,22 +254,38 @@ public abstract class RungeKuttaIntegrator
System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
}
if (switchesHandler.reset(currentT, y) && ! lastStep) {
if (switchesHandler.reset(stepStart, y) && ! lastStep) {
// some switching function has triggered changes that
// invalidate the derivatives, we need to recompute them
equations.computeDerivatives(currentT, y, yDotK[0]);
equations.computeDerivatives(stepStart, y, yDotK[0]);
}
if (needUpdate) {
// a switching function has changed the step
// we need to recompute stepsize
nbStep = Math.max(1l, Math.abs(Math.round((t - currentT) / step)));
h = (t - currentT) / nbStep;
nbStep = Math.max(1l, Math.abs(Math.round((t - stepStart) / step)));
stepSize = (t - stepStart) / nbStep;
i = -1;
}
}
resetInternalState();
}
public double getCurrentStepStart() {
return stepStart;
}
public double getCurrentStepsize() {
return stepSize;
}
/** Reset internal state to dummy values. */
private void resetInternalState() {
stepStart = Double.NaN;
stepSize = Double.NaN;
}
/** Indicator for <i>fsal</i> methods. */
@ -295,4 +312,10 @@ public abstract class RungeKuttaIntegrator
/** Switching functions handler. */
protected SwitchingFunctionsHandler switchesHandler;
/** Current step start time. */
private double stepStart;
/** Current stepsize. */
private double stepSize;
}