commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r919166 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
Date Thu, 04 Mar 2010 20:40:49 GMT
Author: luc
Date: Thu Mar  4 20:40:49 2010
New Revision: 919166

URL: http://svn.apache.org/viewvc?rev=919166&view=rev
Log:
fixed an index error in dy/dy0 computation

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java

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=919166&r1=919165&r2=919166&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
Thu Mar  4 20:40:49 2010
@@ -242,27 +242,29 @@
                     final double[] dFdYi = dFdY[i];
                     for (int j = 0; j < n; ++j) {
                         double s = 0;
-                        int zIndex = n + j;
+                        final int startIndex = n + j;
+                        int zIndex = startIndex;
                         for (int l = 0; l < n; ++l) {
                             s += dFdYi[l] * z[zIndex];
-                            zIndex += l;
+                            zIndex += n;
                         }
-                        zDot[n + i * n + j] = s;
+                        zDot[startIndex + i * n] = s;
                     }
                 }
 
-                // variational equations: d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt
+                // variational equations: from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt
                 for (int i = 0; i < n; ++i) {
                     final double[] dFdYi = dFdY[i];
                     final double[] dFdPi = dFdP[i];
                     for (int j = 0; j < k; ++j) {
                         double s = dFdPi[j];
-                        int zIndex = n * (n + 1)+ j;
+                        final int startIndex = n * (n + 1) + j;
+                        int zIndex = startIndex;
                         for (int l = 0; l < n; ++l) {
                             s += dFdYi[l] * z[zIndex];
                             zIndex += k;
                         }
-                        zDot[n * (n + 1) + i * k + j] = s;
+                        zDot[startIndex + i * k] = s;
                     }
                 }
 
@@ -549,8 +551,19 @@
 
         /** {@inheritDoc} */
         public StepInterpolatorWithJacobians copy() throws DerivativeException {
-            return new StepInterpolatorWrapper(interpolator.copy(),
-                                               y.length, dydy0[0].length);
+            final int n = y.length;
+            final int k = dydp[0].length;
+            StepInterpolatorWrapper copied =
+                new StepInterpolatorWrapper(interpolator.copy(), n, k);
+            System.arraycopy(y,    0, copied.y,    0, n);
+            System.arraycopy(yDot, 0, copied.yDot, 0, n);
+            for (int i = 0; i < n; ++i) {
+                System.arraycopy(dydy0[i], 0, copied.dydy0[i], 0, n);
+            }
+            for (int i = 0; i < n; ++i) {
+                System.arraycopy(dydp[i], 0, copied.dydp[i], 0, k);
+            }
+            return copied;
         }
 
         /** {@inheritDoc} */

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=919166&r1=919165&r2=919166&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
Thu Mar  4 20:40:49 2010
@@ -20,9 +20,8 @@
 import org.apache.commons.math.ode.DerivativeException;
 import org.apache.commons.math.ode.FirstOrderIntegrator;
 import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.jacobians.FirstOrderIntegratorWithJacobians;
-import org.apache.commons.math.ode.jacobians.ParameterizedODEWithJacobians;
 import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator;
+import org.apache.commons.math.optimization.OptimizationException;
 import org.apache.commons.math.stat.descriptive.SummaryStatistics;
 import org.junit.Assert;
 import org.junit.Test;
@@ -35,8 +34,8 @@
         FirstOrderIntegrator integ =
             new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4);
         double hP = 1.0e-12;
-        SummaryStatistics residuals0 = new SummaryStatistics();
-        SummaryStatistics residuals1 = new SummaryStatistics();
+        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 };
@@ -44,13 +43,13 @@
             double[] yP = { 1.3, b + hP };
             brusselator.setParameter(0, b + hP);
             integ.integrate(brusselator, 0, yP, 20.0, yP);
-            residuals0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
-            residuals1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
+            residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
+            residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
         }
-        Assert.assertTrue((residuals0.getMax() - residuals0.getMin()) > 600);
-        Assert.assertTrue(residuals0.getStandardDeviation() > 30);
-        Assert.assertTrue((residuals1.getMax() - residuals1.getMin()) > 800);
-        Assert.assertTrue(residuals1.getStandardDeviation() > 50);
+        Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 600);
+        Assert.assertTrue(residualsP0.getStandardDeviation() > 30);
+        Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 800);
+        Assert.assertTrue(residualsP1.getStandardDeviation() > 50);
     }
 
     @Test
@@ -59,8 +58,8 @@
         FirstOrderIntegrator integ =
             new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
         double hP = 1.0e-12;
-        SummaryStatistics residuals0 = new SummaryStatistics();
-        SummaryStatistics residuals1 = new SummaryStatistics();
+        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 };
@@ -68,17 +67,17 @@
             double[] yP = { 1.3, b + hP };
             brusselator.setParameter(0, b + hP);
             integ.integrate(brusselator, 0, yP, 20.0, yP);
-            residuals0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
-            residuals1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
+            residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
+            residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
         }
-        Assert.assertTrue((residuals0.getMax() - residuals0.getMin()) > 0.02);
-        Assert.assertTrue((residuals0.getMax() - residuals0.getMin()) < 0.03);
-        Assert.assertTrue(residuals0.getStandardDeviation() > 0.003);
-        Assert.assertTrue(residuals0.getStandardDeviation() < 0.004);
-        Assert.assertTrue((residuals1.getMax() - residuals1.getMin()) > 0.04);
-        Assert.assertTrue((residuals1.getMax() - residuals1.getMin()) < 0.05);
-        Assert.assertTrue(residuals1.getStandardDeviation() > 0.006);
-        Assert.assertTrue(residuals1.getStandardDeviation() < 0.007);
+        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.006);
+        Assert.assertTrue(residualsP1.getStandardDeviation() < 0.007);
     }
 
     @Test
@@ -87,8 +86,8 @@
         FirstOrderIntegrator integ =
             new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4);
         double hP = 1.0e-12;
-        SummaryStatistics residuals0 = new SummaryStatistics();
-        SummaryStatistics residuals1 = new SummaryStatistics();
+        SummaryStatistics residualsP0 = new SummaryStatistics();
+        SummaryStatistics residualsP1 = new SummaryStatistics();
         for (double b = 2.88; b < 3.08; b += 0.001) {
             Brusselator brusselator = new Brusselator(b);
             brusselator.setParameter(0, b);
@@ -100,22 +99,22 @@
                 new FirstOrderIntegratorWithJacobians(integ, brusselator, new double[] {
b },
                                                  new double[] { hY, hY }, new double[] {
hP });
             extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
-            residuals0.addValue(dZdP[0][0] - brusselator.dYdP0());
-            residuals1.addValue(dZdP[1][0] - brusselator.dYdP1());
+            residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
+            residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
         }
-        Assert.assertTrue((residuals0.getMax() - residuals0.getMin()) < 0.006);
-        Assert.assertTrue(residuals0.getStandardDeviation() < 0.0009);
-        Assert.assertTrue((residuals1.getMax() - residuals1.getMin()) < 0.006);
-        Assert.assertTrue(residuals1.getStandardDeviation() < 0.0012);
+        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);
     }
 
     @Test
     public void testAnalyticalDifferentiation()
-        throws IntegratorException, DerivativeException {
+        throws IntegratorException, DerivativeException, OptimizationException {
         FirstOrderIntegrator integ =
-            new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4);
-        SummaryStatistics residuals0 = new SummaryStatistics();
-        SummaryStatistics residuals1 = new SummaryStatistics();
+            new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
+        SummaryStatistics residualsP0 = new SummaryStatistics();
+        SummaryStatistics residualsP1 = new SummaryStatistics();
         for (double b = 2.88; b < 3.08; b += 0.001) {
             Brusselator brusselator = new Brusselator(b);
             brusselator.setParameter(0, b);
@@ -125,13 +124,13 @@
             FirstOrderIntegratorWithJacobians extInt =
                 new FirstOrderIntegratorWithJacobians(integ, brusselator);
             extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
-            residuals0.addValue(dZdP[0][0] - brusselator.dYdP0());
-            residuals1.addValue(dZdP[1][0] - brusselator.dYdP1());
-       }
-        Assert.assertTrue((residuals0.getMax() - residuals0.getMin()) < 0.004);
-        Assert.assertTrue(residuals0.getStandardDeviation() < 0.001);
-        Assert.assertTrue((residuals1.getMax() - residuals1.getMin()) < 0.005);
-        Assert.assertTrue(residuals1.getStandardDeviation() < 0.001);
+            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);
     }
 
     private static class Brusselator implements ParameterizedODEWithJacobians {
@@ -172,11 +171,11 @@
         }
 
         public double dYdP0() {
-            return -1087.8787631970476 + (1050.4387741821572 + (-338.90621620263096 + 36.51793006801084
* b) * b) * b;
+            return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327
* b) * b) * b;
         }
 
         public double dYdP1() {
-            return 1499.0904666097015 + (-1434.9574631810726 + (459.71079478756945 - 49.29949940968588
* b) * b) * b;
+            return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082
* b) * b) * b;
         }
 
     };



Mime
View raw message