commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1175409 [2/2] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/exception/util/ main/java/org/apache/commons/math/ode/ main/java/org/apache/commons/math/ode/nonstiff/ main/resources/META-INF/localization/ site/xdoc/ te...
Date Sun, 25 Sep 2011 15:04:40 GMT
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java Sun Sep 25 15:04:39 2011
@@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nons
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
-import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.MultistepIntegrator;
 
 
@@ -86,7 +86,7 @@ public abstract class AdamsIntegrator ex
 
     /** {@inheritDoc} */
     @Override
-    public abstract double integrate(final FirstOrderDifferentialEquations equations,
+    public abstract double integrate(final ExpandableFirstOrderDifferentialEquations equations,
                                      final double t0, final double[] y0,
                                      final double t, final double[] y)
         throws MathIllegalStateException, MathIllegalArgumentException;

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=1175409&r1=1175408&r2=1175409&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 Sun Sep 25 15:04:39 2011
@@ -23,7 +23,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.MathIllegalStateException;
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
 import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
-import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
@@ -206,24 +206,29 @@ public class AdamsMoultonIntegrator exte
 
     /** {@inheritDoc} */
     @Override
-    public double integrate(final FirstOrderDifferentialEquations equations,
-                            final double t0, final double[] y0,
-                            final double t, final double[] y)
+    public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
+                            final double t0, final double[] z0,
+                            final double t, final double[] z)
         throws MathIllegalStateException, MathIllegalArgumentException {
 
-        final int n = y0.length;
-        sanityChecks(equations, t0, y0, t, y);
+        sanityChecks(equations, t0, z0, t, z);
         setEquations(equations);
         resetEvaluations();
         final boolean forward = t > t0;
 
         // initialize working arrays
+        final int totalDim = equations.getDimension();
+        final int mainDim  = equations.getMainSetDimension();
+        final double[] y0  = new double[totalDim];
+        final double[] y   = new double[totalDim];
+        System.arraycopy(z0, 0, y0, 0, mainDim);
+        System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
         if (y != y0) {
-            System.arraycopy(y0, 0, y, 0, n);
+            System.arraycopy(y0, 0, y, 0, totalDim);
         }
-        final double[] yDot = new double[y0.length];
-        final double[] yTmp = new double[y0.length];
-        final double[] predictedScaled = new double[y0.length];
+        final double[] yDot = new double[totalDim];
+        final double[] yTmp = new double[totalDim];
+        final double[] predictedScaled = new double[totalDim];
         Array2DRowRealMatrix nordsieckTmp = null;
 
         // set up two interpolators sharing the integrator arrays
@@ -290,7 +295,7 @@ public class AdamsMoultonIntegrator exte
             updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
 
             // discrete events handling
-            System.arraycopy(yTmp, 0, y, 0, n);
+            System.arraycopy(yTmp, 0, y, 0, totalDim);
             interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
             interpolator.storeTime(stepStart);
             interpolator.shift();
@@ -330,6 +335,10 @@ public class AdamsMoultonIntegrator exte
 
         } while (!isLastStep);
 
+        // dispatch result between main and additional states
+        System.arraycopy(y, 0, z, 0, z.length);
+        equations.setCurrentAdditionalState(y);
+
         final double stopTime  = stepStart;
         stepStart = Double.NaN;
         stepSize  = Double.NaN;

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=1175409&r1=1175408&r2=1175409&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 Sun Sep 25 15:04:39 2011
@@ -23,7 +23,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.NumberIsTooSmallException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
 import org.apache.commons.math.ode.AbstractIntegrator;
-import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.util.FastMath;
 
@@ -42,12 +42,11 @@ import org.apache.commons.math.util.Fast
  * component. The user can also use only two scalar values absTol and
  * 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>
+ * If the Ordinary Differential Equations is an {@link ExpandableFirstOrderDifferentialEquations
+ * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, then
+ * <em>only</em> the {@link ExpandableFirstOrderDifferentialEquations#getMainSet() main 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
@@ -224,18 +223,14 @@ public abstract class AdaptiveStepsizeIn
    * @exception NumberIsTooSmallException if integration span is too small
    */
   @Override
-  protected void sanityChecks(final FirstOrderDifferentialEquations equations,
+  protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations equations,
                               final double t0, final double[] y0,
                               final double t, final double[] y)
       throws DimensionMismatchException, NumberIsTooSmallException {
 
       super.sanityChecks(equations, t0, y0, t, y);
 
-      if (equations instanceof ExtendedFirstOrderDifferentialEquations) {
-          mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension();
-      } else {
-          mainSetDimension = equations.getDimension();
-      }
+      mainSetDimension = equations.getMainSetDimension();
 
       if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
           throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
@@ -248,7 +243,6 @@ public abstract class AdaptiveStepsizeIn
   }
 
   /** Initialize the integration step.
-   * @param equations differential equations set
    * @param forward forward integration indicator
    * @param order order of the method
    * @param scale scaling vector for the state vector (can be shorter than state vector)
@@ -259,8 +253,7 @@ public abstract class AdaptiveStepsizeIn
    * @param yDot1 work array for the first time derivative of y1
    * @return first integration step
    */
-  public double initializeStep(final FirstOrderDifferentialEquations equations,
-                               final boolean forward, final int order, final double[] scale,
+  public double initializeStep(final boolean forward, final int order, final double[] scale,
                                final double t0, final double[] y0, final double[] yDot0,
                                final double[] y1, final double[] yDot1) {
 
@@ -356,7 +349,7 @@ public abstract class AdaptiveStepsizeIn
   }
 
   /** {@inheritDoc} */
-  public abstract double integrate (FirstOrderDifferentialEquations equations,
+  public abstract double integrate (ExpandableFirstOrderDifferentialEquations equations,
                                     double t0, double[] y0,
                                     double t, double[] y)
     throws MathIllegalStateException, MathIllegalArgumentException;

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=1175409&r1=1175408&r2=1175409&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 Sun Sep 25 15:04:39 2011
@@ -19,7 +19,7 @@ package org.apache.commons.math.ode.nons
 
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
-import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
 
@@ -189,24 +189,30 @@ public abstract class EmbeddedRungeKutta
 
   /** {@inheritDoc} */
   @Override
-  public double integrate(final FirstOrderDifferentialEquations equations,
-                          final double t0, final double[] y0,
-                          final double t, final double[] y)
+  public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
+                          final double t0, final double[] z0,
+                          final double t, final double[] z)
       throws MathIllegalStateException, MathIllegalArgumentException {
 
-    sanityChecks(equations, t0, y0, t, y);
+    sanityChecks(equations, t0, z0, t, z);
     setEquations(equations);
     resetEvaluations();
     final boolean forward = t > t0;
 
     // create some internal working arrays
+    final int totalDim = equations.getDimension();
+    final int mainDim  = equations.getMainSetDimension();
+    final double[] y0  = new double[totalDim];
+    final double[] y   = new double[totalDim];
+    System.arraycopy(z0, 0, y0, 0, mainDim);
+    System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
     final int stages = c.length + 1;
     if (y != y0) {
-      System.arraycopy(y0, 0, y, 0, y0.length);
+      System.arraycopy(y0, 0, y, 0, totalDim);
     }
-    final double[][] yDotK = new double[stages][y0.length];
-    final double[] yTmp    = new double[y0.length];
-    final double[] yDotTmp = new double[y0.length];
+    final double[][] yDotK = new double[stages][totalDim];
+    final double[] yTmp    = new double[totalDim];
+    final double[] yDotTmp = new double[totalDim];
 
     // set up an interpolator sharing the integrator arrays
     final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
@@ -248,7 +254,7 @@ public abstract class EmbeddedRungeKutta
                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
               }
           }
-          hNew = initializeStep(equations, forward, getOrder(), scale,
+          hNew = initializeStep(forward, getOrder(), scale,
                                 stepStart, y, yDotK[0], yTmp, yDotK[1]);
           firstTime = false;
         }
@@ -325,6 +331,10 @@ public abstract class EmbeddedRungeKutta
 
     } while (!isLastStep);
 
+    // dispatch result between main and additional states
+    System.arraycopy(y, 0, z, 0, z.length);
+    equations.setCurrentAdditionalState(y);
+
     final double stopTime = stepStart;
     resetInternalState();
     return stopTime;

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=1175409&r1=1175408&r2=1175409&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 Sun Sep 25 15:04:39 2011
@@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nons
 import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
-import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.events.EventHandler;
 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
@@ -541,26 +541,32 @@ public class GraggBulirschStoerIntegrato
 
   /** {@inheritDoc} */
   @Override
-  public double integrate(final FirstOrderDifferentialEquations equations,
-                          final double t0, final double[] y0, final double t, final double[] y)
+  public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
+                          final double t0, final double[] z0, final double t, final double[] z)
       throws MathIllegalStateException, MathIllegalArgumentException {
 
-    sanityChecks(equations, t0, y0, t, y);
+    sanityChecks(equations, t0, z0, t, z);
     setEquations(equations);
     resetEvaluations();
     final boolean forward = t > t0;
 
     // create some internal working arrays
-    final double[] yDot0   = new double[y0.length];
-    final double[] y1      = new double[y0.length];
-    final double[] yTmp    = new double[y0.length];
-    final double[] yTmpDot = new double[y0.length];
+    final int totalDim = equations.getDimension();
+    final int mainDim  = equations.getMainSetDimension();
+    final double[] y0 = new double[totalDim];
+    final double[] y  = new double[totalDim];
+    System.arraycopy(z0, 0, y0, 0, mainDim);
+    System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
+    final double[] yDot0   = new double[totalDim];
+    final double[] y1      = new double[totalDim];
+    final double[] yTmp    = new double[totalDim];
+    final double[] yTmpDot = new double[totalDim];
 
     final double[][] diagonal = new double[sequence.length-1][];
     final double[][] y1Diag = new double[sequence.length-1][];
     for (int k = 0; k < sequence.length-1; ++k) {
-      diagonal[k] = new double[y0.length];
-      y1Diag[k] = new double[y0.length];
+      diagonal[k] = new double[totalDim];
+      y1Diag[k] = new double[totalDim];
     }
 
     final double[][][] fk  = new double[sequence.length][][];
@@ -631,8 +637,7 @@ public class GraggBulirschStoerIntegrato
         }
 
         if (firstTime) {
-          hNew = initializeStep(equations, forward,
-                                2 * targetIter + 1, scale,
+          hNew = initializeStep(forward, 2 * targetIter + 1, scale,
                                 stepStart, y, yDot0, yTmp, yTmpDot);
         }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java Sun Sep 25 15:04:39 2011
@@ -21,7 +21,7 @@ package org.apache.commons.math.ode.nons
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
 import org.apache.commons.math.ode.AbstractIntegrator;
-import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
 
@@ -90,17 +90,23 @@ public abstract class RungeKuttaIntegrat
   }
 
   /** {@inheritDoc} */
-  public double integrate(final FirstOrderDifferentialEquations equations,
-                          final double t0, final double[] y0,
-                          final double t, final double[] y)
+  public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
+                          final double t0, final double[] z0,
+                          final double t, final double[] z)
       throws MathIllegalStateException, MathIllegalArgumentException {
 
-    sanityChecks(equations, t0, y0, t, y);
+    sanityChecks(equations, t0, z0, t, z);
     setEquations(equations);
     resetEvaluations();
     final boolean forward = t > t0;
 
     // create some internal working arrays
+    final int totalDim = equations.getDimension();
+    final int mainDim  = equations.getMainSetDimension();
+    final double[] y0  = new double[totalDim];
+    final double[] y   = new double[totalDim];
+    System.arraycopy(z0, 0, y0, 0, mainDim);
+    System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
     final int stages = c.length + 1;
     if (y != y0) {
       System.arraycopy(y0, 0, y, 0, y0.length);
@@ -179,6 +185,10 @@ public abstract class RungeKuttaIntegrat
 
     } while (!isLastStep);
 
+    // dispatch result between main and additional states
+    System.arraycopy(y, 0, z, 0, z.length);
+    equations.setCurrentAdditionalState(y);
+
     final double stopTime = stepStart;
     stepStart = Double.NaN;
     stepSize  = Double.NaN;

Modified: commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties?rev=1175409&r1=1175408&r2=1175409&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties (original)
+++ commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties Sun Sep 25 15:04:39 2011
@@ -279,6 +279,9 @@ UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JA
 UNABLE_TO_SOLVE_SINGULAR_PROBLEM = r\u00e9solution impossible : probl\u00e8me singulier
 UNBOUNDED_SOLUTION = solution non born\u00e9e
 UNKNOWN_MODE = mode {0} inconnu, modes connus : {1} ({2}), {3} ({4}), {5} ({6}), {7} ({8}), {9} ({10}) et {11} ({12})
+UNKNOWN_ADDITIONAL_EQUATION = \u00e9quation additionnelle inconnue
+UNKNOWN_PARAMETER = param\u00e8tre {0} inconnu
+UNMATCHED_ODE_IN_EXTENDED_SET = l''\u00e9quation diff\u00e9rentielle ne correspond pas \u00e0 l''\u00e9quation principale du jeu \u00e9tendu
 CANNOT_PARSE_AS_TYPE = cha\u00eene {0} non analysable (\u00e0 partir de la position {1}) en un objet de type {2}
 CANNOT_PARSE = cha\u00eene {0} non analysable (\u00e0 partir de la position {1})
 UNPARSEABLE_3D_VECTOR = vecteur 3D non analysable : "{0}"

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=1175409&r1=1175408&r2=1175409&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun Sep 25 15:04:39 2011
@@ -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="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-380" due-to="Pascal Parraud">
+        Completely revamped the computation of Jacobians in ODE. This computation is now
+        included in the mainstream class hierarchy, not in a separate package anymore,
+        and it allows adding other types of equations to a main ODE, not only variational
+        equations for Jacobians computation.
+      </action>
       <action dev="gregs" type="update" issue="MATH-607">
         SimpleRegression implements UpdatingMultipleLinearRegression interface.
       </action>

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java?rev=1175409&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java Sun Sep 25 15:04:39 2011
@@ -0,0 +1,598 @@
+/*
+ * 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;
+
+import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator;
+import org.apache.commons.math.stat.descriptive.SummaryStatistics;
+import org.apache.commons.math.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class JacobianMatricesTest {
+
+    @Test
+    public void testLowAccuracyExternalDifferentiation() {
+        // this test does not really test FirstOrderIntegratorWithJacobians,
+        // it only shows that WITHOUT this class, attempting to recover
+        // the jacobians from external differentiation on simple integration
+        // results with low accuracy gives very poor results. In fact,
+        // the curves dy/dp = g(b) when b varies from 2.88 to 3.08 are
+        // essentially noise.
+        // This test is taken from Hairer, Norsett and Wanner book
+        // 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, 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();
+        for (double b = 2.88; b < 3.08; b += 0.001) {
+            Brusselator brusselator = new Brusselator(b);
+            double[] y = { 1.3, b };
+            integ.integrate(brusselator, 0, y, 20.0, y);
+            double[] yP = { 1.3, b + hP };
+            integ.integrate(brusselator, 0, yP, 20.0, yP);
+            residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
+            residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
+        }
+        Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 500);
+        Assert.assertTrue(residualsP0.getStandardDeviation() > 30);
+        Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 700);
+        Assert.assertTrue(residualsP1.getStandardDeviation() > 40);
+    }
+
+    @Test
+    public void testHighAccuracyExternalDifferentiation() {
+        FirstOrderIntegrator integ =
+            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();
+        for (double b = 2.88; b < 3.08; b += 0.001) {
+            ParamBrusselator brusselator = new ParamBrusselator(b);
+            double[] y = { 1.3, b };
+            integ.integrate(brusselator, 0, y, 20.0, y);
+            double[] yP = { 1.3, b + hP };
+            brusselator.setParameter("b", b + hP);
+            integ.integrate(brusselator, 0, yP, 20.0, yP);
+            residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
+            residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
+        }
+        Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02);
+        Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03);
+        Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003);
+        Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004);
+        Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04);
+        Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
+        Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007);
+        Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008);
+    }
+
+    @Test
+    public void testInternalDifferentiation() {
+        ExpandableFirstOrderIntegrator integ =
+            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;
+        double hY = 1.0e-12;
+        SummaryStatistics residualsP0 = new SummaryStatistics();
+        SummaryStatistics residualsP1 = new SummaryStatistics();
+        for (double b = 2.88; b < 3.08; b += 0.001) {
+            ParamBrusselator brusselator = new ParamBrusselator(b);
+            brusselator.setParameter(ParamBrusselator.B, b);
+            double[] z = { 1.3, b };
+            double[][] dZdZ0 = new double[2][2];
+            double[]   dZdP  = new double[2];
+            ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(brusselator);
+            JacobianMatrices jacob = new JacobianMatrices(efode, brusselator);
+            jacob.setParameterizedODE(brusselator);
+            jacob.selectParameters(ParamBrusselator.B);
+            jacob.setMainStateSteps(new double[] { hY, hY });
+            jacob.setParameterStep(ParamBrusselator.B, hP);
+            jacob.setInitialParameterJacobian(ParamBrusselator.B, new double[] { 0.0, 1.0 });
+            integ.setMaxEvaluations(5000);
+            integ.integrate(efode, 0, z, 20.0, z);
+            jacob.getCurrentMainSetJacobian(dZdZ0);
+            jacob.getCurrentParameterJacobian(ParamBrusselator.B, dZdP);
+//            Assert.assertEquals(5000, integ.getMaxEvaluations());
+//            Assert.assertTrue(integ.getEvaluations() > 1500);
+//            Assert.assertTrue(integ.getEvaluations() < 2100);
+//            Assert.assertEquals(4 * integ.getEvaluations(), integ.getEvaluations());
+            residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
+            residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
+        }
+        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() {
+        ExpandableFirstOrderIntegrator integ =
+            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) {
+            Brusselator brusselator = new Brusselator(b);
+            double[] z = { 1.3, b };
+            double[][] dZdZ0 = new double[2][2];
+            double[]   dZdP  = new double[2];
+            ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(brusselator);
+            JacobianMatrices jacob = new JacobianMatrices(efode, brusselator);
+            jacob.setParameterJacobianProvider(brusselator);
+            jacob.selectParameters(Brusselator.B);
+            jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 });
+            integ.setMaxEvaluations(5000);
+            integ.integrate(efode, 0, z, 20.0, z);
+            jacob.getCurrentMainSetJacobian(dZdZ0);
+            jacob.getCurrentParameterJacobian(Brusselator.B, dZdP);
+//            Assert.assertEquals(5000, integ.getMaxEvaluations());
+//            Assert.assertTrue(integ.getEvaluations() > 350);
+//            Assert.assertTrue(integ.getEvaluations() < 510);
+            residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
+            residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
+        }
+        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() {
+
+        ExpandableFirstOrderIntegrator integ =
+            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);
+
+        ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(circle);
+        JacobianMatrices jacob = new JacobianMatrices(efode, circle);
+        jacob.setParameterJacobianProvider(circle);
+        jacob.selectParameters(Circle.CX, Circle.CY, Circle.OMEGA);
+        jacob.setInitialMainStateJacobian(circle.exactDyDy0(0));
+        jacob.setInitialParameterJacobian(Circle.CX, circle.exactDyDcx(0));
+        jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
+        jacob.setInitialParameterJacobian(Circle.OMEGA, circle.exactDyDom(0));
+        integ.setMaxEvaluations(5000);
+
+        double t = 18 * FastMath.PI;
+        integ.integrate(efode, 0, y, t, y);
+        for (int i = 0; i < y.length; ++i) {
+            Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
+        }
+
+        double[][] dydy0 = new double[2][2];
+        jacob.getCurrentMainSetJacobian(dydy0);
+        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-9);
+            }
+        }
+        double[] dydcx = new double[2];
+        jacob.getCurrentParameterJacobian(Circle.CX, dydcx);
+        for (int i = 0; i < dydcx.length; ++i) {
+            Assert.assertEquals(circle.exactDyDcx(t)[i], dydcx[i], 1.0e-7);
+        }
+        double[] dydcy = new double[2];
+        jacob.getCurrentParameterJacobian(Circle.CY, dydcy);
+        for (int i = 0; i < dydcy.length; ++i) {
+            Assert.assertEquals(circle.exactDyDcy(t)[i], dydcy[i], 1.0e-7);
+        }
+        double[] dydom = new double[2];
+        jacob.getCurrentParameterJacobian(Circle.OMEGA, dydom);
+        for (int i = 0; i < dydom.length; ++i) {
+            Assert.assertEquals(circle.exactDyDom(t)[i], dydom[i], 1.0e-7);
+        }
+    }
+
+    @Test
+    public void testParameterizable() {
+
+        ExpandableFirstOrderIntegrator integ =
+            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 };
+        ParameterizedCircle pcircle = new ParameterizedCircle(y, 1.0, 1.0, 0.1);
+//        pcircle.setParameter(ParameterizedCircle.CX, 1.0);
+//        pcircle.setParameter(ParameterizedCircle.CY, 1.0);
+//        pcircle.setParameter(ParameterizedCircle.OMEGA, 0.1);
+
+        ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(pcircle);
+
+        double hP = 1.0e-12;
+        double hY = 1.0e-12;
+
+        JacobianMatrices jacob = new JacobianMatrices(efode, pcircle);
+        jacob.setParameterJacobianProvider(pcircle);
+        jacob.setParameterizedODE(pcircle);
+        jacob.selectParameters(Circle.CX, Circle.OMEGA);
+        jacob.setMainStateSteps(new double[] { hY, hY });
+        jacob.setParameterStep(Circle.OMEGA, hP);
+        jacob.setInitialMainStateJacobian(pcircle.exactDyDy0(0));
+        jacob.setInitialParameterJacobian(Circle.CX, pcircle.exactDyDcx(0));
+//        jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
+        jacob.setInitialParameterJacobian(Circle.OMEGA, pcircle.exactDyDom(0));
+        integ.setMaxEvaluations(50000);
+
+        double t = 18 * FastMath.PI;
+        integ.integrate(efode, 0, y, t, y);
+        for (int i = 0; i < y.length; ++i) {
+            Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9);
+        }
+
+        double[][] dydy0 = new double[2][2];
+        jacob.getCurrentMainSetJacobian(dydy0);
+        for (int i = 0; i < dydy0.length; ++i) {
+            for (int j = 0; j < dydy0[i].length; ++j) {
+                Assert.assertEquals(pcircle.exactDyDy0(t)[i][j], dydy0[i][j], 5.0e-4);
+            }
+        }
+
+        double[] dydp0 = new double[2];
+        jacob.getCurrentParameterJacobian(Circle.CX, dydp0);
+        for (int i = 0; i < dydp0.length; ++i) {
+            Assert.assertEquals(pcircle.exactDyDcx(t)[i], dydp0[i], 5.0e-4);
+        }
+
+        double[] dydp1 = new double[2];
+        jacob.getCurrentParameterJacobian(Circle.OMEGA, dydp1);
+        for (int i = 0; i < dydp1.length; ++i) {
+            Assert.assertEquals(pcircle.exactDyDom(t)[i], dydp1[i], 1.0e-2);
+        }
+    }
+
+    private static class Brusselator extends AbstractParameterizable
+        implements MainStateJacobianProvider, ParameterJacobianProvider {
+
+        public static final String B = "b";
+
+        private double b;
+
+        public Brusselator(double b) {
+            super(B);
+            this.b = b;
+        }
+
+        public int getDimension() {
+            return 2;
+        }
+
+        public void computeDerivatives(double t, double[] y, double[] yDot) {
+            double prod = y[0] * y[0] * y[1];
+            yDot[0] = 1 + prod - (b + 1) * y[0];
+            yDot[1] = b * y[0] - prod;
+        }
+
+        public void computeMainStateJacobian(double t, double[] y, double[] yDot,
+                                             double[][] dFdY) {
+            double p = 2 * y[0] * y[1];
+            double y02 = y[0] * y[0];
+            dFdY[0][0] = p - (1 + b);
+            dFdY[0][1] = y02;
+            dFdY[1][0] = b - p;
+            dFdY[1][1] = -y02;
+        }
+
+        public void computeParameterJacobian(double t, double[] y, double[] yDot,
+                                             String paramName, double[] dFdP) {
+            complainIfNotSupported(paramName);
+            dFdP[0] = -y[0];
+            dFdP[1] = y[0];
+        }
+
+        public double dYdP0() {
+            return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
+        }
+
+        public double dYdP1() {
+            return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
+        }
+
+    }
+
+    private static class ParamBrusselator extends AbstractParameterizable
+        implements FirstOrderDifferentialEquations, ParameterizedODE {
+
+        public static final String B = "b";
+
+        private double b;
+
+        public ParamBrusselator(double b) {
+            super(B);
+            this.b = b;
+        }
+
+        public int getDimension() {
+            return 2;
+        }
+
+        /** {@inheritDoc} */
+        public double getParameter(final String name)
+            throws IllegalArgumentException {
+            complainIfNotSupported(name);
+            return b;
+        }
+
+        /** {@inheritDoc} */
+        public void setParameter(final String name, final double value)
+            throws IllegalArgumentException {
+            complainIfNotSupported(name);
+            b = value;
+        }
+
+        public void computeDerivatives(double t, double[] y, double[] yDot) {
+            double prod = y[0] * y[0] * y[1];
+            yDot[0] = 1 + prod - (b + 1) * y[0];
+            yDot[1] = b * y[0] - prod;
+        }
+
+        public double dYdP0() {
+            return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
+        }
+
+        public double dYdP1() {
+            return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
+        }
+
+    }
+
+    /** ODE representing a point moving on a circle with provided center and angular rate. */
+    private static class Circle extends AbstractParameterizable
+        implements MainStateJacobianProvider, ParameterJacobianProvider {
+
+        public static final String CX = "cx";
+        public static final String CY = "cy";
+        public static final String OMEGA = "omega";
+
+        private final double[] y0;
+        private double cx;
+        private double cy;
+        private double omega;
+
+        public Circle(double[] y0, double cx, double cy, double omega) {
+            super(CX,CY,OMEGA);
+            this.y0    = y0.clone();
+            this.cx    = cx;
+            this.cy    = cy;
+            this.omega = omega;
+        }
+
+        public int getDimension() {
+            return 2;
+        }
+
+        public void computeDerivatives(double t, double[] y, double[] yDot) {
+            yDot[0] = omega * (cy - y[1]);
+            yDot[1] = omega * (y[0] - cx);
+        }
+
+        public void computeMainStateJacobian(double t, double[] y,
+                                             double[] yDot, double[][] dFdY) {
+            dFdY[0][0] = 0;
+            dFdY[0][1] = -omega;
+            dFdY[1][0] = omega;
+            dFdY[1][1] = 0;
+        }
+
+        public void computeParameterJacobian(double t, double[] y, double[] yDot,
+                                             String paramName, double[] dFdP) {
+            complainIfNotSupported(paramName);
+            if (paramName.equals(CX)) {
+                dFdP[0] = 0;
+                dFdP[1] = -omega;
+            } else if (paramName.equals(CY)) {
+                dFdP[0] = omega;
+                dFdP[1] = 0;
+            }  else {
+                dFdP[0] = cy - y[1];
+                dFdP[1] = y[0] - cx;
+            }           
+        }
+
+        public double[] exactY(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            double dx0 = y0[0] - cx;
+            double dy0 = y0[1] - cy;
+            return new double[] {
+                cx + cos * dx0 - sin * dy0,
+                cy + sin * dx0 + cos * dy0
+            };
+        }
+
+        public double[][] exactDyDy0(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            return new double[][] {
+                { cos, -sin },
+                { sin,  cos }
+            };
+        }
+
+        public double[] exactDyDcx(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            return new double[] {1 - cos, -sin};
+        }
+
+        public double[] exactDyDcy(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            return new double[] {sin, 1 - cos};
+        }
+
+        public double[] exactDyDom(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            double dx0 = y0[0] - cx;
+            double dy0 = y0[1] - cy;
+            return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
+        }
+
+        public double[][] exactDyDp(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            double dx0 = y0[0] - cx;
+            double dy0 = y0[1] - cy;
+            return new double[][] {
+                { 1 - cos, sin,    -t * (sin * dx0 + cos * dy0) },
+                { -sin,    1 - cos, t * (cos * dx0 - sin * dy0) }
+            };
+        }
+
+        public double[] exactYDot(double t) {
+            double oCos = omega * FastMath.cos(omega * t);
+            double oSin = omega * FastMath.sin(omega * t);
+            double dx0 = y0[0] - cx;
+            double dy0 = y0[1] - cy;
+            return new double[] {
+                -oSin * dx0 - oCos * dy0,
+                 oCos * dx0 - oSin * dy0
+            };
+        }
+
+        public double[][] exactDyDy0Dot(double t) {
+            double oCos = omega * FastMath.cos(omega * t);
+            double oSin = omega * FastMath.sin(omega * t);
+            return new double[][] {
+                { -oSin, -oCos },
+                {  oCos, -oSin }
+            };
+        }
+
+        public double[][] exactDyDpDot(double t) {
+            double cos  = FastMath.cos(omega * t);
+            double sin  = FastMath.sin(omega * t);
+            double oCos = omega * cos;
+            double oSin = omega * sin;
+            double dx0  = y0[0] - cx;
+            double dy0  = y0[1] - cy;
+            return new double[][] {
+                {  oSin, oCos, -sin * dx0 - cos * dy0 - t * ( oCos * dx0 - oSin * dy0) },
+                { -oCos, oSin,  cos * dx0 - sin * dy0 + t * (-oSin * dx0 - oCos * dy0) }
+            };
+        }
+
+    }
+
+    /** ODE representing a point moving on a circle with provided center and angular rate. */
+    private static class ParameterizedCircle extends AbstractParameterizable
+        implements FirstOrderDifferentialEquations, ParameterJacobianProvider, ParameterizedODE {
+
+        public static final String CX = "cx";
+        public static final String CY = "cy";
+        public static final String OMEGA = "omega";
+
+        private final double[] y0;
+        private double cx;
+        private double cy;
+        private double omega;
+
+        public ParameterizedCircle(double[] y0, double cx, double cy, double omega) {
+            super(CX,CY,OMEGA);
+            this.y0    = y0.clone();
+            this.cx    = cx;
+            this.cy    = cy;
+            this.omega = omega;
+        }
+
+        public int getDimension() {
+            return 2;
+        }
+
+        public void computeDerivatives(double t, double[] y, double[] yDot) {
+            yDot[0] = omega * (cy - y[1]);
+            yDot[1] = omega * (y[0] - cx);
+        }
+
+        public void computeParameterJacobian(double t, double[] y, double[] yDot,
+                                             String paramName, double[] dFdP)
+            throws IllegalArgumentException {
+            if (paramName.equals(CX)) {
+                dFdP[0] = 0;
+                dFdP[1] = -omega;
+            } else {
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public double getParameter(final String name)
+            throws IllegalArgumentException {
+            if (name.equals(CY)) {
+                return cy;
+            } else if (name.equals(OMEGA)) {
+                return omega;
+            } else {
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public void setParameter(final String name, final double value)
+            throws IllegalArgumentException {
+            if (name.equals(CY)) {
+                cy = value;
+            } else if (name.equals(OMEGA)) {
+                omega = value;
+            } else {
+                throw new IllegalArgumentException();
+            }
+        }
+
+        public double[] exactY(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            double dx0 = y0[0] - cx;
+            double dy0 = y0[1] - cy;
+            return new double[] {
+                cx + cos * dx0 - sin * dy0,
+                cy + sin * dx0 + cos * dy0
+            };
+        }
+
+        public double[][] exactDyDy0(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            return new double[][] {
+                { cos, -sin },
+                { sin,  cos }
+            };
+        }
+
+        public double[] exactDyDcx(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            return new double[] {1 - cos, -sin};
+        }
+
+        public double[] exactDyDcy(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            return new double[] {sin, 1 - cos};
+        }
+
+        public double[] exactDyDom(double t) {
+            double cos = FastMath.cos(omega * t);
+            double sin = FastMath.sin(omega * t);
+            double dx0 = y0[0] - cx;
+            double dy0 = y0[1] - cy;
+            return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
+        }
+
+    }
+
+}

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

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

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java?rev=1175409&r1=1175408&r2=1175409&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java Sun Sep 25 15:04:39 2011
@@ -256,8 +256,7 @@ public class GraggBulirschStoerIntegrato
   }
 
   @Test
-  public void testUnstableDerivative()
-    {
+  public void testUnstableDerivative() {
     final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
     FirstOrderIntegrator integ =
       new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);



Mime
View raw message