RootFinding: default values, newton-raphson initial guess optional (separate function)

src/Numerics/RootFinding/Brent.cs

@@ -42,11 +42,11 @@ namespace MathNet.Numerics.RootFinding
         /// <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="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, int maxIterations)
+        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))

src/Numerics/RootFinding/NewtonRaphson.cs

@@ -42,14 +42,33 @@ namespace MathNet.Numerics.RootFinding
         /// <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="initialGuess">Initial guess of the root.</param>
         /// <param name="lowerBound">The low value of the range where the root is supposed to be. Aborts if it leaves the interval.</param>
         /// <param name="upperBound">The high value of the range where the root is supposed to be. Aborts if it leaves the interval.</param>
-        /// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Example: 1e-14.</param>
-        /// <param name="maxIterations">Maximum number of iterations. Example: 100.</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, Func<double, double> df, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
+        {
+            double root;
+            if (TryFindRoot(f, df, 0.5 * (lowerBound + upperBound), lowerBound, upperBound, accuracy, maxIterations, out root))
+            {
+                return root;
+            }
+            throw new NonConvergenceException("The algorithm failed or has exceeded the number of iterations allowed. Consider to use RobustNewtonRaphson instead.");
+        }
+
+        /// <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="initialGuess">Initial guess of the root.</param>
+        /// <param name="lowerBound">The low value of the range where the root is supposed to be. Aborts if it leaves the interval. Default MinValue.</param>
+        /// <param name="upperBound">The high value of the range where the root is supposed to be. Aborts if it leaves the interval. Default MaxValue.</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, Func<double, double> df, double initialGuess, double lowerBound, double upperBound, double accuracy, int maxIterations)
+        public static double FindRootNearGuess(Func<double, double> f, Func<double, double> df, double initialGuess, double lowerBound = double.MinValue, double upperBound = double.MaxValue, double accuracy = 1e-8, int maxIterations = 100)
         {
             double root;
             if (TryFindRoot(f, df, initialGuess, lowerBound, upperBound, accuracy, maxIterations, out root))

src/Numerics/RootFinding/RobustNewtonRaphson.cs

@@ -43,12 +43,12 @@ namespace MathNet.Numerics.RootFinding
         /// <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. Example: 1e-14.</param>
-        /// <param name="maxIterations">Maximum number of iterations. Example: 100.</param>
-        /// <param name="subdivision">How many parts an interval should be split into for zero crossing scanning in case of lacking bracketing. Example: 20.</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="subdivision">How many parts an interval should be split into for zero crossing scanning in case of lacking bracketing. Default 20.</param>
         /// <returns>Returns the root with the specified accuracy.</returns>
         /// <exception cref="NonConvergenceException"></exception>
-        public static double FindRoot(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 = 1e-8, int maxIterations = 100, int subdivision = 20)
         {
             double root;
             if (TryFindRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, subdivision, out root))

src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs

@@ -43,21 +43,44 @@ 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(NewtonRaphson.FindRoot(f1, df1, 0.5, -5, 6, 1e-14, 100)));
-            Assert.AreEqual(-2, NewtonRaphson.FindRoot(f1, df1, -3.0, -5, -1, 1e-14, 100));
-            Assert.AreEqual(2, NewtonRaphson.FindRoot(f1, df1, 2.5, 1, 4, 1e-14, 100));
-            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 0.6, -5, 6, 1e-14, 100)));
-            Assert.AreEqual(-2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -3, -5, -1, 1e-14, 100));
-            Assert.AreEqual(2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 2.5, 1, 4, 1e-14, 100));
+            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -5, 6, 1e-14)));
+            Assert.AreEqual(-2, NewtonRaphson.FindRoot(f1, df1, -5, -1, 1e-14));
+            Assert.AreEqual(2, NewtonRaphson.FindRoot(f1, df1, 1, 4, 1e-14));
+            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, 6, 1e-14)));
+            Assert.AreEqual(-2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14));
+            Assert.AreEqual(2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14));
 
             // 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(NewtonRaphson.FindRoot(f2, df2, 0.5, -5, 6, 1e-14, 100)));
-            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, -0.75, -5, 3.5, 1e-14, 100));
-            Assert.AreEqual(4, NewtonRaphson.FindRoot(f2, df2, 4.1, 3.2, 5, 1e-14, 100));
-            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 3, 2.1, 3.9, 0.001, 50), 0.001);
-            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.75, 2.1, 3.4, 0.001, 50), 0.001);
+            Assert.AreEqual(0, f2(NewtonRaphson.FindRoot(f2, df2, -5, 6, 1e-14)));
+            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, -5, 3.5, 1e-14));
+            Assert.AreEqual(4, NewtonRaphson.FindRoot(f2, df2, 3.2, 5, 1e-14));
+            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.1, 3.9, 0.001, 50), 0.001);
+            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.1, 3.4, 0.001, 50), 0.001);
+        }
+
+        [Test]
+        public void MultipleRootsNearGuess()
+        {
+            // Roots at -2, 2
+            Func<double, double> f1 = x => x * x - 4;
+            Func<double, double> df1 = x => 2 * x;
+            Assert.AreEqual(0, f1(NewtonRaphson.FindRootNearGuess(f1, df1, 0.5, accuracy:1e-14)));
+            Assert.AreEqual(-2, NewtonRaphson.FindRootNearGuess(f1, df1, -3.0, accuracy: 1e-14));
+            Assert.AreEqual(2, NewtonRaphson.FindRootNearGuess(f1, df1, 2.5, accuracy: 1e-14));
+            Assert.AreEqual(0, f1(NewtonRaphson.FindRootNearGuess(x => -f1(x), x => -df1(x), 0.6, accuracy: 1e-14)));
+            Assert.AreEqual(-2, NewtonRaphson.FindRootNearGuess(x => -f1(x), x => -df1(x), -3, accuracy: 1e-14));
+            Assert.AreEqual(2, NewtonRaphson.FindRootNearGuess(x => -f1(x), x => -df1(x), 2.5, accuracy: 1e-14));
+
+            // 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(NewtonRaphson.FindRootNearGuess(f2, df2, 0.5, accuracy: 1e-14)));
+            Assert.AreEqual(3, NewtonRaphson.FindRootNearGuess(f2, df2, -0.75, accuracy: 1e-14));
+            Assert.AreEqual(4, NewtonRaphson.FindRootNearGuess(f2, df2, 4.1, accuracy: 1e-14));
+            Assert.AreEqual(3, NewtonRaphson.FindRootNearGuess(f2, df2, 3, accuracy: 0.001, maxIterations: 50), 0.001);
+            Assert.AreEqual(3, NewtonRaphson.FindRootNearGuess(f2, df2, 2.75, accuracy: 0.001, maxIterations: 50), 0.001);
         }
 
         [Test]
@@ -65,7 +88,7 @@ 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(NewtonRaphson.FindRoot(f1, df1, 0.5, -5, 6, 1e-14, 100)));
+            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -5, 6, 1e-14)));
             //Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, 1, -2, 4, 1e-14, 100)));
         }
 
@@ -75,15 +98,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
             // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5
             Func<double, double> f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3);
             Func<double, double> df1 = x => Evaluate.Polynomial(x, 5, 8, 9);
-            Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -1.5, -2, -1, 1e-10, 100), 1e-6);
-            Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, 0, -5, 5, 1e-10, 100), 1e-6);
+            Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10), 1e-6);
+            Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10), 1e-6);
 
             // real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50
             Func<double, double> f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2);
             Func<double, double> df2 = x => Evaluate.Polynomial(x, -50, 8, 6);
-            Assert.AreEqual(-6.1466562197069, NewtonRaphson.FindRoot(f2, df2, -6.5, -8, -5, 1e-10, 100), 1e-6);
-            Assert.AreEqual(0.12124737195841, NewtonRaphson.FindRoot(f2, df2, 0, -1, 1, 1e-10, 100), 1e-6);
-            Assert.AreEqual(4.0254088477485, NewtonRaphson.FindRoot(f2, df2, 4, 3, 5, 1e-10, 100), 1e-6);
+            Assert.AreEqual(-6.1466562197069, NewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10), 1e-6);
+            Assert.AreEqual(0.12124737195841, NewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10), 1e-6);
+            Assert.AreEqual(4.0254088477485, NewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10), 1e-6);
         }
 
         [Test]
@@ -91,7 +114,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
         {
             Func<double, double> f1 = x => x * x + 4;
             Func<double, double> df1 = x => 2 * x;
-            Assert.Throws<NonConvergenceException>(() => NewtonRaphson.FindRoot(f1, df1, 0, -5, 5, 1e-14, 50));
+            Assert.Throws<NonConvergenceException>(() => NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50));
         }
     }
 }