RootFinding: Update Brent to TryFindRoot pattern as well

src/Numerics/RootFinding/Algorithms/Brent.cs

@@ -46,14 +46,37 @@ namespace MathNet.Numerics.RootFinding.Algorithms
         /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press
         /// </remarks>
         /// <exception cref="NonConvergenceException"></exception>
-        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
+        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy, int maxIterations)
+        {
+            double root;
+            if (TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out root))
+            {
+                return root;
+            }
+            throw new NonConvergenceException("The algorithm has 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>
+        /// <remarks>
+        /// Algorithm by by Brent, Van Wijngaarden, Dekker et al.
+        /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press
+        /// </remarks>
+        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);
-            double root = upperBound;
             double froot = fmax;
             double d = 0.0, e = 0.0;
 
+            root = upperBound;
+
             for (int i = 0; i <= maxIterations; i++)
             {
                 // adjust bounds
@@ -79,7 +102,7 @@ namespace MathNet.Numerics.RootFinding.Algorithms
                 double xMid = (upperBound - root)/2.0;
                 if (Math.Abs(xMid) <= xAcc1 && froot.AlmostEqualWithAbsoluteError(0, froot, accuracy))
                 {
-                    return root;
+                    return true;
                 }
 
                 if (Math.Abs(e) >= xAcc1 && Math.Abs(fmin) > Math.Abs(froot))
@@ -139,7 +162,7 @@ namespace MathNet.Numerics.RootFinding.Algorithms
                 froot = f(root);
             }
 
-            throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+            return false;
         }
 
         /// <summary>Helper method useful for preventing rounding errors.</summary>

src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs

@@ -35,12 +35,19 @@ namespace MathNet.Numerics.RootFinding.Algorithms
     public static class HybridNewtonRaphson
     {
         /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <param name="f">The function to find roots from.</param>
+        /// <param name="df">The first derivative of 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.</param>
+        /// <returns>Returns the root with the specified accuracy.</returns>
         /// <remarks>Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.</remarks>
         /// <exception cref="NonConvergenceException"></exception>
-        public static double FindSingleRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision)
+        public static double FindRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision)
         {
             double root;
-            if (TryFindSingleRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, subdivision, out root))
+            if (TryFindRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, subdivision, out root))
             {
                 return root;
             }
@@ -48,8 +55,16 @@ namespace MathNet.Numerics.RootFinding.Algorithms
         }
 
         /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <param name="f">The function to find roots from.</param>
+        /// <param name="df">The first derivative of 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.</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>
         /// <remarks>Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.</remarks>
-        public static bool TryFindSingleRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
+        public static bool TryFindRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
         {
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
@@ -148,7 +163,7 @@ namespace MathNet.Numerics.RootFinding.Algorithms
             var zeroCrossings = ZeroCrossingBracketing.FindIntervalsWithin(f, lowerBound, upperBound, subdivision);
             foreach (Tuple<double, double> bounds in zeroCrossings)
             {
-                if (TryFindSingleRoot(f, df, bounds.Item1, bounds.Item2, accuracy, maxIterations, subdivision, out root))
+                if (TryFindRoot(f, df, bounds.Item1, bounds.Item2, accuracy, maxIterations, subdivision, out root))
                 {
                     return true;
                 }

src/Numerics/RootFinding/FloatingPointRoots.cs

@@ -43,14 +43,16 @@ namespace MathNet.Numerics.RootFinding
         public static double OfFunctionAndDerivative(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy = 1e-8)
         {
             double root;
-            if (HybridNewtonRaphson.TryFindSingleRoot(f, df, lowerBound, upperBound, accuracy, 100, 20, out root))
+            if (HybridNewtonRaphson.TryFindRoot(f, df, lowerBound, upperBound, accuracy, 100, 20, out root))
+            {
+                return root;
+            }
+            if (Brent.TryFindRoot(f, lowerBound, upperBound, accuracy, 100, out root))
             {
                 return root;
             }
 
-            return Brent.FindRoot(f, lowerBound, upperBound, accuracy, 100);
-
-            //throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+            throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
         }
     }
 }

src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs

@@ -44,21 +44,21 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
             // Roots at -2, 2
             Func<double, double> f1 = x => x * x - 4;
             Func<double, double> df1 = x => 2 * x;
-            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
-            Assert.AreEqual(-2, HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, -1, 1e-14, 100, 20));
-            Assert.AreEqual(2, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 4, 1e-14, 100, 20));
-            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100, 20)));
-            Assert.AreEqual(-2, HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100, 20));
-            Assert.AreEqual(2, HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100, 20));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(-2, HybridNewtonRaphson.FindRoot(f1, df1, -5, -1, 1e-14, 100, 20));
+            Assert.AreEqual(2, HybridNewtonRaphson.FindRoot(f1, df1, 1, 4, 1e-14, 100, 20));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(-2, HybridNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100, 20));
+            Assert.AreEqual(2, HybridNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100, 20));
 
             // Roots at 3, 4
             Func<double, double> f2 = x => (x - 3) * (x - 4);
             Func<double, double> df2 = x => 2 * x - 7;
-            Assert.AreEqual(0, f2(HybridNewtonRaphson.FindSingleRoot(f2, df2, -5, 5, 1e-14, 100, 20)));
-            Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, -5, 3.5, 1e-14, 100, 20));
-            Assert.AreEqual(4, HybridNewtonRaphson.FindSingleRoot(f2, df2, 3.2, 5, 1e-14, 100, 20));
-            Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2.1, 3.9, 0.001, 50, 20), 0.001);
-            Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2.1, 3.4, 0.001, 50, 20), 0.001);
+            Assert.AreEqual(0, f2(HybridNewtonRaphson.FindRoot(f2, df2, -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(3, HybridNewtonRaphson.FindRoot(f2, df2, -5, 3.5, 1e-14, 100, 20));
+            Assert.AreEqual(4, HybridNewtonRaphson.FindRoot(f2, df2, 3.2, 5, 1e-14, 100, 20));
+            Assert.AreEqual(3, HybridNewtonRaphson.FindRoot(f2, df2, 2.1, 3.9, 0.001, 50, 20), 0.001);
+            Assert.AreEqual(3, HybridNewtonRaphson.FindRoot(f2, df2, 2.1, 3.4, 0.001, 50, 20), 0.001);
         }
 
         [Test]
@@ -66,8 +66,8 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
         {
             Func<double, double> f1 = x => x * x * x - 2 * x + 2;
             Func<double, double> df1 = x => 3 * x * x - 2;
-            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
-            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -2, 4, 1e-14, 100, 20)));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindRoot(f1, df1, -2, 4, 1e-14, 100, 20)));
         }
 
         [Test]
@@ -75,22 +75,22 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
         {
             Func<double, double> f1 = x => 1/(x - 2) + 2;
             Func<double, double> df1 = x => -1/(x*x - 4*x + 4);
-            Assert.AreEqual(1.5, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 2, 1e-14, 100, 20));
-            Assert.AreEqual(1.5, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 6, 1e-14, 100, 20));
+            Assert.AreEqual(1.5, HybridNewtonRaphson.FindRoot(f1, df1, 1, 2, 1e-14, 100, 20));
+            Assert.AreEqual(1.5, HybridNewtonRaphson.FindRoot(f1, df1, 1, 6, 1e-14, 100, 20));
             Assert.AreEqual(1.5, FloatingPointRoots.OfFunctionAndDerivative(f1, df1, 1, 6));
 
             Func<double, double> f2 = x => -1/(x - 2) + 2;
             Func<double, double> df2 = x => 1/(x*x - 4*x + 4);
-            Assert.AreEqual(2.5, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2, 3, 1e-14, 100, 20));
-            Assert.AreEqual(2.5, HybridNewtonRaphson.FindSingleRoot(f2, df2, -2, 3, 1e-14, 100, 20));
+            Assert.AreEqual(2.5, HybridNewtonRaphson.FindRoot(f2, df2, 2, 3, 1e-14, 100, 20));
+            Assert.AreEqual(2.5, HybridNewtonRaphson.FindRoot(f2, df2, -2, 3, 1e-14, 100, 20));
             Assert.AreEqual(2.5, FloatingPointRoots.OfFunctionAndDerivative(f2, df2, -2, 3));
 
             Func<double, double> f3 = x => 1/(x - 2) + x + 2;
             Func<double, double> df3 = x => -1/(x*x - 4*x + 4) + 1;
-            Assert.AreEqual(-Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, -2, -1, 1e-14, 100, 20), 1e-14);
-            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, 1, 1.99, 1e-14, 100, 20));
-            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, -1.5, 1.99, 1e-14, 100, 20));
-            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, 1, 6, 1e-14, 100, 20));
+            Assert.AreEqual(-Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f3, df3, -2, -1, 1e-14, 100, 20), 1e-14);
+            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f3, df3, 1, 1.99, 1e-14, 100, 20));
+            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f3, df3, -1.5, 1.99, 1e-14, 100, 20));
+            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f3, df3, 1, 6, 1e-14, 100, 20));
         }
 
         [Test]
@@ -98,7 +98,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
         {
             Func<double, double> f1 = x => x * x + 4;
             Func<double, double> df1 = x => 2 * x;
-            Assert.Throws<NonConvergenceException>(() => HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 50, 20));
+            Assert.Throws<NonConvergenceException>(() => HybridNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50, 20));
         }
     }
 }