From 5b5f84db290450f81cc6b93cad7e7900ef9a6cf2 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Tue, 2 Apr 2013 19:04:13 +0000 Subject: [PATCH] Fixed wrong array dimensions in secondary equations handling. JIRA: MATH-961 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1463684 13f79535-47bb-0310-9956-ffa450edef68 --- src/changes/changes.xml | 3 + .../commons/math3/ode/AbstractIntegrator.java | 35 ++++++-- .../math3/ode/ContinuousOutputModel.java | 16 ++++ .../math3/ode/MultistepIntegrator.java | 80 ++++++++++--------- .../nonstiff/AdamsBashforthIntegrator.java | 10 ++- .../ode/nonstiff/AdamsMoultonIntegrator.java | 10 ++- 6 files changed, 111 insertions(+), 43 deletions(-) diff --git a/src/changes/changes.xml b/src/changes/changes.xml index c12a2a6cc..68eee9520 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -78,6 +78,9 @@ Users are encouraged to upgrade to this version as this release not 2. A few methods in the FastMath class are in fact slower that their counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901). "> + + Fixed wrong array dimensions in secondary equations handling in some cases. + Fixed missing side effects of secondary equations on main state in Ordinary Differential Equations integration. diff --git a/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java b/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java index 6e0237f9f..809230a1d 100644 --- a/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java @@ -206,6 +206,22 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { this.expandable = equations; } + /** Get the differential equations to integrate. + * @return differential equations to integrate + * @since 3.2 + */ + protected ExpandableStatefulODE getExpandable() { + return expandable; + } + + /** Get the evaluations counter. + * @return evaluations counter + * @since 3.2 + */ + protected Incrementor getEvaluationsCounter() { + return evaluations; + } + /** {@inheritDoc} */ public double integrate(final FirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) @@ -340,11 +356,19 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { // get state at event time interpolator.setInterpolatedTime(eventT); - final double[] eventY = interpolator.getInterpolatedState().clone(); + final double[] eventYPrimary = interpolator.getInterpolatedState().clone(); + final double[] eventYComplete = new double[y.length]; + expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(), + eventYComplete); + int index = 0; + for (EquationsMapper secondary : expandable.getSecondaryMappers()) { + secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++), + eventYComplete); + } // advance all event states to current time for (final EventState state : eventsStates) { - state.stepAccepted(eventT, eventY); + state.stepAccepted(eventT, eventYPrimary); isLastStep = isLastStep || state.stop(); } @@ -355,18 +379,19 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { if (isLastStep) { // the event asked to stop integration - System.arraycopy(eventY, 0, y, 0, y.length); + System.arraycopy(eventYComplete, 0, y, 0, y.length); return eventT; } boolean needReset = false; for (final EventState state : eventsStates) { - needReset = needReset || state.reset(eventT, eventY); + needReset = needReset || state.reset(eventT, eventYComplete); } if (needReset) { // some event handler has triggered changes that // invalidate the derivatives, we need to recompute them - System.arraycopy(eventY, 0, y, 0, y.length); + interpolator.setInterpolatedTime(eventT); + System.arraycopy(eventYComplete, 0, y, 0, y.length); computeDerivatives(eventT, y, yDot); resetOccurred = true; return eventT; diff --git a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java index 670b1830d..a5f5faf31 100644 --- a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java +++ b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java @@ -333,11 +333,27 @@ public class ContinuousOutputModel * Get the state vector of the interpolated point. * @return state vector at time {@link #getInterpolatedTime} * @exception MaxCountExceededException if the number of functions evaluations is exceeded + * @see #getInterpolatedSecondaryState(int) */ public double[] getInterpolatedState() throws MaxCountExceededException { return steps.get(index).getInterpolatedState(); } + /** Get the interpolated secondary state corresponding to the secondary equations. + * @param secondaryStateIndex index of the secondary set, as returned by {@link + * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( + * org.apache.commons.math3.ode.SecondaryEquations) + * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)} + * @return interpolated secondary state at the current interpolation date + * @see #getInterpolatedState() + * @since 3.2 + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + */ + public double[] getInterpolatedSecondaryState(final int secondaryStateIndex) + throws MaxCountExceededException { + return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex); + } + /** Compare a step interval and a double. * @param time point to locate * @param interval step interval diff --git a/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java b/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java index ff8d69120..a729441ef 100644 --- a/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java @@ -228,10 +228,31 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { // start integration, expecting a InitializationCompletedMarkerException try { - starter.integrate(new CountingDifferentialEquations(y0.length), - t0, y0, t, new double[y0.length]); + + if (starter instanceof AbstractIntegrator) { + ((AbstractIntegrator) starter).integrate(getExpandable(), t); + } else { + starter.integrate(new FirstOrderDifferentialEquations() { + + /** {@inheritDoc} */ + public int getDimension() { + return getExpandable().getTotalDimension(); + } + + /** {@inheritDoc} */ + public void computeDerivatives(double t, double[] y, double[] yDot) { + getExpandable().computeDerivatives(t, y, yDot); + } + + }, t0, y0, t, new double[y0.length]); + } + } catch (InitializationCompletedMarkerException icme) { // NOPMD // this is the expected nominal interruption of the start integrator + + // count the evaluations used by the starter + getEvaluationsCounter().incrementCount(starter.getEvaluations()); + } // remove the specific step handler @@ -353,20 +374,33 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { // first step, we need to store also the beginning of the step interpolator.setInterpolatedTime(prev); t[0] = prev; - System.arraycopy(interpolator.getInterpolatedState(), 0, - y[0], 0, y[0].length); - System.arraycopy(interpolator.getInterpolatedDerivatives(), 0, - yDot[0], 0, yDot[0].length); + final ExpandableStatefulODE expandable = getExpandable(); + final EquationsMapper primary = expandable.getPrimaryMapper(); + primary.insertEquationData(interpolator.getInterpolatedState(), y[count]); + primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]); + int index = 0; + for (final EquationsMapper secondary : expandable.getSecondaryMappers()) { + secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]); + secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]); + ++index; + } } // store the end of the step ++count; interpolator.setInterpolatedTime(curr); t[count] = curr; - System.arraycopy(interpolator.getInterpolatedState(), 0, - y[count], 0, y[count].length); - System.arraycopy(interpolator.getInterpolatedDerivatives(), 0, - yDot[count], 0, yDot[count].length); + + final ExpandableStatefulODE expandable = getExpandable(); + final EquationsMapper primary = expandable.getPrimaryMapper(); + primary.insertEquationData(interpolator.getInterpolatedState(), y[count]); + primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]); + int index = 0; + for (final EquationsMapper secondary : expandable.getSecondaryMappers()) { + secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]); + secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]); + ++index; + } if (count == t.length - 1) { @@ -411,30 +445,4 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { } - /** Wrapper for differential equations, ensuring start evaluations are counted. */ - private class CountingDifferentialEquations implements FirstOrderDifferentialEquations { - - /** Dimension of the problem. */ - private final int dimension; - - /** Simple constructor. - * @param dimension dimension of the problem - */ - public CountingDifferentialEquations(final int dimension) { - this.dimension = dimension; - } - - /** {@inheritDoc} */ - public void computeDerivatives(double t, double[] y, double[] dot) - throws MaxCountExceededException, DimensionMismatchException { - MultistepIntegrator.this.computeDerivatives(t, y, dot); - } - - /** {@inheritDoc} */ - public int getDimension() { - return dimension; - } - - } - } diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java index 4b6710be5..d15dc7391 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java @@ -22,6 +22,7 @@ import org.apache.commons.math3.exception.MaxCountExceededException; import org.apache.commons.math3.exception.NoBracketingException; import org.apache.commons.math3.exception.NumberIsTooSmallException; import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.ode.EquationsMapper; import org.apache.commons.math3.ode.ExpandableStatefulODE; import org.apache.commons.math3.ode.sampling.NordsieckStepInterpolator; import org.apache.commons.math3.util.FastMath; @@ -255,7 +256,14 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator { final double stepEnd = stepStart + stepSize; interpolator.shift(); interpolator.setInterpolatedTime(stepEnd); - System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length); + final ExpandableStatefulODE expandable = getExpandable(); + final EquationsMapper primary = expandable.getPrimaryMapper(); + primary.insertEquationData(interpolator.getInterpolatedState(), y); + int index = 0; + for (final EquationsMapper secondary : expandable.getSecondaryMappers()) { + secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y); + ++index; + } // evaluate the derivative computeDerivatives(stepEnd, y, yDot); diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java index a58716528..5d5a7a7d4 100644 --- a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java @@ -25,6 +25,7 @@ import org.apache.commons.math3.exception.NoBracketingException; import org.apache.commons.math3.exception.NumberIsTooSmallException; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.RealMatrixPreservingVisitor; +import org.apache.commons.math3.ode.EquationsMapper; import org.apache.commons.math3.ode.ExpandableStatefulODE; import org.apache.commons.math3.ode.sampling.NordsieckStepInterpolator; import org.apache.commons.math3.util.FastMath; @@ -250,7 +251,14 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator { // predict a first estimate of the state at step end (P in the PECE sequence) final double stepEnd = stepStart + stepSize; interpolator.setInterpolatedTime(stepEnd); - System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length); + final ExpandableStatefulODE expandable = getExpandable(); + final EquationsMapper primary = expandable.getPrimaryMapper(); + primary.insertEquationData(interpolator.getInterpolatedState(), yTmp); + int index = 0; + for (final EquationsMapper secondary : expandable.getSecondaryMappers()) { + secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), yTmp); + ++index; + } // evaluate a first estimate of the derivative (first E in the PECE sequence) computeDerivatives(stepEnd, yTmp, yDot);