Merge pull request #959 from jkalias/robust-newton-raphson-root-on-the-bound

Robust Newton Raphson root on the bound

Author: cdrnet

src/Numerics.Tests/RootFindingTests/RobustNewtonRaphsonTest.cs

@@ -115,6 +115,58 @@ public void Cubic()
             Assert.AreEqual(4.0254088477485, RobustNewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10, 100, 20), 1e-6);
         }
 
+        [Test]
+        public void InfinitelyManyRoots()
+        {
+            // degenerate case with infinitely many roots
+            Func<double, double> f1 = x => 0.0;
+            Func<double, double> df1 = x => 0.0;
+            Assert.AreEqual(-50, RobustNewtonRaphson.FindRoot(f1, df1, -200, 100), 1e-6);
+        }
+
+        [Test]
+        public void InfinitelyManyRootsWithGivenAccuracy()
+        {
+            // degenerate case with infinitely many roots
+            Func<double, double> f1 = x => 1e-10;
+            Func<double, double> df1 = x => 1e-6;
+            Assert.AreEqual(-50, RobustNewtonRaphson.FindRoot(f1, df1, -200, 100, 1e-8), 1e-6);
+        }
+
+        [Test]
+        public void BoundsNearMaxValueNoOverflow()
+        {
+            Func<double, double> f1 = x => 0.0;
+            Func<double, double> df1 = x => 0.0;
+            Assert.AreEqual(double.MaxValue, RobustNewtonRaphson.FindRoot(f1, df1, double.MaxValue - 2, double.MaxValue - 1), 1e-6);
+        }
+
+        [Test]
+        public void BoundsNearMinValueNoUnderflow()
+        {
+            Func<double, double> f1 = x => 0.0;
+            Func<double, double> df1 = x => 0.0;
+            Assert.AreEqual(double.MinValue, RobustNewtonRaphson.FindRoot(f1, df1, double.MinValue + 1, double.MinValue + 2), 1e-6);
+        }
+
+        [Test]
+        public void InfiniteLowerBound()
+        {
+            Func<double, double> f1 = x => 0.0;
+            Func<double, double> df1 = x => 0.0;
+            Assert.That(() => RobustNewtonRaphson.FindRoot(f1, df1, double.NegativeInfinity, 5), Throws.TypeOf<ArgumentOutOfRangeException>());
+            Assert.That(() => RobustNewtonRaphson.FindRoot(f1, df1, double.PositiveInfinity, 5), Throws.TypeOf<ArgumentOutOfRangeException>());
+        }
+
+        [Test]
+        public void InfiniteUpperBound()
+        {
+            Func<double, double> f1 = x => 0.0;
+            Func<double, double> df1 = x => 0.0;
+            Assert.That(() => RobustNewtonRaphson.FindRoot(f1, df1, -5, double.NegativeInfinity), Throws.TypeOf<ArgumentOutOfRangeException>());
+            Assert.That(() => RobustNewtonRaphson.FindRoot(f1, df1, -5, double.PositiveInfinity), Throws.TypeOf<ArgumentOutOfRangeException>());
+        }
+
         [Test]
         public void NoRoot()
         {

src/Numerics/RootFinding/RobustNewtonRaphson.cs

@@ -60,8 +60,8 @@ public static double FindRoot(Func<double, double> f, Func<double, double> df, d
         /// <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="lowerBound">The low value of the range where the root is supposed to be (finite number).</param>
+        /// <param name="upperBound">The high value of the range where the root is supposed to be (finite number).</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. Must be greater than 0.</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>
@@ -74,6 +74,23 @@ public static bool TryFindRoot(Func<double, double> f, Func<double, double> df,
                 throw new ArgumentOutOfRangeException(nameof(accuracy), "Must be greater than zero.");
             }
 
+            if (double.IsInfinity(lowerBound))
+            {
+                throw new ArgumentOutOfRangeException(nameof(lowerBound), "Must be a finite number.");
+            }
+
+            if (double.IsInfinity(upperBound))
+            {
+                throw new ArgumentOutOfRangeException(nameof(upperBound), "Must be a finite number.");
+            }
+
+            root = lowerBound + 0.5 * (upperBound - lowerBound);
+            double fx = f(root);
+            if (Math.Abs(fx) < accuracy)
+            {
+                return true;
+            }
+
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
 
@@ -89,13 +106,6 @@ public static bool TryFindRoot(Func<double, double> f, Func<double, double> df,
                 return true;
             }
 
-            root = 0.5*(lowerBound + upperBound);
-            double fx = f(root);
-            if (fx == 0.0)
-            {
-                return true;
-            }
-
             double lastStep = Math.Abs(upperBound - lowerBound);
             for (int i = 0; i < maxIterations; i++)
             {