Added lower and upper boundaries to finite differences.
When a function is defined only on an interval, finite differences should not attempt to compute sample points outside this interval. This case is now detected properly and the sample is shifted if needed. This may result in some loss of accuracy as the formula is not centered anymore, but at least we should not get catastrophic errors or exceptions. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1416249 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
77da5588a6
commit
8f98daf191
|
|
@ -25,6 +25,7 @@ import org.apache.commons.math3.exception.MathIllegalArgumentException;
|
|||
import org.apache.commons.math3.exception.NotPositiveException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
|
||||
/** Univariate functions differentiator using finite differences.
|
||||
* <p>
|
||||
|
|
@ -58,8 +59,8 @@ import org.apache.commons.math3.exception.NumberIsTooSmallException;
|
|||
* <ul>
|
||||
* <li>step size = 0.25, second order derivative error about 9.97e-10</li>
|
||||
* <li>step size = 0.25, fourth order derivative error about 5.43e-8</li>
|
||||
* <li>step size = 1.0e-6, second order derivative error about 56.25</li>
|
||||
* <li>step size = 1.0e-6, fourth order derivative error about 2.47e+14</li>
|
||||
* <li>step size = 1.0e-6, second order derivative error about 148</li>
|
||||
* <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li>
|
||||
* </ul>
|
||||
* This example shows that the small step size is really bad, even simply
|
||||
* for second order derivative!
|
||||
|
|
@ -78,10 +79,19 @@ public class FiniteDifferencesDifferentiator
|
|||
private final int nbPoints;
|
||||
|
||||
/** Step size. */
|
||||
private double stepSize;
|
||||
private final double stepSize;
|
||||
|
||||
/** Half sample span. */
|
||||
private final double halfSampleSpan;
|
||||
|
||||
/** Lower bound for independent variable. */
|
||||
private final double tMin;
|
||||
|
||||
/** Upper bound for independent variable. */
|
||||
private final double tMax;
|
||||
|
||||
/**
|
||||
* Build a differentiator with number of points and step size.
|
||||
* Build a differentiator with number of points and step size when independent variable is unbounded.
|
||||
* <p>
|
||||
* Beware that wrong settings for the finite differences differentiator
|
||||
* can lead to highly unstable and inaccurate results, especially for
|
||||
|
|
@ -96,6 +106,41 @@ public class FiniteDifferencesDifferentiator
|
|||
*/
|
||||
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize)
|
||||
throws NotPositiveException, NumberIsTooSmallException {
|
||||
this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
|
||||
}
|
||||
|
||||
/**
|
||||
* Build a differentiator with number of points and step size when independent variable is bounded.
|
||||
* <p>
|
||||
* When the independent variable is bounded (tLower < t < tUpper), the sampling
|
||||
* points used for differentiation will be adapted to ensure the constraint holds
|
||||
* even near the boundaries. This means the sample will not be centered anymore in
|
||||
* these cases. At an extreme case, computing derivatives exactly at the lower bound
|
||||
* will lead the sample to be entirely on the right side of the derivation point.
|
||||
* </p>
|
||||
* <p>
|
||||
* Note that the boundaries are considered to be excluded for function evaluation.
|
||||
* </p>
|
||||
* <p>
|
||||
* Beware that wrong settings for the finite differences differentiator
|
||||
* can lead to highly unstable and inaccurate results, especially for
|
||||
* high derivation orders. Using very small step sizes is often a
|
||||
* <em>bad</em> idea.
|
||||
* </p>
|
||||
* @param nbPoints number of points to use
|
||||
* @param stepSize step size (gap between each point)
|
||||
* @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
|
||||
* if there are no lower bounds)
|
||||
* @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
|
||||
* if there are no upper bounds)
|
||||
* @exception NotPositiveException if {@code stepsize <= 0} (note that
|
||||
* {@link NotPositiveException} extends {@link NumberIsTooSmallException})
|
||||
* @exception NumberIsTooSmallException {@code nbPoint <= 1}
|
||||
* @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
|
||||
*/
|
||||
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
|
||||
final double tLower, final double tUpper)
|
||||
throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
|
||||
|
||||
if (nbPoints <= 1) {
|
||||
throw new NumberIsTooSmallException(stepSize, 1, false);
|
||||
|
|
@ -107,6 +152,14 @@ public class FiniteDifferencesDifferentiator
|
|||
}
|
||||
this.stepSize = stepSize;
|
||||
|
||||
halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
|
||||
if (2 * halfSampleSpan >= tUpper - tLower) {
|
||||
throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
|
||||
}
|
||||
final double safety = FastMath.ulp(halfSampleSpan);
|
||||
this.tMin = tLower + halfSampleSpan + safety;
|
||||
this.tMax = tUpper - halfSampleSpan - safety;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -126,14 +179,19 @@ public class FiniteDifferencesDifferentiator
|
|||
}
|
||||
|
||||
/**
|
||||
* Evaluate derivatives from a centered sample.
|
||||
* @param t central value and derivatives
|
||||
* @param y function values at {@code t + stepSize * (i - 0.5 * (nbPoints - 1))}
|
||||
* Evaluate derivatives from a sample.
|
||||
* <p>
|
||||
* Evaluation is done using divided differences.
|
||||
* </p>
|
||||
* @param te evaluation abscissa value and derivatives
|
||||
* @param t0 first sample point abscissa
|
||||
* @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)}
|
||||
* @return value and derivatives at {@code t}
|
||||
* @exception NumberIsTooLargeException if the requested derivation order
|
||||
* is larger or equal to the number of points
|
||||
*/
|
||||
private DerivativeStructure evaluate(final DerivativeStructure t, final double[] y)
|
||||
private DerivativeStructure evaluate(final DerivativeStructure t, final double t0,
|
||||
final double[] y)
|
||||
throws NumberIsTooLargeException {
|
||||
|
||||
// create divided differences diagonal arrays
|
||||
|
|
@ -154,16 +212,23 @@ public class FiniteDifferencesDifferentiator
|
|||
}
|
||||
|
||||
// evaluate interpolation polynomial (represented by top diagonal) at t
|
||||
final int order = t.getOrder();
|
||||
final int parameters = t.getFreeParameters();
|
||||
final int order = t.getOrder();
|
||||
final int parameters = t.getFreeParameters();
|
||||
final double[] derivatives = t.getAllDerivatives();
|
||||
final double dt0 = t.getValue() - t0;
|
||||
DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0);
|
||||
DerivativeStructure monomial = new DerivativeStructure(parameters, order, 1.0);
|
||||
DerivativeStructure monomial = null;
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
if (i == 0) {
|
||||
// start with monomial(t) = 1
|
||||
monomial = new DerivativeStructure(parameters, order, 1.0);
|
||||
} else {
|
||||
// monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1))
|
||||
derivatives[0] = dt0 - (i - 1) * stepSize;
|
||||
final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives);
|
||||
monomial = monomial.multiply(deltaX);
|
||||
}
|
||||
interpolation = interpolation.add(monomial.multiply(top[i]));
|
||||
derivatives[0] = stepSize * (0.5 * (nbPoints - 1) - i);
|
||||
final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives);
|
||||
monomial = monomial.multiply(deltaX);
|
||||
}
|
||||
|
||||
return interpolation;
|
||||
|
|
@ -193,16 +258,17 @@ public class FiniteDifferencesDifferentiator
|
|||
throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
|
||||
}
|
||||
|
||||
// get sample points centered around t value
|
||||
final double t0 = t.getValue();
|
||||
// compute sample position, trying to be centered if possible
|
||||
final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
|
||||
|
||||
// compute sample points
|
||||
final double[] y = new double[nbPoints];
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
final double xi = t0 + stepSize * (i - 0.5 * (nbPoints - 1));
|
||||
y[i] = function.value(xi);
|
||||
y[i] = function.value(t0 + i * stepSize);
|
||||
}
|
||||
|
||||
// evaluate derivatives
|
||||
return evaluate(t, y);
|
||||
return evaluate(t, t0, y);
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -232,12 +298,13 @@ public class FiniteDifferencesDifferentiator
|
|||
throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
|
||||
}
|
||||
|
||||
// get sample points centered around t value
|
||||
final double t0 = t.getValue();
|
||||
// compute sample position, trying to be centered if possible
|
||||
final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
|
||||
|
||||
// compute sample points
|
||||
double[][] y = null;
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
final double xi = t0 + stepSize * (i - 0.5 * (nbPoints - 1));
|
||||
final double[] v = function.value(xi);
|
||||
final double[] v = function.value(t0 + i * stepSize);
|
||||
if (i == 0) {
|
||||
y = new double[v.length][nbPoints];
|
||||
}
|
||||
|
|
@ -249,7 +316,7 @@ public class FiniteDifferencesDifferentiator
|
|||
// evaluate derivatives
|
||||
final DerivativeStructure[] value = new DerivativeStructure[y.length];
|
||||
for (int j = 0; j < value.length; ++j) {
|
||||
value[j] = evaluate(t, y[j]);
|
||||
value[j] = evaluate(t, t0, y[j]);
|
||||
}
|
||||
|
||||
return value;
|
||||
|
|
@ -282,12 +349,13 @@ public class FiniteDifferencesDifferentiator
|
|||
throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
|
||||
}
|
||||
|
||||
// get sample points centered around t value
|
||||
final double t0 = t.getValue();
|
||||
// compute sample position, trying to be centered if possible
|
||||
final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
|
||||
|
||||
// compute sample points
|
||||
double[][][] y = null;
|
||||
for (int i = 0; i < nbPoints; ++i) {
|
||||
final double xi = t0 + stepSize * (i - 0.5 * (nbPoints - 1));
|
||||
final double[][] v = function.value(xi);
|
||||
final double[][] v = function.value(t0 + i * stepSize);
|
||||
if (i == 0) {
|
||||
y = new double[v.length][v[0].length][nbPoints];
|
||||
}
|
||||
|
|
@ -302,7 +370,7 @@ public class FiniteDifferencesDifferentiator
|
|||
final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length];
|
||||
for (int j = 0; j < value.length; ++j) {
|
||||
for (int k = 0; k < y[j].length; ++k) {
|
||||
value[j][k] = evaluate(t, y[j][k]);
|
||||
value[j][k] = evaluate(t, t0, y[j][k]);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -87,9 +87,9 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
});
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
DerivativeStructure y = f.value(new DerivativeStructure(1, 2, 0, x));
|
||||
Assert.assertEquals(2 - 3 * x, y.getValue(), 1.0e-20);
|
||||
Assert.assertEquals("" + (2 - 3 * x - y.getValue()), 2 - 3 * x, y.getValue(), 2.0e-15);
|
||||
Assert.assertEquals(-3.0, y.getPartialDerivative(1), 4.0e-13);
|
||||
Assert.assertEquals( 0.0, y.getPartialDerivative(2), 5.0e-11);
|
||||
Assert.assertEquals( 0.0, y.getPartialDerivative(2), 9.0e-11);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -101,7 +101,7 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
UnivariateDifferentiableFunction f =
|
||||
differentiator.differentiate(gaussian);
|
||||
double[] expectedError = new double[] {
|
||||
2.776e-17, 1.742e-15, 2.385e-13, 1.329e-11, 2.668e-9, 8.873e-8
|
||||
6.939e-18, 1.284e-15, 2.477e-13, 1.168e-11, 2.840e-9, 7.971e-8
|
||||
};
|
||||
double[] maxError = new double[expectedError.length];
|
||||
for (double x = -10; x < 10; x += 0.1) {
|
||||
|
|
@ -152,7 +152,7 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
// the 1.0e-6 step size is far too small for finite differences in the quintic on this abscissa range for 7 points
|
||||
// the errors are huge!
|
||||
final double[] expectedBad = new double[] {
|
||||
1.792e-22, 6.926e-5, 56.25, 1.783e8, 2.468e14, 3.056e20, 5.857e26
|
||||
2.910e-11, 2.087e-5, 147.7, 3.820e7, 6.354e14, 6.548e19, 1.543e27
|
||||
};
|
||||
|
||||
for (int i = 0; i < maxErrorGood.length; ++i) {
|
||||
|
|
@ -201,6 +201,87 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
f.value(new DerivativeStructure(1, 3, 0, 1.0));
|
||||
}
|
||||
|
||||
@Test(expected=NumberIsTooLargeException.class)
|
||||
public void testTooLargeStep() {
|
||||
new FiniteDifferencesDifferentiator(3, 2.5, 0.0, 1.0);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testBounds() {
|
||||
|
||||
final double slope = 2.5;
|
||||
UnivariateFunction f = new UnivariateFunction() {
|
||||
public double value(double x) {
|
||||
if (x < 0) {
|
||||
throw new NumberIsTooSmallException(x, 0, true);
|
||||
} else if (x > 1) {
|
||||
throw new NumberIsTooLargeException(x, 1, true);
|
||||
} else {
|
||||
return slope * x;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
UnivariateDifferentiableFunction missingBounds =
|
||||
new FiniteDifferencesDifferentiator(3, 0.1).differentiate(f);
|
||||
UnivariateDifferentiableFunction properlyBounded =
|
||||
new FiniteDifferencesDifferentiator(3, 0.1, 0.0, 1.0).differentiate(f);
|
||||
DerivativeStructure tLow = new DerivativeStructure(1, 1, 0, 0.05);
|
||||
DerivativeStructure tHigh = new DerivativeStructure(1, 1, 0, 0.95);
|
||||
|
||||
try {
|
||||
// here, we did not set the bounds, so the differences are evaluated out of domain
|
||||
// using f(-0.05), f(0.05), f(0.15)
|
||||
missingBounds.value(tLow);
|
||||
Assert.fail("an exception should have been thrown");
|
||||
} catch (NumberIsTooSmallException nse) {
|
||||
Assert.assertEquals(-0.05, nse.getArgument().doubleValue(), 1.0e-10);
|
||||
} catch (Exception e) {
|
||||
Assert.fail("wrong exception caught: " + e.getClass().getName());
|
||||
}
|
||||
|
||||
try {
|
||||
// here, we did not set the bounds, so the differences are evaluated out of domain
|
||||
// using f(0.85), f(0.95), f(1.05)
|
||||
missingBounds.value(tHigh);
|
||||
Assert.fail("an exception should have been thrown");
|
||||
} catch (NumberIsTooLargeException nle) {
|
||||
Assert.assertEquals(1.05, nle.getArgument().doubleValue(), 1.0e-10);
|
||||
} catch (Exception e) {
|
||||
Assert.fail("wrong exception caught: " + e.getClass().getName());
|
||||
}
|
||||
|
||||
// here, we did set the bounds, so evaluations are done within domain
|
||||
// using f(0.0), f(0.1), f(0.2)
|
||||
Assert.assertEquals(slope, properlyBounded.value(tLow).getPartialDerivative(1), 1.0e-10);
|
||||
|
||||
// here, we did set the bounds, so evaluations are done within domain
|
||||
// using f(0.8), f(0.9), f(1.0)
|
||||
Assert.assertEquals(slope, properlyBounded.value(tHigh).getPartialDerivative(1), 1.0e-10);
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testBoundedSqrt() {
|
||||
|
||||
UnivariateFunctionDifferentiator differentiator =
|
||||
new FiniteDifferencesDifferentiator(9, 1.0 / 32, 0.0, Double.POSITIVE_INFINITY);
|
||||
UnivariateDifferentiableFunction sqrt = differentiator.differentiate(new UnivariateFunction() {
|
||||
public double value(double x) {
|
||||
return FastMath.sqrt(x);
|
||||
}
|
||||
});
|
||||
|
||||
// we are able to compute derivative near 0, but the accuracy is much poorer there
|
||||
DerivativeStructure t001 = new DerivativeStructure(1, 1, 0, 0.01);
|
||||
Assert.assertEquals(0.5 / FastMath.sqrt(t001.getValue()), sqrt.value(t001).getPartialDerivative(1), 1.6);
|
||||
DerivativeStructure t01 = new DerivativeStructure(1, 1, 0, 0.1);
|
||||
Assert.assertEquals(0.5 / FastMath.sqrt(t01.getValue()), sqrt.value(t01).getPartialDerivative(1), 7.0e-3);
|
||||
DerivativeStructure t03 = new DerivativeStructure(1, 1, 0, 0.3);
|
||||
Assert.assertEquals(0.5 / FastMath.sqrt(t03.getValue()), sqrt.value(t03).getPartialDerivative(1), 2.1e-7);
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testVectorFunction() {
|
||||
|
||||
|
|
@ -219,12 +300,12 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
DerivativeStructure[] y = f.value(new DerivativeStructure(1, 2, 0, x));
|
||||
double cos = FastMath.cos(x);
|
||||
double sin = FastMath.sin(x);
|
||||
Assert.assertEquals(cos, y[0].getValue(), 2.0e-16);
|
||||
Assert.assertEquals(sin, y[1].getValue(), 2.0e-16);
|
||||
Assert.assertEquals(-sin, y[0].getPartialDerivative(1), 5.0e-14);
|
||||
Assert.assertEquals( cos, y[1].getPartialDerivative(1), 5.0e-14);
|
||||
Assert.assertEquals(-cos, y[0].getPartialDerivative(2), 6.0e-12);
|
||||
Assert.assertEquals(-sin, y[1].getPartialDerivative(2), 6.0e-12);
|
||||
Assert.assertEquals( cos, y[0].getValue(), 7.0e-16);
|
||||
Assert.assertEquals( sin, y[1].getValue(), 7.0e-16);
|
||||
Assert.assertEquals(-sin, y[0].getPartialDerivative(1), 6.0e-14);
|
||||
Assert.assertEquals( cos, y[1].getPartialDerivative(1), 6.0e-14);
|
||||
Assert.assertEquals(-cos, y[0].getPartialDerivative(2), 2.0e-11);
|
||||
Assert.assertEquals(-sin, y[1].getPartialDerivative(2), 2.0e-11);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -253,7 +334,7 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
double cosh = FastMath.cosh(x);
|
||||
double sinh = FastMath.sinh(x);
|
||||
Assert.assertEquals(cos, y[0][0].getValue(), 7.0e-18);
|
||||
Assert.assertEquals(sin, y[0][1].getValue(), 7.0e-18);
|
||||
Assert.assertEquals(sin, y[0][1].getValue(), 6.0e-17);
|
||||
Assert.assertEquals(cosh, y[1][0].getValue(), 3.0e-16);
|
||||
Assert.assertEquals(sinh, y[1][1].getValue(), 3.0e-16);
|
||||
Assert.assertEquals(-sin, y[0][0].getPartialDerivative(1), 2.0e-14);
|
||||
|
|
@ -276,7 +357,7 @@ public class FiniteDifferencesDifferentiatorTest {
|
|||
UnivariateDifferentiableFunction f =
|
||||
differentiator.differentiate(sine);
|
||||
double[] expectedError = new double[] {
|
||||
1.110e-16, 2.66e-12, 4.803e-9, 5.486e-5
|
||||
6.696e-16, 1.371e-12, 2.007e-8, 1.754e-5
|
||||
};
|
||||
double[] maxError = new double[expectedError.length];
|
||||
for (double x = -2; x < 2; x += 0.1) {
|
||||
|
|
|
|||
Loading…
Reference in New Issue