From 92f005ad13c989d70da393072c4f3d66798e4a1b Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Fri, 23 Jul 2010 08:59:37 +0000 Subject: [PATCH] Added a feature allowing error estimation to be computed only on a subset of Ordinary Differential Equations, considered as the main set, the remaining equations being considered only as an extension set that should not influence the ODE integration algorithm JIRA: MATH-388 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@967007 13f79535-47bb-0310-9956-ffa450edef68 --- ...tendedFirstOrderDifferentialEquations.java | 66 +++++++++++++++++++ .../commons/math/ode/MultistepIntegrator.java | 7 +- .../FirstOrderIntegratorWithJacobians.java | 11 +++- .../math/ode/jacobians/ParameterizedODE.java | 3 +- .../nonstiff/AdamsBashforthIntegrator.java | 4 +- .../ode/nonstiff/AdamsMoultonIntegrator.java | 18 ++--- .../nonstiff/AdaptiveStepsizeIntegrator.java | 38 ++++++++--- .../nonstiff/DormandPrince54Integrator.java | 4 +- .../nonstiff/DormandPrince853Integrator.java | 4 +- .../EmbeddedRungeKuttaIntegrator.java | 2 +- .../GraggBulirschStoerIntegrator.java | 14 ++-- .../GraggBulirschStoerStepInterpolator.java | 4 +- .../ode/nonstiff/HighamHall54Integrator.java | 4 +- src/site/xdoc/changes.xml | 6 ++ ...FirstOrderIntegratorWithJacobiansTest.java | 56 ++++++++-------- 15 files changed, 171 insertions(+), 70 deletions(-) create mode 100644 src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java diff --git a/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java b/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java new file mode 100644 index 000000000..94da709f4 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java @@ -0,0 +1,66 @@ +/* + * 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.math.ode; + + +/** This interface represents a first order differential equations set + * with a main set of equations and an extension set. + * + *

+ * This interface is a simple extension on the {@link + * FirstOrderDifferentialEquations} that allows to identify which part + * of a complete set of differential equations correspond to the main + * set and which part correspond to the extension set. + *

+ *

+ * One typical use case is the computation of Jacobians. The main + * set of equations correspond to the raw ode, and we add to this set + * another bunch of equations which represent the jacobians of the + * main set. In that case, we want the integrator to use only + * the main set to estimate the errors and hence the step sizes. It should + * not use the additional equations in this computation. If the + * complete ode implements this interface, the {@link FirstOrderIntegrator + * integrator} will be able to know where the main set ends and where the + * extended set begins. + *

+ *

+ * We consider that the main set always corresponds to the first equations + * and the extended set to the last equations. + *

+ * + * @see FirstOrderDifferentialEquations + * + * @version $Revision$ $Date$ + * @since 2.2 + */ + +public interface ExtendedFirstOrderDifferentialEquations extends FirstOrderDifferentialEquations { + + /** Return the dimension of the main set of equations. + *

+ * The main set of equations represent the first part of an ODE state. + * The error estimations and adaptive step size computation should be + * done on this first part only, not on the final part of the state + * which represent an extension set of equations which are considered + * secondary. + *

+ * @return dimension of the main set of equations, must be lesser than or + * equal to the {@link #getDimension() total dimension} + */ + int getMainSetDimension(); + +} \ No newline at end of file diff --git a/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java b/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java index 5ec05c675..836b897ea 100644 --- a/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java @@ -378,7 +378,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { } /** Wrapper for differential equations, ensuring start evaluations are counted. */ - private class CountingDifferentialEquations implements FirstOrderDifferentialEquations { + private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations { /** Dimension of the problem. */ private final int dimension; @@ -400,6 +400,11 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { public int getDimension() { return dimension; } + + /** {@inheritDoc} */ + public int getMainSetDimension() { + return mainSetDimension; + } } } diff --git a/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java b/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java index f8b6c2cfc..15a8d8712 100644 --- a/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java +++ b/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java @@ -27,6 +27,7 @@ import java.util.Collection; import org.apache.commons.math.MathRuntimeException; import org.apache.commons.math.MaxEvaluationsExceededException; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; @@ -354,7 +355,7 @@ public class FirstOrderIntegratorWithJacobians { } /** Wrapper class used to map state and jacobians into compound state. */ - private class MappingWrapper implements FirstOrderDifferentialEquations { + private class MappingWrapper implements ExtendedFirstOrderDifferentialEquations { /** Current state. */ private final double[] y; @@ -388,6 +389,11 @@ public class FirstOrderIntegratorWithJacobians { return n * (1 + n + k); } + /** {@inheritDoc} */ + public int getMainSetDimension() { + return ode.getDimension(); + } + /** {@inheritDoc} */ public void computeDerivatives(final double t, final double[] z, final double[] zDot) throws DerivativeException { @@ -442,8 +448,7 @@ public class FirstOrderIntegratorWithJacobians { } /** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */ - private class FiniteDifferencesWrapper - implements ODEWithJacobians { + private class FiniteDifferencesWrapper implements ODEWithJacobians { /** Raw ODE without jacobians computation. */ private final ParameterizedODE ode; diff --git a/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java b/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java index a3ad45477..6c0020945 100644 --- a/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java +++ b/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java @@ -29,8 +29,7 @@ import org.apache.commons.math.ode.FirstOrderDifferentialEquations; * @since 2.1 */ -public interface ParameterizedODE - extends FirstOrderDifferentialEquations { +public interface ParameterizedODE extends FirstOrderDifferentialEquations { /** Get the number of parameters. * @return number of parameters diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java index 6ee9bd56c..2e9ccafd3 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java @@ -235,7 +235,7 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator { // evaluate error using the last term of the Taylor expansion error = 0; - for (int i = 0; i < y0.length; ++i) { + for (int i = 0; i < mainSetDimension; ++i) { final double yScale = Math.abs(y[i]); final double tol = (vecAbsoluteTolerance == null) ? (scalAbsoluteTolerance + scalRelativeTolerance * yScale) : @@ -243,7 +243,7 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator { final double ratio = nordsieck.getEntry(lastRow, i) / tol; error += ratio * ratio; } - error = Math.sqrt(error / y0.length); + error = Math.sqrt(error / mainSetDimension); if (error <= 1.0) { diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java index e0e2f0d2c..684a2e2c1 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java @@ -419,7 +419,7 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator { } /** - * End visiting te Nordsieck vector. + * End visiting the Nordsieck vector. *

The correction is used to control stepsize. So its amplitude is * considered to be an error, which must be normalized according to * error control settings. If the normalized value is greater than 1, @@ -432,15 +432,17 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator { double error = 0; for (int i = 0; i < after.length; ++i) { after[i] += previous[i] + scaled[i]; - final double yScale = Math.max(Math.abs(previous[i]), Math.abs(after[i])); - final double tol = (vecAbsoluteTolerance == null) ? - (scalAbsoluteTolerance + scalRelativeTolerance * yScale) : - (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale); - final double ratio = (after[i] - before[i]) / tol; - error += ratio * ratio; + if (i < mainSetDimension) { + final double yScale = Math.max(Math.abs(previous[i]), Math.abs(after[i])); + final double tol = (vecAbsoluteTolerance == null) ? + (scalAbsoluteTolerance + scalRelativeTolerance * yScale) : + (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale); + final double ratio = (after[i] - before[i]) / tol; + error += ratio * ratio; + } } - return Math.sqrt(error / after.length); + return Math.sqrt(error / mainSetDimension); } } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java index 9811c1497..65348069e 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java @@ -19,6 +19,7 @@ package org.apache.commons.math.ode.nonstiff; import org.apache.commons.math.ode.AbstractIntegrator; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.util.LocalizedFormats; @@ -36,14 +37,22 @@ import org.apache.commons.math.util.LocalizedFormats; * where absTol_i is the absolute tolerance for component i of the * state vector and relTol_i is the relative tolerance for the same * component. The user can also use only two scalar values absTol and - * relTol which will be used for all components.

+ * relTol which will be used for all components. + *

+ * + *

If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations + * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, + * then only the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension() + * main set} part of the state vector is used for stepsize control, not the complete + * state vector. + *

* *

If the estimated error for ym+1 is such that * 

  * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
  *
* - * (where n is the state vector dimension) then the step is accepted, + * (where n is the main set dimension) then the step is accepted, * otherwise the step is rejected and a new attempt is made with a new * stepsize.

* @@ -76,6 +85,9 @@ public abstract class AdaptiveStepsizeIntegrator /** Maximal step. */ private final double maxStep; + /** Main set dimension. */ + protected int mainSetDimension; + /** Build an integrator with the given stepsize bounds. * The default step handler does nothing. * @param name name of the method @@ -171,14 +183,20 @@ public abstract class AdaptiveStepsizeIntegrator super.sanityChecks(equations, t0, y0, t, y); - if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length)) { - throw new IntegratorException( - LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, y0.length, vecAbsoluteTolerance.length); + if (equations instanceof ExtendedFirstOrderDifferentialEquations) { + mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension(); + } else { + mainSetDimension = equations.getDimension(); } - if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != y0.length)) { + if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) { throw new IntegratorException( - LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, y0.length, vecRelativeTolerance.length); + LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecAbsoluteTolerance.length); + } + + if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension)) { + throw new IntegratorException( + LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecRelativeTolerance.length); } } @@ -187,7 +205,7 @@ public abstract class AdaptiveStepsizeIntegrator * @param equations differential equations set * @param forward forward integration indicator * @param order order of the method - * @param scale scaling vector for the state vector + * @param scale scaling vector for the state vector (can be shorter than state vector) * @param t0 start time * @param y0 state vector at t0 * @param yDot0 first time derivative of y0 @@ -213,7 +231,7 @@ public abstract class AdaptiveStepsizeIntegrator double ratio; double yOnScale2 = 0; double yDotOnScale2 = 0; - for (int j = 0; j < y0.length; ++j) { + for (int j = 0; j < scale.length; ++j) { ratio = y0[j] / scale[j]; yOnScale2 += ratio * ratio; ratio = yDot0[j] / scale[j]; @@ -234,7 +252,7 @@ public abstract class AdaptiveStepsizeIntegrator // estimate the second derivative of the solution double yDDotOnScale = 0; - for (int j = 0; j < y0.length; ++j) { + for (int j = 0; j < scale.length; ++j) { ratio = (yDot1[j] - yDot0[j]) / scale[j]; yDDotOnScale += ratio * ratio; } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java index 7ed175a6a..49d161b80 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java @@ -135,7 +135,7 @@ public class DormandPrince54Integrator extends EmbeddedRungeKuttaIntegrator { double error = 0; - for (int j = 0; j < y0.length; ++j) { + for (int j = 0; j < mainSetDimension; ++j) { final double errSum = E1 * yDotK[0][j] + E3 * yDotK[2][j] + E4 * yDotK[3][j] + E5 * yDotK[4][j] + E6 * yDotK[5][j] + E7 * yDotK[6][j]; @@ -149,7 +149,7 @@ public class DormandPrince54Integrator extends EmbeddedRungeKuttaIntegrator { } - return Math.sqrt(error / y0.length); + return Math.sqrt(error / mainSetDimension); } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java index 2d6b17e29..ebcc9914d 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java @@ -249,7 +249,7 @@ public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator { double error1 = 0; double error2 = 0; - for (int j = 0; j < y0.length; ++j) { + for (int j = 0; j < mainSetDimension; ++j) { final double errSum1 = E1_01 * yDotK[0][j] + E1_06 * yDotK[5][j] + E1_07 * yDotK[6][j] + E1_08 * yDotK[7][j] + E1_09 * yDotK[8][j] + E1_10 * yDotK[9][j] + @@ -274,7 +274,7 @@ public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator { den = 1.0; } - return Math.abs(h) * error1 / Math.sqrt(y0.length * den); + return Math.abs(h) * error1 / Math.sqrt(mainSetDimension * den); } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java index e03be9ed0..2c0f3a8fe 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java @@ -242,7 +242,7 @@ public abstract class EmbeddedRungeKuttaIntegrator } if (firstTime) { - final double[] scale = new double[y0.length]; + final double[] scale = new double[mainSetDimension]; if (vecAbsoluteTolerance == null) { for (int i = 0; i < scale.length; ++i) { scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * Math.abs(y[i]); diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java index 90c1073d3..5d8045f3c 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java @@ -429,7 +429,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { /** Update scaling array. * @param y1 first state vector to use for scaling * @param y2 second state vector to use for scaling - * @param scale scaling array to update + * @param scale scaling array to update (can be shorter than state) */ private void rescale(final double[] y1, final double[] y2, final double[] scale) { if (vecAbsoluteTolerance == null) { @@ -451,7 +451,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { * @param y0 initial value of the state vector at t0 * @param step global step * @param k iteration number (from 0 to sequence.length - 1) - * @param scale scaling array + * @param scale scaling array (can be shorter than state) * @param f placeholder where to put the state vector derivatives at each substep * (element 0 already contains initial derivative) * @param yMiddle placeholder where to put the state vector at the middle of the step @@ -500,12 +500,12 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { // stability check if (performTest && (j <= maxChecks) && (k < maxIter)) { double initialNorm = 0.0; - for (int l = 0; l < y0.length; ++l) { + for (int l = 0; l < scale.length; ++l) { final double ratio = f[0][l] / scale[l]; initialNorm += ratio * ratio; } double deltaNorm = 0.0; - for (int l = 0; l < y0.length; ++l) { + for (int l = 0; l < scale.length; ++l) { final double ratio = (f[j+1][l] - f[0][l]) / scale[l]; deltaNorm += ratio * ratio; } @@ -607,7 +607,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { } // initial scaling - final double[] scale = new double[y0.length]; + final double[] scale = new double[mainSetDimension]; rescale(y, y, scale); // initial order selection @@ -709,11 +709,11 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { // estimate the error at the end of the step. error = 0; - for (int j = 0; j < y0.length; ++j) { + for (int j = 0; j < mainSetDimension; ++j) { final double e = Math.abs(y1[j] - y1Diag[0][j]) / scale[j]; error += e * e; } - error = Math.sqrt(error / y0.length); + error = Math.sqrt(error / mainSetDimension); if ((error > 1.0e15) || ((k > 1) && (error > maxError))) { // error is too big, we reduce the global step diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java index 737131c2b..e91002b81 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java @@ -297,11 +297,11 @@ class GraggBulirschStoerStepInterpolator public double estimateError(final double[] scale) { double error = 0; if (currentDegree >= 5) { - for (int i = 0; i < currentState.length; ++i) { + for (int i = 0; i < scale.length; ++i) { final double e = polynoms[currentDegree][i] / scale[i]; error += e * e; } - error = Math.sqrt(error / currentState.length) * errfac[currentDegree-5]; + error = Math.sqrt(error / scale.length) * errfac[currentDegree - 5]; } return error; } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java index c78afeadd..8203d554f 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java @@ -108,7 +108,7 @@ public class HighamHall54Integrator extends EmbeddedRungeKuttaIntegrator { double error = 0; - for (int j = 0; j < y0.length; ++j) { + for (int j = 0; j < mainSetDimension; ++j) { double errSum = STATIC_E[0] * yDotK[0][j]; for (int l = 1; l < STATIC_E.length; ++l) { errSum += STATIC_E[l] * yDotK[l][j]; @@ -123,7 +123,7 @@ public class HighamHall54Integrator extends EmbeddedRungeKuttaIntegrator { } - return Math.sqrt(error / y0.length); + return Math.sqrt(error / mainSetDimension); } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 441ed69f1..2a1336f33 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,12 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Added a feature allowing error estimation to be computed only on a subset of + Ordinary Differential Equations, considered as the main set, the remaining equations + being considered only as an extension set that should not influence the ODE integration + algorithm + Fixed bug in precondition check (method "setMicrosphereElements"). diff --git a/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java b/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java index 436ac9bfd..efbfb8f49 100644 --- a/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java +++ b/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java @@ -40,7 +40,7 @@ public class FirstOrderIntegratorWithJacobiansTest { // Solving Ordinary Differential Equations I (Nonstiff problems), // the curves dy/dp = g(b) are in figure 6.5 FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); double hP = 1.0e-12; SummaryStatistics residualsP0 = new SummaryStatistics(); SummaryStatistics residualsP1 = new SummaryStatistics(); @@ -64,7 +64,7 @@ public class FirstOrderIntegratorWithJacobiansTest { public void testHighAccuracyExternalDifferentiation() throws IntegratorException, DerivativeException { FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); double hP = 1.0e-12; SummaryStatistics residualsP0 = new SummaryStatistics(); SummaryStatistics residualsP1 = new SummaryStatistics(); @@ -92,7 +92,7 @@ public class FirstOrderIntegratorWithJacobiansTest { public void testInternalDifferentiation() throws IntegratorException, DerivativeException { FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); double hP = 1.0e-12; SummaryStatistics residualsP0 = new SummaryStatistics(); SummaryStatistics residualsP1 = new SummaryStatistics(); @@ -109,23 +109,23 @@ public class FirstOrderIntegratorWithJacobiansTest { extInt.setMaxEvaluations(5000); extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP); Assert.assertEquals(5000, extInt.getMaxEvaluations()); - Assert.assertTrue(extInt.getEvaluations() > 2000); - Assert.assertTrue(extInt.getEvaluations() < 2500); + Assert.assertTrue(extInt.getEvaluations() > 1500); + Assert.assertTrue(extInt.getEvaluations() < 2100); Assert.assertEquals(4 * integ.getEvaluations(), extInt.getEvaluations()); residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0()); residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1()); } - Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.006); - Assert.assertTrue(residualsP0.getStandardDeviation() < 0.0009); - Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.009); - Assert.assertTrue(residualsP1.getStandardDeviation() < 0.0014); + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01); } @Test public void testAnalyticalDifferentiation() throws IntegratorException, DerivativeException { FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); SummaryStatistics residualsP0 = new SummaryStatistics(); SummaryStatistics residualsP1 = new SummaryStatistics(); for (double b = 2.88; b < 3.08; b += 0.001) { @@ -139,22 +139,22 @@ public class FirstOrderIntegratorWithJacobiansTest { extInt.setMaxEvaluations(5000); extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP); Assert.assertEquals(5000, extInt.getMaxEvaluations()); - Assert.assertTrue(extInt.getEvaluations() > 510); - Assert.assertTrue(extInt.getEvaluations() < 610); + Assert.assertTrue(extInt.getEvaluations() > 350); + Assert.assertTrue(extInt.getEvaluations() < 510); Assert.assertEquals(integ.getEvaluations(), extInt.getEvaluations()); residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0()); residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1()); } - Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.004); - Assert.assertTrue(residualsP0.getStandardDeviation() < 0.0008); - Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.005); - Assert.assertTrue(residualsP1.getStandardDeviation() < 0.0010); + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01); } @Test public void testFinalResult() throws IntegratorException, DerivativeException { FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); double[] y = new double[] { 0.0, 1.0 }; Circle circle = new Circle(y, 1.0, 1.0, 0.1); double[][] dydy0 = new double[2][2]; @@ -164,16 +164,16 @@ public class FirstOrderIntegratorWithJacobiansTest { new FirstOrderIntegratorWithJacobians(integ, circle); extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp); for (int i = 0; i < y.length; ++i) { - Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-10); + Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9); } for (int i = 0; i < dydy0.length; ++i) { for (int j = 0; j < dydy0[i].length; ++j) { - Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-10); + Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9); } } for (int i = 0; i < dydp.length; ++i) { for (int j = 0; j < dydp[i].length; ++j) { - Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-8); + Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-7); } } } @@ -181,7 +181,7 @@ public class FirstOrderIntegratorWithJacobiansTest { @Test public void testStepHandlerResult() throws IntegratorException, DerivativeException { FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); double[] y = new double[] { 0.0, 1.0 }; final Circle circle = new Circle(y, 1.0, 1.0, 0.1); double[][] dydy0 = new double[2][2]; @@ -207,16 +207,16 @@ public class FirstOrderIntegratorWithJacobiansTest { Assert.assertEquals(interpolator.getPreviousTime(), extInt.getCurrentStepStart(), 1.0e-10); Assert.assertTrue(extInt.getCurrentSignedStepsize() < 0.5); for (int i = 0; i < y.length; ++i) { - Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-10); + Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9); } for (int i = 0; i < dydy0.length; ++i) { for (int j = 0; j < dydy0[i].length; ++j) { - Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-10); + Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9); } } for (int i = 0; i < dydp.length; ++i) { for (int j = 0; j < dydp[i].length; ++j) { - Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-8); + Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 3.0e-8); } } @@ -225,16 +225,16 @@ public class FirstOrderIntegratorWithJacobiansTest { double[][] dydpDot = interpolator.getInterpolatedDyDpDot(); for (int i = 0; i < yDot.length; ++i) { - Assert.assertEquals(circle.exactYDot(t)[i], yDot[i], 1.0e-11); + Assert.assertEquals(circle.exactYDot(t)[i], yDot[i], 1.0e-10); } for (int i = 0; i < dydy0Dot.length; ++i) { for (int j = 0; j < dydy0Dot[i].length; ++j) { - Assert.assertEquals(circle.exactDyDy0Dot(t)[i][j], dydy0Dot[i][j], 1.0e-11); + Assert.assertEquals(circle.exactDyDy0Dot(t)[i][j], dydy0Dot[i][j], 1.0e-10); } } for (int i = 0; i < dydpDot.length; ++i) { for (int j = 0; j < dydpDot[i].length; ++j) { - Assert.assertEquals(circle.exactDyDpDot(t)[i][j], dydpDot[i][j], 1.0e-9); + Assert.assertEquals(circle.exactDyDpDot(t)[i][j], dydpDot[i][j], 3.0e-9); } } } @@ -245,7 +245,7 @@ public class FirstOrderIntegratorWithJacobiansTest { @Test public void testEventHandler() throws IntegratorException, DerivativeException { FirstOrderIntegrator integ = - new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10); + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); double[] y = new double[] { 0.0, 1.0 }; final Circle circle = new Circle(y, 1.0, 1.0, 0.1); double[][] dydy0 = new double[2][2];