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

Added hyperbolic trigonometric functions and inverses to DSCompiler. 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1372902 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2012-08-14 14:25:50 +00:00
parent 6f859c5f5b
commit ff43dac3d2
3 changed files with 426 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1328,6 +1328,284 @@ public class DSCompiler {
}
/** Compute hyperbolic cosine of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param result array where result must be stored (for
* hyperbolic cosine the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void cosh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.cosh(operand[operandOffset]);
if (order > 0) {
function[1] = FastMath.sinh(operand[operandOffset]);
for (int i = 2; i <= order; ++i) {
function[i] = function[i - 2];
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute hyperbolic sine of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param result array where result must be stored (for
* hyperbolic sine the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void sinh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.sinh(operand[operandOffset]);
if (order > 0) {
function[1] = FastMath.cosh(operand[operandOffset]);
for (int i = 2; i <= order; ++i) {
function[i] = function[i - 2];
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute hyperbolic tangent of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param result array where result must be stored (for
* hyperbolic tangent the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void tanh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
final double[] function = new double[1 + order];
final double t = FastMath.tanh(operand[operandOffset]);
function[0] = t;
if (order > 0) {
// the nth order derivative of tanh has the form:
// dn(tanh(x)/dxn = P_n(tanh(x))
// where P_n(t) is a degree n+1 polynomial with same parity as n+1
// P_0(t) = t, P_1(t) = 1 - t^2, P_2(x) = -2 t (1 - t^2) ...
// the general recurrence relation for P_n is:
// P_n(x) = (1-t^2) P_(n-1)'(t)
// 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 + 2];
p[1] = 1;
final double t2 = t * t;
for (int n = 1; n <= order; ++n) {
// update and evaluate polynomial P_n(t)
double v = 0;
p[n + 1] = -n * p[n];
for (int k = n + 1; k >= 0; k -= 2) {
v = v * t2 + p[k];
if (k > 2) {
p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3];
} else if (k == 2) {
p[0] = p[1];
}
}
if ((n & 0x1) == 0) {
v *= t;
}
function[n] = v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute inverse hyperbolic cosine of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param result array where result must be stored (for
* inverse hyperbolic cosine the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void acosh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.acosh(x);
if (order > 0) {
// the nth order derivative of acosh has the form:
// dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((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) = (x^2-1) 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 / (x2 - 1);
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] = (1 - n) * p[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
v = v * x2 + p[k];
if (k > 2) {
p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3];
} else if (k == 2) {
p[0] = -p[1];
}
}
if ((n & 0x1) == 0) {
v *= x;
}
coeff *= f;
function[n] = coeff * v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute inverse hyperbolic sine of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param result array where result must be stored (for
* inverse hyperbolic sine the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void asinh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.asinh(x);
if (order > 0) {
// the nth order derivative of asinh has the form:
// dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((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) = (x^2+1) 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] = (1 - n) * 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] + (k - 2 * n) * p[k - 3];
} else if (k == 2) {
p[0] = p[1];
}
}
if ((n & 0x1) == 0) {
v *= x;
}
coeff *= f;
function[n] = coeff * v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute inverse hyperbolic tangent of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param result array where result must be stored (for
* inverse hyperbolic tangent the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void atanh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.atanh(x);
if (order > 0) {
// the nth order derivative of atanh has the form:
// dn(atanh(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] + (2 * n - k + 1) * q[k - 3];
} else if (k == 2) {
q[0] = q[1];
}
}
if ((n & 0x1) == 0) {
v *= x;
}
coeff *= f;
function[n] = coeff * v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute composition of a derivative structure by a function.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array

56
src/main/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructure.java
View File

@ -505,7 +505,7 @@ public class DerivativeStructure implements FieldElement<DerivativeStructure>, S
}
/** Arc tangent operation.
* @return tan(this)
* @return atan(this)
*/
public DerivativeStructure atan() {
final DerivativeStructure result = new DerivativeStructure(compiler);
@ -527,6 +527,60 @@ public class DerivativeStructure implements FieldElement<DerivativeStructure>, S
return result;
}
/** Hyperbolic cosine operation.
* @return cosh(this)
*/
public DerivativeStructure cosh() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.cosh(data, 0, result.data, 0);
return result;
}
/** Hyperbolic sine operation.
* @return sinh(this)
*/
public DerivativeStructure sinh() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.sinh(data, 0, result.data, 0);
return result;
}
/** Hyperbolic tangent operation.
* @return tanh(this)
*/
public DerivativeStructure tanh() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.tanh(data, 0, result.data, 0);
return result;
}
/** Inverse hyperbolic cosine operation.
* @return acosh(this)
*/
public DerivativeStructure acosh() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.acosh(data, 0, result.data, 0);
return result;
}
/** Inverse hyperbolic sine operation.
* @return asin(this)
*/
public DerivativeStructure asinh() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.asinh(data, 0, result.data, 0);
return result;
}
/** Inverse hyperbolic tangent operation.
* @return atanh(this)
*/
public DerivativeStructure atanh() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.atanh(data, 0, result.data, 0);
return result;
}
/** Evaluate Taylor expansion a derivative structure.
* @param delta parameters offsets (&Delta;x, &Delta;y, ...)
* @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...

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

@ -568,6 +568,99 @@ public class DerivativeStructureTest {
}
}
@Test
public void testSinhDefinition() {
double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 2.0e-15, 6.0e-15 };
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 sinh1 = dsX.exp().subtract(dsX.exp().reciprocal()).multiply(0.5);
DerivativeStructure sinh2 = dsX.sinh();
DerivativeStructure zero = sinh1.subtract(sinh2);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testCoshDefinition() {
double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 5.0e-16, 2.0e-15, 6.0e-15 };
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 cosh1 = dsX.exp().add(dsX.exp().reciprocal()).multiply(0.5);
DerivativeStructure cosh2 = dsX.cosh();
DerivativeStructure zero = cosh1.subtract(cosh2);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testTanhDefinition() {
double[] epsilon = new double[] { 3.0e-16, 5.0e-16, 7.0e-16, 3.0e-15, 2.0e-14 };
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 tanh1 = dsX.exp().subtract(dsX.exp().reciprocal()).divide(dsX.exp().add(dsX.exp().reciprocal()));
DerivativeStructure tanh2 = dsX.tanh();
DerivativeStructure zero = tanh1.subtract(tanh2);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testSinhAsinh() {
double[] epsilon = new double[] { 3.0e-16, 3.0e-16, 4.0e-16, 7.0e-16, 3.0e-15, 8.0e-15 };
for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure rebuiltX = dsX.sinh().asinh();
DerivativeStructure zero = rebuiltX.subtract(dsX);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testCoshAcosh() {
double[] epsilon = new double[] { 2.0e-15, 1.0e-14, 2.0e-13, 6.0e-12, 3.0e-10, 2.0e-8 };
for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure rebuiltX = dsX.cosh().acosh();
DerivativeStructure zero = rebuiltX.subtract(dsX);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testTanhAtanh() {
double[] epsilon = new double[] { 3.0e-16, 2.0e-16, 7.0e-16, 4.0e-15, 3.0e-14, 4.0e-13 };
for (int maxOrder = 0; maxOrder < 6; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure rebuiltX = dsX.tanh().atanh();
DerivativeStructure zero = rebuiltX.subtract(dsX);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testCompositionOneVariableY() {
double epsilon = 1.0e-13;