Author: luc
Date: Fri Aug 10 18:33:06 2012
New Revision: 1371805
URL: http://svn.apache.org/viewvc?rev=1371805&view=rev
Log:
Completed support fo asin, acos and atan in DSCompiler.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java?rev=1371805&r1=1371804&r2=1371805&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java
Fri Aug 10 18:33:06 2012
@@ -1097,10 +1097,39 @@ public class DSCompiler {
final double x = operand[operandOffset];
function[0] = FastMath.acos(x);
if (order > 0) {
- function[1] = -1.0 / FastMath.sqrt(1 - x * x);
- for (int i = 2; i <= order; ++i) {
- // TODO compute higher order derivatives
- function[i] = Double.NaN;
+ // the nth order derivative of acos has the form:
+ // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
+ // where P_n(x) is a degree n-1 polynomial with same parity as n-1
+ // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ...
+ // the general recurrence relation for P_n is:
+ // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
+ // as per polynomial parity, we can store coefficients of both
P_(n-1) and P_n in the same array
+ final double[] p = new double[order];
+ p[0] = -1;
+ final double x2 = x * x;
+ final double f = 1.0 / (1 - x2);
+ double coeff = FastMath.sqrt(f);
+ function[1] = coeff * p[0];
+ for (int n = 2; n <= order; ++n) {
+
+ // update and evaluate polynomial P_n(x)
+ double v = 0;
+ p[n - 1] = (n - 1) * p[n - 2];
+ for (int k = n - 1; k >= 0; k -= 2) {
+ v = v * x2 + p[k];
+ if (k > 2) {
+ p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
+ } else if (k == 2) {
+ p[0] = p[1];
+ }
+ }
+ if ((n & 0x1) == 0) {
+ v *= x;
+ }
+
+ coeff *= f;
+ function[n] = coeff * v;
+
}
}
@@ -1125,10 +1154,39 @@ public class DSCompiler {
final double x = operand[operandOffset];
function[0] = FastMath.asin(x);
if (order > 0) {
- function[1] = 1.0 / FastMath.sqrt(1 - x * x);
- for (int i = 2; i <= order; ++i) {
- // TODO compute higher order derivatives
- function[i] = Double.NaN;
+ // the nth order derivative of asin has the form:
+ // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
+ // where P_n(x) is a degree n-1 polynomial with same parity as n-1
+ // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
+ // the general recurrence relation for P_n is:
+ // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
+ // as per polynomial parity, we can store coefficients of both
P_(n-1) and P_n in the same array
+ final double[] p = new double[order];
+ p[0] = 1;
+ final double x2 = x * x;
+ final double f = 1.0 / (1 - x2);
+ double coeff = FastMath.sqrt(f);
+ function[1] = coeff * p[0];
+ for (int n = 2; n <= order; ++n) {
+
+ // update and evaluate polynomial P_n(x)
+ double v = 0;
+ p[n - 1] = (n - 1) * p[n - 2];
+ for (int k = n - 1; k >= 0; k -= 2) {
+ v = v * x2 + p[k];
+ if (k > 2) {
+ p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
+ } else if (k == 2) {
+ p[0] = p[1];
+ }
+ }
+ if ((n & 0x1) == 0) {
+ v *= x;
+ }
+
+ coeff *= f;
+ function[n] = coeff * v;
+
}
}
@@ -1153,10 +1211,39 @@ public class DSCompiler {
final double x = operand[operandOffset];
function[0] = FastMath.atan(x);
if (order > 0) {
- function[1] = 1.0 / (1 + x * x);
- for (int i = 2; i <= order; ++i) {
- // TODO compute higher order derivatives
- function[i] = Double.NaN;
+ // the nth order derivative of atan has the form:
+ // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n
+ // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
+ // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ...
+ // the general recurrence relation for Q_n is:
+ // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x)
+ // as per polynomial parity, we can store coefficients of both
Q_(n-1) and Q_n in the same array
+ final double[] q = new double[order];
+ q[0] = 1;
+ final double x2 = x * x;
+ final double f = 1.0 / (1 + x2);
+ double coeff = f;
+ function[1] = coeff * q[0];
+ for (int n = 2; n <= order; ++n) {
+
+ // update and evaluate polynomial Q_n(x)
+ double v = 0;
+ q[n - 1] = -n * q[n - 2];
+ for (int k = n - 1; k >= 0; k -= 2) {
+ v = v * x2 + q[k];
+ if (k > 2) {
+ q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k
- 3];
+ } else if (k == 2) {
+ q[0] = q[1];
+ }
+ }
+ if ((n & 0x1) == 0) {
+ v *= x;
+ }
+
+ coeff *= f;
+ function[n] = coeff * v;
+
}
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java?rev=1371805&r1=1371804&r2=1371805&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
Fri Aug 10 18:33:06 2012
@@ -460,6 +460,51 @@ public class DerivativeStructureTest {
}
@Test
+ public void testSinAsin() {
+ double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 3.0e-15, 2.0e-14,
4.0e-13 };
+ for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+ for (double x = 0.1; x < 1.2; x += 0.001) {
+ DerivativeStructure dsX = new DerivativeStructure(1, maxOrder,
0, x);
+ DerivativeStructure rebuiltX = dsX.sin().asin();
+ DerivativeStructure zero = rebuiltX.subtract(dsX);
+ for (int n = 0; n <= maxOrder; ++n) {
+ Assert.assertEquals(0.0, zero.getPartialDerivative(n),
epsilon[n]);
+ }
+ }
+ }
+ }
+
+ @Test
+ public void testCosAcos() {
+ double[] epsilon = new double[] { 6.0e-16, 6.0e-15, 2.0e-13, 4.0e-12,
2.0e-10 };
+ for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+ for (double x = 0.1; x < 1.2; x += 0.001) {
+ DerivativeStructure dsX = new DerivativeStructure(1, maxOrder,
0, x);
+ DerivativeStructure rebuiltX = dsX.cos().acos();
+ DerivativeStructure zero = rebuiltX.subtract(dsX);
+ for (int n = 0; n <= maxOrder; ++n) {
+ Assert.assertEquals(0.0, zero.getPartialDerivative(n),
epsilon[n]);
+ }
+ }
+ }
+ }
+
+ @Test
+ public void testTanAtan() {
+ double[] epsilon = new double[] { 6.0e-17, 2.0e-16, 2.0e-15, 4.0e-14,
2.0e-12 };
+ for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+ for (double x = 0.1; x < 1.2; x += 0.001) {
+ DerivativeStructure dsX = new DerivativeStructure(1, maxOrder,
0, x);
+ DerivativeStructure rebuiltX = dsX.tan().atan();
+ DerivativeStructure zero = rebuiltX.subtract(dsX);
+ for (int n = 0; n <= maxOrder; ++n) {
+ Assert.assertEquals(0.0, zero.getPartialDerivative(n),
epsilon[n]);
+ }
+ }
+ }
+ }
+
+ @Test
public void testTangentDefinition() {
double[] epsilon = new double[] { 5.0e-16, 2.0e-15, 3.0e-14, 5.0e-13,
2.0e-11 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
