Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

fixed an index error in dy/dy0 computation 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@919166 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2010-03-04 20:40:49 +00:00
parent 5256967cf7
commit 21c37897f0
2 changed files with 63 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
View File

@ -242,27 +242,29 @@ public class FirstOrderIntegratorWithJacobians {
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 @@ public class FirstOrderIntegratorWithJacobians {
/** {@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} */

85
src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
View File

@ -20,9 +20,8 @@ package org.apache.commons.math.ode.jacobians;
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 @@ public class FirstOrderIntegratorWithJacobiansTest {
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 class FirstOrderIntegratorWithJacobiansTest {
}
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;
}
};