RobustNewtonRaphson returns the mean value of the search range for degenerate cases, in which the function is tautologically 0 everywhere.

Author: jkalias

src/Numerics.Tests/RootFindingTests/RobustNewtonRaphsonTest.cs

@@ -115,6 +115,24 @@ 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 NoRoot()
         {

src/Numerics/RootFinding/RobustNewtonRaphson.cs

@@ -77,22 +77,22 @@ public static bool TryFindRoot(Func<double, double> f, Func<double, double> df,
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
 
-            if (Math.Abs(fmin) < accuracy)
+            root = 0.5 * (lowerBound + upperBound);
+            double fx = f(root);
+            if (Math.Abs(fx) < accuracy)
             {
-                root = lowerBound;
                 return true;
             }
 
-            if (Math.Abs(fmax) < accuracy)
+            if (Math.Abs(fmin) < accuracy)
             {
-                root = upperBound;
+                root = lowerBound;
                 return true;
             }
 
-            root = 0.5*(lowerBound + upperBound);
-            double fx = f(root);
-            if (fx == 0.0)
+            if (Math.Abs(fmax) < accuracy)
             {
+                root = upperBound;
                 return true;
             }