commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r967007 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/ main/java/org/apache/commons/math/ode/jacobians/ main/java/org/apache/commons/math/ode/nonstiff/ site/xdoc/ test/java/org/apache/commons/math/ode/jacobians/
Date Fri, 23 Jul 2010 08:59:38 GMT
Author: luc
Date: Fri Jul 23 08:59:37 2010
New Revision: 967007

URL: http://svn.apache.org/viewvc?rev=967007&view=rev
Log:
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

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
  (with props)
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java?rev=967007&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
Fri Jul 23 08:59:37 2010
@@ -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.
+ *
+ * <p>
+ * 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.
+ * </p>
+ * <p>
+ * 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 <em>only</em>
+ * the main set to estimate the errors and hence the step sizes. It should
+ * <em>not</em> 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.
+ * </p>
+ * <p>
+ * We consider that the main set always corresponds to the first equations
+ * and the extended set to the last equations.
+ * </p>
+ *
+ * @see FirstOrderDifferentialEquations
+ *
+ * @version $Revision$ $Date$
+ * @since 2.2
+ */
+
+public interface ExtendedFirstOrderDifferentialEquations extends FirstOrderDifferentialEquations
{
+
+    /** Return the dimension of the main set of equations.
+     * <p>
+     * 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.
+     * </p>
+     * @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

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java
Fri Jul 23 08:59:37 2010
@@ -378,7 +378,7 @@ public abstract class MultistepIntegrato
     }
 
     /** 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 MultistepIntegrato
         public int getDimension() {
             return dimension;
         }
+
+        /** {@inheritDoc} */
+        public int getMainSetDimension() {
+            return mainSetDimension;
+        }
     }
 
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
Fri Jul 23 08:59:37 2010
@@ -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 FirstOrderIntegratorWithJac
     }
 
     /** 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;
@@ -389,6 +390,11 @@ public class FirstOrderIntegratorWithJac
         }
 
         /** {@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 FirstOrderIntegratorWithJac
     }
 
     /** 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;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java
Fri Jul 23 08:59:37 2010
@@ -29,8 +29,7 @@ import org.apache.commons.math.ode.First
  * @since 2.1
  */
 
-public interface ParameterizedODE
-    extends FirstOrderDifferentialEquations {
+public interface ParameterizedODE extends FirstOrderDifferentialEquations {
 
     /** Get the number of parameters.
      * @return number of parameters

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
Fri Jul 23 08:59:37 2010
@@ -235,7 +235,7 @@ public class AdamsBashforthIntegrator ex
 
                 // 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 ex
                     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) {
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
Fri Jul 23 08:59:37 2010
@@ -419,7 +419,7 @@ public class AdamsMoultonIntegrator exte
         }
 
         /**
-         * End visiting te Nordsieck vector.
+         * End visiting the Nordsieck vector.
          * <p>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 exte
             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);
 
         }
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
Fri Jul 23 08:59:37 2010
@@ -19,6 +19,7 @@ package org.apache.commons.math.ode.nons
 
 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.Loca
  * 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.</p>
+ * relTol which will be used for all components.
+ * </p>
+ *
+ * <p>If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations
+ * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE},
+ * then <em>only</em> the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension()
+ * main set} part of the state vector is used for stepsize control, not the complete
+ * state vector.
+ * </p>
  *
  * <p>If the estimated error for ym+1 is such that
  * <pre>
  * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
  * </pre>
  *
- * (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.</p>
  *
@@ -76,6 +85,9 @@ public abstract class AdaptiveStepsizeIn
     /** 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 AdaptiveStepsizeIn
 
       super.sanityChecks(equations, t0, y0, t, y);
 
-      if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length))
{
+      if (equations instanceof ExtendedFirstOrderDifferentialEquations) {
+          mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension();
+      } else {
+          mainSetDimension = equations.getDimension();
+      }
+
+      if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension))
{
           throw new IntegratorException(
-                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, y0.length, vecAbsoluteTolerance.length);
+                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecAbsoluteTolerance.length);
       }
 
-      if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != y0.length))
{
+      if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension))
{
           throw new IntegratorException(
-                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, y0.length, vecRelativeTolerance.length);
+                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecRelativeTolerance.length);
       }
 
   }
@@ -187,7 +205,7 @@ public abstract class AdaptiveStepsizeIn
    * @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 AdaptiveStepsizeIn
     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 AdaptiveStepsizeIn
 
     // 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;
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.java
Fri Jul 23 08:59:37 2010
@@ -135,7 +135,7 @@ public class DormandPrince54Integrator e
 
     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 e
 
     }
 
-    return Math.sqrt(error / y0.length);
+    return Math.sqrt(error / mainSetDimension);
 
   }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.java
Fri Jul 23 08:59:37 2010
@@ -249,7 +249,7 @@ public class DormandPrince853Integrator 
     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 
       den = 1.0;
     }
 
-    return Math.abs(h) * error1 / Math.sqrt(y0.length * den);
+    return Math.abs(h) * error1 / Math.sqrt(mainSetDimension * den);
 
   }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
Fri Jul 23 08:59:37 2010
@@ -242,7 +242,7 @@ public abstract class EmbeddedRungeKutta
         }
 
         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]);

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
Fri Jul 23 08:59:37 2010
@@ -429,7 +429,7 @@ public class GraggBulirschStoerIntegrato
   /** 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 GraggBulirschStoerIntegrato
    * @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 GraggBulirschStoerIntegrato
       // 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 GraggBulirschStoerIntegrato
     }
 
     // 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 GraggBulirschStoerIntegrato
 
             // 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

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
Fri Jul 23 08:59:37 2010
@@ -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;
   }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.java
Fri Jul 23 08:59:37 2010
@@ -108,7 +108,7 @@ public class HighamHall54Integrator exte
 
     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 exte
 
     }
 
-    return Math.sqrt(error / y0.length);
+    return Math.sqrt(error / mainSetDimension);
 
   }
 

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Jul 23 08:59:37 2010
@@ -52,6 +52,12 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="2.2" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="MATH-388">
+        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 
+      </action>
       <action dev="erans" type="fix" issue="MATH-382">
         Fixed bug in precondition check (method "setMicrosphereElements").
       </action>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java?rev=967007&r1=967006&r2=967007&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
Fri Jul 23 08:59:37 2010
@@ -40,7 +40,7 @@ public class FirstOrderIntegratorWithJac
         // 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 FirstOrderIntegratorWithJac
     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 FirstOrderIntegratorWithJac
     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 FirstOrderIntegratorWithJac
             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 FirstOrderIntegratorWithJac
             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 FirstOrderIntegratorWithJac
             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 FirstOrderIntegratorWithJac
     @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 FirstOrderIntegratorWithJac
                 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 FirstOrderIntegratorWithJac
                 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 FirstOrderIntegratorWithJac
     @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];



Mime
View raw message