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

Fixed some issues in nth root derivatives at 0. 

The current behavior is correct with respect to infinities and NaN being
appropriately generated, but in some cases counter-intuitive.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1391451 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2012-09-28 13:28:42 +00:00
parent e1ed0b96ba
commit 6129654bc2
4 changed files with 75 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/changes/changes.xml
View File

@ -52,6 +52,9 @@ If the output is not quite correct, check for invisible trailing spaces!
<body>
<release version="3.1" date="TBD" description="
">
<action dev="luc" type="fix" >
Fixed some issues in nth root derivatives at 0.
</action>
<action dev="tn" type="fix" issue="MATH-848">
Fixed transformation to a Schur matrix for certain input matrices.
</action>

11
src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java
View File

@ -952,15 +952,18 @@ public class DSCompiler {
double[] function = new double[1 + order];
double xk;
if (n == 2) {
xk = FastMath.sqrt(operand[operandOffset]);
function[0] = FastMath.sqrt(operand[operandOffset]);
xk = 0.5 / function[0];
} else if (n == 3) {
xk = FastMath.cbrt(operand[operandOffset]);
function[0] = FastMath.cbrt(operand[operandOffset]);
xk = 1.0 / (3.0 * function[0] * function[0]);
} else {
xk = FastMath.pow(operand[operandOffset], 1.0 / n);
function[0] = FastMath.pow(operand[operandOffset], 1.0 / n);
xk = 1.0 / (n * FastMath.pow(function[0], n - 1));
}
final double nReciprocal = 1.0 / n;
final double xReciprocal = 1.0 / operand[operandOffset];
for (int i = 0; i <= order; ++i) {
for (int i = 1; i <= order; ++i) {
function[i] = xk;
xk *= xReciprocal * (nReciprocal - i);
}

65
src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
View File

@ -337,6 +337,69 @@ public class DerivativeStructureTest {
}
}
@Test
public void testRootNSingularity() {
for (int n = 2; n < 10; ++n) {
for (int maxOrder = 0; maxOrder < 12; ++maxOrder) {
DerivativeStructure dsZero = new DerivativeStructure(1, maxOrder, 0, 0.0);
DerivativeStructure rootN = dsZero.rootN(n);
Assert.assertEquals(0.0, rootN.getValue(), 1.0e-20);
if (maxOrder > 0) {
Assert.assertTrue(Double.isInfinite(rootN.getPartialDerivative(1)));
Assert.assertTrue(rootN.getPartialDerivative(1) > 0);
for (int order = 2; order <= maxOrder; ++order) {
// the following checks shows a LIMITATION of the current implementation
// we have no way to tell dsZero is a pure linear variable x = 0
// we only say: "dsZero is a structure with value = 0.0,
// first derivative = 1.0, second and higher derivatives = 0.0".
// Function composition rule for second derivatives is:
// d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
// when function f is the nth root and x = 0 we have:
// f(0) = 0, f'(0) = +infinity, f''(0) = -infinity (and higher
// derivatives keep switching between +infinity and -infinity)
// so given that in our case dsZero represents g, we have g(x) = 0,
// g'(x) = 1 and g''(x) = 0
// applying the composition rules gives:
// d2[f(g(x))]/dx2 = f''(g(x)) * [g'(x)]^2 + f'(g(x)) * g''(x)
// = -infinity * 1^2 + +infinity * 0
// = -infinity + NaN
// = NaN
// if we knew dsZero is really the x variable and not the identity
// function applied to x, we would not have computed f'(g(x)) * g''(x)
// and we would have found that the result was -infinity and not NaN
Assert.assertTrue(Double.isNaN(rootN.getPartialDerivative(order)));
}
}
// the following shows that the limitation explained above is NOT a bug...
// if we set up the higher order derivatives for g appropriately, we do
// compute the higher order derivatives of the composition correctly
double[] gDerivatives = new double[ 1 + maxOrder];
gDerivatives[0] = 0.0;
for (int k = 1; k <= maxOrder; ++k) {
gDerivatives[k] = FastMath.pow(-1.0, k + 1);
}
DerivativeStructure correctRoot = new DerivativeStructure(1, maxOrder, gDerivatives).rootN(n);
Assert.assertEquals(0.0, correctRoot.getValue(), 1.0e-20);
if (maxOrder > 0) {
Assert.assertTrue(Double.isInfinite(correctRoot.getPartialDerivative(1)));
Assert.assertTrue(correctRoot.getPartialDerivative(1) > 0);
for (int order = 2; order <= maxOrder; ++order) {
Assert.assertTrue(Double.isInfinite(correctRoot.getPartialDerivative(order)));
if ((order % 2) == 0) {
Assert.assertTrue(correctRoot.getPartialDerivative(order) < 0);
} else {
Assert.assertTrue(correctRoot.getPartialDerivative(order) > 0);
}
}
}
}
}
}
@Test
public void testSqrtPow2() {
double[] epsilon = new double[] { 1.0e-16, 3.0e-16, 2.0e-15, 6.0e-14, 6.0e-12 };
@ -662,7 +725,7 @@ public class DerivativeStructureTest {
@Test
public void testAtan2() {
double[] epsilon = new double[] { 5.0e-16, 3.0e-15, 2.0e-14, 1.0e-12, 8.0e-11 };
double[] epsilon = new double[] { 5.0e-16, 3.0e-15, 2.2e-14, 1.0e-12, 8.0e-11 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = -1.7; x < 2; x += 0.2) {
DerivativeStructure dsX = new DerivativeStructure(2, maxOrder, 0, x);

2
src/test/java/org/apache/commons/math3/analysis/function/SqrtTest.java
View File

@ -65,7 +65,7 @@ public class SqrtTest {
Assert.assertEquals(-0.1901814435781826783, s.getPartialDerivative(2), 1.0e-16);
Assert.assertEquals(0.23772680447272834785, s.getPartialDerivative(3), 1.0e-16);
Assert.assertEquals(-0.49526417598485072465, s.getPartialDerivative(4), 1.0e-16);
Assert.assertEquals(1.4445205132891479465, s.getPartialDerivative(5), 3.0e-16);
Assert.assertEquals(1.4445205132891479465, s.getPartialDerivative(5), 5.0e-16);
}
}