Apache/commons-math
1
0
Fork 0
mirror of https://github.com/apache/commons-math.git synced 2025-02-06 10:09:26 +00:00
Code Issues Packages Projects Releases Wiki Activity

Replaced size adjustment of all steps of fixed steps Runge-Kutta integrators by a truncation of the last step only.

Browse Source 
JIRA: MATH-214

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@675551 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-07-10 12:40:44 +00:00
parent f6c3d3f104
commit b0b7c6ef38
8 changed files with 187 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
View File

@ -22,6 +22,7 @@ import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
@ -106,19 +107,20 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
}
interpolator.storeTime(t0);
// recompute the step
long nbStep = Math.max(1l, Math.abs(Math.round((t - t0) / step)));
boolean lastStep = false;
// set up integration control objects
stepStart = t0;
stepSize = (t - t0) / nbStep;
stepSize = step;
for (StepHandler handler : stepHandlers) {
handler.reset();
}
for (long i = 0; ! lastStep; ++i) {
CombinedEventsManager manager = addEndTimeChecker(t, eventsHandlersManager);
boolean lastStep = false;
// main integration loop
while (!lastStep) {
interpolator.shift();
boolean needUpdate = false;
for (boolean loop = true; loop;) {
// first stage
@ -148,11 +150,10 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
yTmp[j] = y[j] + stepSize * sum;
}
// Discrete events handling
// discrete events handling
interpolator.storeTime(stepStart + stepSize);
if (eventsHandlersManager.evaluateStep(interpolator)) {
needUpdate = true;
stepSize = eventsHandlersManager.getEventTime() - stepStart;
if (manager.evaluateStep(interpolator)) {
stepSize = manager.getEventTime() - stepStart;
} else {
loop = false;
}
@ -162,12 +163,8 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
// the step has been accepted
final double nextStep = stepStart + stepSize;
System.arraycopy(yTmp, 0, y, 0, y0.length);
eventsHandlersManager.stepAccepted(nextStep, y);
if (eventsHandlersManager.stop()) {
lastStep = true;
} else {
lastStep = (i == (nbStep - 1));
}
manager.stepAccepted(nextStep, y);
lastStep = manager.stop();
// provide the step data to the step handler
interpolator.storeTime(nextStep);
@ -176,19 +173,14 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
}
stepStart = nextStep;
if (eventsHandlersManager.reset(stepStart, y) && ! lastStep) {
if (manager.reset(stepStart, y) && ! lastStep) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to recompute them
equations.computeDerivatives(stepStart, y, yDotK[0]);
}
if (needUpdate) {
// an event handler has changed the step
// we need to recompute stepsize
nbStep = Math.max(1l, Math.abs(Math.round((t - stepStart) / step)));
stepSize = (t - stepStart) / nbStep;
i = -1;
}
// make sure step size is set to default before next step
stepSize = step;
}

4
src/site/xdoc/changes.xml
View File

@ -39,6 +39,10 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="fix" issue="MATH-214" >
Replaced size adjustment of all steps of fixed steps Runge-Kutta integrators by
a truncation of the last step only.
</action>
<action dev="luc" type="update">
The ODE integrators now support several step handlers at once, instead of just one.
This is more consistent with event handlers management.

31
src/test/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
View File

@ -20,6 +20,7 @@ package org.apache.commons.math.ode.nonstiff;
import junit.framework.*;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.EventHandler;
@ -191,6 +192,36 @@ public class ClassicalRungeKuttaIntegratorTest
private TestProblem3 pb;
}
public void testStepSize()
throws DerivativeException, IntegratorException {
final double step = 1.23456;
FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
integ.addStepHandler(new StepHandler() {
private static final long serialVersionUID = 0L;
public void handleStep(StepInterpolator interpolator, boolean isLast) {
if (! isLast) {
assertEquals(step,
interpolator.getCurrentTime() - interpolator.getPreviousTime(),
1.0e-12);
}
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
}
});
integ.integrate(new FirstOrderDifferentialEquations() {
private static final long serialVersionUID = 0L;
public void computeDerivatives(double t, double[] y, double[] dot) {
dot[0] = 1.0;
}
public int getDimension() {
return 1;
}
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
public static Test suite() {
return new TestSuite(ClassicalRungeKuttaIntegratorTest.class);
}

33
src/test/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java
View File

@ -20,10 +20,13 @@ package org.apache.commons.math.ode.nonstiff;
import junit.framework.*;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.nonstiff.EulerIntegrator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;
public class EulerIntegratorTest
extends TestCase {
@ -124,6 +127,36 @@ public class EulerIntegratorTest
}
public void testStepSize()
throws DerivativeException, IntegratorException {
final double step = 1.23456;
FirstOrderIntegrator integ = new EulerIntegrator(step);
integ.addStepHandler(new StepHandler() {
private static final long serialVersionUID = 0L;
public void handleStep(StepInterpolator interpolator, boolean isLast) {
if (! isLast) {
assertEquals(step,
interpolator.getCurrentTime() - interpolator.getPreviousTime(),
1.0e-12);
}
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
}
});
integ.integrate(new FirstOrderDifferentialEquations() {
private static final long serialVersionUID = 0L;
public void computeDerivatives(double t, double[] y, double[] dot) {
dot[0] = 1.0;
}
public int getDimension() {
return 1;
}
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
public static Test suite() {
return new TestSuite(EulerIntegratorTest.class);
}

31
src/test/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
View File

@ -20,6 +20,7 @@ package org.apache.commons.math.ode.nonstiff;
import junit.framework.*;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.EventHandler;
@ -181,6 +182,36 @@ public class GillIntegratorTest
private TestProblem3 pb;
}
public void testStepSize()
throws DerivativeException, IntegratorException {
final double step = 1.23456;
FirstOrderIntegrator integ = new GillIntegrator(step);
integ.addStepHandler(new StepHandler() {
private static final long serialVersionUID = 0L;
public void handleStep(StepInterpolator interpolator, boolean isLast) {
if (! isLast) {
assertEquals(step,
interpolator.getCurrentTime() - interpolator.getPreviousTime(),
1.0e-12);
}
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
}
});
integ.integrate(new FirstOrderDifferentialEquations() {
private static final long serialVersionUID = 0L;
public void computeDerivatives(double t, double[] y, double[] dot) {
dot[0] = 1.0;
}
public int getDimension() {
return 1;
}
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
public static Test suite() {
return new TestSuite(GillIntegratorTest.class);
}

33
src/test/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java
View File

@ -20,10 +20,13 @@ package org.apache.commons.math.ode.nonstiff;
import junit.framework.*;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.nonstiff.MidpointIntegrator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;
public class MidpointIntegratorTest
extends TestCase {
@ -124,6 +127,36 @@ public class MidpointIntegratorTest
}
public void testStepSize()
throws DerivativeException, IntegratorException {
final double step = 1.23456;
FirstOrderIntegrator integ = new MidpointIntegrator(step);
integ.addStepHandler(new StepHandler() {
private static final long serialVersionUID = 0L;
public void handleStep(StepInterpolator interpolator, boolean isLast) {
if (! isLast) {
assertEquals(step,
interpolator.getCurrentTime() - interpolator.getPreviousTime(),
1.0e-12);
}
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
}
});
integ.integrate(new FirstOrderDifferentialEquations() {
private static final long serialVersionUID = 0L;
public void computeDerivatives(double t, double[] y, double[] dot) {
dot[0] = 1.0;
}
public int getDimension() {
return 1;
}
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
public static Test suite() {
return new TestSuite(MidpointIntegratorTest.class);
}

7
src/test/org/apache/commons/math/ode/nonstiff/StepInterpolatorAbstractTest.java
View File

@ -48,6 +48,10 @@ public class StepInterpolatorAbstractTest extends TestCase {
final double h = 0.001 * (interpolator.getCurrentTime() - interpolator.getPreviousTime());
final double t = interpolator.getCurrentTime() - 300 * h;
if (Math.abs(h) < 10 * Math.ulp(t)) {
return;
}
interpolator.setInterpolatedTime(t - 4 * h);
final double[] yM4h = interpolator.getInterpolatedState().clone();
interpolator.setInterpolatedTime(t - 3 * h);
@ -73,6 +77,9 @@ public class StepInterpolatorAbstractTest extends TestCase {
32 * (yP3h[i] - yM3h[i]) +
-168 * (yP2h[i] - yM2h[i]) +
672 * (yP1h[i] - yM1h[i])) / (840 * h);
if (Math.abs(approYDot - yDot[i]) >= threshold) {
System.out.println("gotcha!");
}
assertEquals(approYDot, yDot[i], threshold);
}

31
src/test/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
View File

@ -20,6 +20,7 @@ package org.apache.commons.math.ode.nonstiff;
import junit.framework.*;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.EventHandler;
@ -178,6 +179,36 @@ public class ThreeEighthesIntegratorTest
}
public void testStepSize()
throws DerivativeException, IntegratorException {
final double step = 1.23456;
FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step);
integ.addStepHandler(new StepHandler() {
private static final long serialVersionUID = 0L;
public void handleStep(StepInterpolator interpolator, boolean isLast) {
if (! isLast) {
assertEquals(step,
interpolator.getCurrentTime() - interpolator.getPreviousTime(),
1.0e-12);
}
}
public boolean requiresDenseOutput() {
return false;
}
public void reset() {
}
});
integ.integrate(new FirstOrderDifferentialEquations() {
private static final long serialVersionUID = 0L;
public void computeDerivatives(double t, double[] y, double[] dot) {
dot[0] = 1.0;
}
public int getDimension() {
return 1;
}
}, 0.0, new double[] { 0.0 }, 5.0, new double[1]);
}
public static Test suite() {
return new TestSuite(ThreeEighthesIntegratorTest.class);
}