RootFinding: protect newton-raphson better against singularities, mor…

…e tests
ArtyomBaranovskiy · May 4, 2013 · ccf89f8 · ccf89f8
1 parent 7ed6c6f
commit ccf89f8
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 2 deletions.
diff --git a/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs b/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
@@ -129,7 +129,7 @@ public static bool TryFindRoot(Func<double, double> f, Func<double, double> df,
                     continue;
                 }
 
-                if (Math.Abs(step) < accuracy)
+                if (Math.Abs(step) < accuracy && Math.Abs(fx) < accuracy)
                 {
                     return true;
                 }
@@ -149,7 +149,7 @@ public static bool TryFindRoot(Func<double, double> f, Func<double, double> df,
                     lowerBound = root;
                     fmin = fx;
                 }
-                else if (Math.Sign(fmin) != Math.Sign(fmax))
+                else if (Math.Sign(fmin) != Math.Sign(fmax) && Math.Abs(fx) < accuracy)
                 {
                     return true;
                 }

diff --git a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
@@ -91,6 +91,13 @@ public void Pole()
             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));
+
+            Func<double, double> f4 = x => 1/(2 - x) - x + 6;
+            Func<double, double> df4 = x => 1/(x*x - 4*x + 4) - 1;
+            Assert.AreEqual(4 + Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f4, df4, 5, 6, 1e-14, 100, 20), 1e-14);
+            Assert.AreEqual(4 - Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f4, df4, 2.01, 3, 1e-14, 100, 20));
+            Assert.AreEqual(4 - Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f4, df4, 2.01, 5, 1e-14, 100, 20));
+            Assert.AreEqual(4 - Math.Sqrt(3), HybridNewtonRaphson.FindRoot(f4, df4, -2, 4, 1e-14, 100, 20));
         }
 
         [Test]