RootFinding: applying common pattern also to bisection algorithm

src/Numerics/RootFinding/Bisection.cs

@@ -29,71 +29,109 @@
 // </copyright>
 
 using System;
-using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.RootFinding
 {
     /// <summary>
-    /// Bisection root-finding algorithm.
+    /// Bisection root-finding algorithm without any recovery measures in case of lacking bracketing.
     /// </summary>
     public static class Bisection
     {
         /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <param name="f">The function to find roots from.</param>
+        /// <param name="guessLowerBound">Guess for the low value of the range where the root is supposed to be. Will be expanded if needed.</param>
+        /// <param name="guessUpperBound">Guess for the  high value of the range where the root is supposed to be. Will be expanded if needed.</param>
+        /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8.</param>
+        /// <param name="maxIterations">Maximum number of iterations. Default 100.</param>
+        /// <param name="expandFactor">Factor at which to expand the bounds, if needed. Default 1.6.</param>
+        /// <param name="maxExpandIteratons">Maximum number of expand iterations. Default 100.</param>
+        /// <returns>Returns the root with the specified accuracy.</returns>
         /// <exception cref="NonConvergenceException"></exception>
-        public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, double expandFactor = 1.6, int maxExpandIteratons = 100)
+        public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, int maxIterations = 100, double expandFactor = 1.6, int maxExpandIteratons = 100)
         {
             ZeroCrossingBracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons);
-            return FindRoot(f, guessLowerBound, guessUpperBound, accuracy);
+            return FindRoot(f, guessLowerBound, guessUpperBound, accuracy, maxIterations);
         }
 
         /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <param name="f">The function to find roots from.</param>
+        /// <param name="lowerBound">The low value of the range where the root is supposed to be.</param>
+        /// <param name="upperBound">The high value of the range where the root is supposed to be.</param>
+        /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8.</param>
+        /// <param name="maxIterations">Maximum number of iterations. Default 100.</param>
+        /// <returns>Returns the root with the specified accuracy.</returns>
         /// <exception cref="NonConvergenceException"></exception>
-        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-8)
+        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
+        {
+            double root;
+            if (TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out root))
+            {
+                return root;
+            }
+            throw new NonConvergenceException("The algorithm has failed or exceeded the number of iterations allowed");
+        }
+
+        /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <param name="f">The function to find roots from.</param>
+        /// <param name="lowerBound">The low value of the range where the root is supposed to be.</param>
+        /// <param name="upperBound">The high value of the range where the root is supposed to be.</param>
+        /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.</param>
+        /// <param name="maxIterations">Maximum number of iterations. Usually 100.</param>
+        /// <param name="root">The root that was found, if any. Undefined if the function returns false.</param>
+        /// <returns>True if a root with the specified accuracy was found, else false.</returns>
+        public static bool TryFindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy, int maxIterations, out double root)
         {
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
 
+            // already there?
             if (Math.Abs(fmin) < accuracy)
             {
-                return lowerBound;
+                root = lowerBound;
+                return true;
             }
             if (Math.Abs(fmax) < accuracy)
             {
-                return upperBound;
+                root = upperBound;
+                return true;
             }
 
+            root = 0.5*(lowerBound + upperBound);
+
+            // bad bracketing?
             if (Math.Sign(fmin) == Math.Sign(fmax))
             {
-                throw new NonConvergenceException(Resources.RootMustBeBracketedByBounds);
+                return false;
             }
 
-            while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision)
+            for (int i = 0; i <= maxIterations; i++)
             {
-                double midpoint = 0.5*(upperBound + lowerBound);
+                if (Math.Abs(fmax - fmin) < 0.5*accuracy && upperBound.AlmostEqual(lowerBound))
-                if ((midpoint == lowerBound) || (midpoint == upperBound))
                 {
-                    return midpoint;
+                    return true;
                 }
 
-                double midval = f(midpoint);
+                double midval = f(root);
 
                 if (Math.Sign(midval) == Math.Sign(fmin))
                 {
-                    lowerBound = midpoint;
+                    lowerBound = root;
                     fmin = midval;
                 }
                 else if (Math.Sign(midval) == Math.Sign(fmax))
                 {
-                    upperBound = midpoint;
+                    upperBound = root;
                     fmax = midval;
                 }
                 else
                 {
-                    return midpoint;
+                    return true;
                 }
+
+                root = 0.5*(lowerBound + upperBound);
             }
 
-            return 0.5*(lowerBound + upperBound);
+            return false;
         }
     }
 }

src/UnitTests/RootFindingTests/BisectionTest.cs

@@ -78,7 +78,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
         public void Oneeq1()
         {
             Func<double, double> f1 = z => 8 * Math.Pow((4 - z) * z, 2) / (Math.Pow(6 - 3 * z, 2) * (2 - z)) - 0.186;
-            double x = Bisection.FindRoot(f1, 0.1, 0.9);
+            double x = Bisection.FindRoot(f1, 0.1, 0.9, accuracy: 1e-9, maxIterations: 80);
             Assert.AreEqual(0.277759543089215, x, 1e-9);
             Assert.AreEqual(0, f1(x), 1e-16);
         }