RootFinding: separate bracketing from algorithm, unify error behavior

src/Numerics/RootFinding/Algorithms/Bisection.cs

@@ -34,8 +34,14 @@ namespace MathNet.Numerics.RootFinding.Algorithms
 {
     public static class Bisection
     {
+        public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-5, double expandFactor = 1.6, int maxExpandIteratons = 100)
+        {
+            Bracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons);
+            return FindRoot(f, guessLowerBound, guessUpperBound, accuracy);
+        }
+
         /// <summary>Find a solution of the equation f(x)=0.</summary>
-        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double fTolerance = 1e-5, double xTolerance = 1e-5, double lowerExpansionFactor = -1.0, double upperExpansionFactor = -1.0, int maxExpansionSteps = 10)
+        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-5)
         {
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
@@ -48,32 +54,12 @@ namespace MathNet.Numerics.RootFinding.Algorithms
             ValidateEvaluation(fmin, lowerBound);
             ValidateEvaluation(fmax, upperBound);
 
-            if (Math.Sign(fmin) == Math.Sign(fmax) && lowerExpansionFactor <= 1.0 && upperExpansionFactor <= 1.0)
-                throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case.");
-
-            int expansionSteps = 0;
-            while (Math.Sign(fmin) == Math.Sign(fmax) && expansionSteps < maxExpansionSteps)
+            if (Math.Sign(fmin) == Math.Sign(fmax))
             {
-                double range = upperBound - lowerBound;
-                if (upperExpansionFactor <= 0.0 || (lowerExpansionFactor > 0.0 && Math.Abs(fmin) < Math.Abs(fmax)))
-                {
-                    lowerBound = upperBound - lowerExpansionFactor * range;
-                    fmin = f(lowerBound);
-                    ValidateEvaluation(fmin, lowerBound);
-                }
-                else
-                {
-                    upperBound = lowerBound + upperExpansionFactor * range;
-                    fmax = f(upperBound);
-                    ValidateEvaluation(fmax, upperBound);
-                }
-                expansionSteps += 1;
+                throw new NonConvergenceException("Bounds do not necessarily span a root.");
             }
 
-            if (expansionSteps == maxExpansionSteps)
-                throw new NonConvergenceException();
-
-            while (Math.Abs(fmax - fmin) > 0.5 * fTolerance || Math.Abs(upperBound - lowerBound) > 0.5 * xTolerance)
+            while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision)
             {
                 double midpoint = 0.5*(upperBound + lowerBound);
                 double midval = f(midpoint);

src/Numerics/RootFinding/Algorithms/Brent.cs

@@ -77,7 +77,7 @@ namespace MathNet.Numerics.RootFinding.Algorithms
                 // convergence check
                 double xAcc1 = 2.0*Precision.DoubleMachinePrecision*Math.Abs(root) + 0.5*accuracy;
                 double xMid = (upperBound - root)/2.0;
-                if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualWithAbsoluteError(0, froot, accuracy))
+                if (Math.Abs(xMid) <= xAcc1 && froot.AlmostEqualWithAbsoluteError(0, froot, accuracy))
                 {
                     return root;
                 }

src/Numerics/RootFinding/Bracketing.cs

@@ -37,21 +37,21 @@ namespace MathNet.Numerics.RootFinding
     {
         /// <summary>Detect a range containing at least one root.</summary>
         /// <param name="f">The function to detect roots from.</param>
-        /// <param name="xmin">Lower value of the range.</param>
-        /// <param name="xmax">Upper value of the range</param>
+        /// <param name="lowerBound">Lower value of the range.</param>
+        /// <param name="upperBound">Upper value of the range</param>
         /// <param name="factor">The growing factor of research. Usually 1.6.</param>
         /// <param name="maxIterations">Maximum number of iterations. Usually 50.</param>
         /// <returns>True if the bracketing operation succeeded, false otherwise.</returns>
         /// <remarks>This iterative methods stops when two values with opposite signs are found.</remarks>
-        public static bool SearchOutward(Func<double, double> f, ref double xmin, ref double xmax, double factor = 1.6, int maxIterations = 50)
+        public static bool Expand(Func<double, double> f, ref double lowerBound, ref double upperBound, double factor = 1.6, int maxIterations = 50)
         {
-            if (xmin >= xmax)
+            if (lowerBound >= upperBound)
             {
-                throw new ArgumentOutOfRangeException("xmax", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin"));
+                throw new ArgumentOutOfRangeException("upperBound", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin"));
             }
 
-            double fmin = f(xmin);
-            double fmax = f(xmax);
+            double fmin = f(lowerBound);
+            double fmax = f(upperBound);
 
             for (int i = 0; i < maxIterations; i++)
             {
@@ -62,13 +62,13 @@ namespace MathNet.Numerics.RootFinding
 
                 if (Math.Abs(fmin) < Math.Abs(fmax))
                 {
-                    xmin += factor*(xmin - xmax);
-                    fmin = f(xmin);
+                    lowerBound += factor*(lowerBound - upperBound);
+                    fmin = f(lowerBound);
                 }
                 else
                 {
-                    xmax += factor*(xmax - xmin);
-                    fmax = f(xmax);
+                    upperBound += factor*(upperBound - lowerBound);
+                    fmax = f(upperBound);
                 }
             }
 

src/UnitTests/RootFindingTests/BisectionTest.cs

@@ -40,15 +40,27 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
         [Test]
         public void MultipleRoots()
         {
-            var f1 = new Func<double, double>(x => (x - 3)*(x - 4));
-            double r1 = Bisection.FindRoot(f1, 2.1, 3.9, 0.001, 0.001);
-            Assert.That(Math.Abs(f1(r1)), Is.LessThan(0.001));
-            Assert.That(Math.Abs(r1 - 3.0), Is.LessThan(0.001));
+            // Roots at -2, 2
+            Func<double, double> f1 = x => x * x - 4;
+            Assert.AreEqual(0, f1(Bisection.FindRoot(f1, 0, 5, 1e-14)));
+            Assert.AreEqual(0, f1(Bisection.FindRootExpand(f1, 3, 4, 1e-14)));
+            Assert.AreEqual(-2, Bisection.FindRoot(f1, -5, -1, 1e-14));
+            Assert.AreEqual(2, Bisection.FindRoot(f1, 1, 4, 1e-14));
 
-            var f2 = new Func<double, double>(x => (x - 3)*(x - 4));
-            double r2 = Bisection.FindRoot(f1, 2.1, 3.4, 0.001, 0.001);
-            Assert.That(Math.Abs(f2(r2)), Is.LessThan(0.001));
-            Assert.That(Math.Abs(r2 - 3.0), Is.LessThan(0.001));
+            // Roots at 3, 4
+            Func<double, double> f2 = x => (x - 3) * (x - 4);
+            Assert.AreEqual(0, f2(Bisection.FindRoot(f2, 3.5, 5, 1e-14)), 1e-14);
+            Assert.AreEqual(3, Bisection.FindRoot(f2, -5, 3.5, 1e-14));
+            Assert.AreEqual(4, Bisection.FindRoot(f2, 3.2, 5, 1e-14));
+            Assert.AreEqual(3, Bisection.FindRoot(f2, 2.1, 3.9, 0.001), 0.001);
+            Assert.AreEqual(3, Bisection.FindRoot(f2, 2.1, 3.4, 0.001), 0.001);
+        }
+
+        [Test]
+        public void NoRoot()
+        {
+            Func<double, double> f1 = x => x * x + 4;
+            Assert.Throws<NonConvergenceException>(() => Bisection.FindRoot(f1, -5, 5, 1e-14));
         }
     }
 }

src/UnitTests/RootFindingTests/BrentTest.cs

@@ -51,6 +51,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
             Assert.AreEqual(0, f2(Brent.FindRoot(f2, -5, 5, 1e-14, 100)));
             Assert.AreEqual(3, Brent.FindRoot(f2, -5, 3.5, 1e-14, 100));
             Assert.AreEqual(4, Brent.FindRoot(f2, 3.2, 5, 1e-14, 100));
+            Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.9, 0.001, 50), 0.001);
+            Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001);
+        }
+
+        [Test]
+        public void NoRoot()
+        {
+            Func<double, double> f1 = x => x * x + 4;
+            Assert.Throws<NonConvergenceException>(() => Brent.FindRoot(f1, -5, 5, 1e-14, 50));
         }
     }
 }