commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1416249 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java test/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiatorTest.java
Date Sun, 02 Dec 2012 20:23:51 GMT
Author: luc
Date: Sun Dec  2 20:23:49 2012
New Revision: 1416249

URL: http://svn.apache.org/viewvc?rev=1416249&view=rev
Log:
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.

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiatorTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java?rev=1416249&r1=1416248&r2=1416249&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiator.java
Sun Dec  2 20:23:49 2012
@@ -25,6 +25,7 @@ import org.apache.commons.math3.exceptio
 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.exceptio
  * <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 FiniteDifferencesDifferenti
     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 FiniteDifferencesDifferenti
      */
     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 &lt; t &lt; 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 FiniteDifferencesDifferenti
         }
         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 FiniteDifferencesDifferenti
     }
 
     /**
-     * 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 FiniteDifferencesDifferenti
         }
 
         // 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 FiniteDifferencesDifferenti
                     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 FiniteDifferencesDifferenti
                     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 FiniteDifferencesDifferenti
                 // 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 FiniteDifferencesDifferenti
                     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 FiniteDifferencesDifferenti
                 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]);
                     }
                 }
 

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiatorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiatorTest.java?rev=1416249&r1=1416248&r2=1416249&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiatorTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/FiniteDifferencesDifferentiatorTest.java
Sun Dec  2 20:23:49 2012
@@ -87,9 +87,9 @@ public class FiniteDifferencesDifferenti
                 });
         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 FiniteDifferencesDifferenti
         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 FiniteDifferencesDifferenti
         // 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 FiniteDifferencesDifferenti
         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 FiniteDifferencesDifferenti
             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 FiniteDifferencesDifferenti
             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 FiniteDifferencesDifferenti
         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) {



Mime
View raw message