diff --git a/src/Numerics/RootFinding/Algorithms/Brent.cs b/src/Numerics/RootFinding/Algorithms/Brent.cs
index 75f789d8..32f982a2 100644
--- a/src/Numerics/RootFinding/Algorithms/Brent.cs
+++ b/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
///
///
- public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
+ public static double FindRoot(Func 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");
+ }
+
+ /// Find a solution of the equation f(x)=0.
+ /// The function to find roots from.
+ /// The low value of the range where the root is supposed to be.
+ /// The high value of the range where the root is supposed to be.
+ /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.
+ /// Maximum number of iterations. Usually 100.
+ /// The root that was found, if any. Undefined if the function returns false.
+ /// True if a root with the specified accuracy was found, else false.
+ ///
+ /// 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
+ ///
+ public static bool TryFindRoot(Func 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;
}
/// Helper method useful for preventing rounding errors.
diff --git a/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs b/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
index b066833f..3fd9bd82 100644
--- a/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
+++ b/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
@@ -35,12 +35,19 @@ namespace MathNet.Numerics.RootFinding.Algorithms
public static class HybridNewtonRaphson
{
/// Find a solution of the equation f(x)=0.
+ /// The function to find roots from.
+ /// The first derivative of the function to find roots from.
+ /// The low value of the range where the root is supposed to be.
+ /// The high value of the range where the root is supposed to be.
+ /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.
+ /// Maximum number of iterations.
+ /// Returns the root with the specified accuracy.
/// Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.
///
- public static double FindSingleRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision)
+ public static double FindRoot(Func f, Func 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
}
/// Find a solution of the equation f(x)=0.
+ /// The function to find roots from.
+ /// The first derivative of the function to find roots from.
+ /// The low value of the range where the root is supposed to be.
+ /// The high value of the range where the root is supposed to be.
+ /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.
+ /// Maximum number of iterations.
+ /// The root that was found, if any. Undefined if the function returns false.
+ /// True if a root with the specified accuracy was found, else false.
/// Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.
- public static bool TryFindSingleRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
+ public static bool TryFindRoot(Func f, Func 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 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;
}
diff --git a/src/Numerics/RootFinding/FloatingPointRoots.cs b/src/Numerics/RootFinding/FloatingPointRoots.cs
index 4fa14d8c..79324a31 100644
--- a/src/Numerics/RootFinding/FloatingPointRoots.cs
+++ b/src/Numerics/RootFinding/FloatingPointRoots.cs
@@ -43,14 +43,16 @@ namespace MathNet.Numerics.RootFinding
public static double OfFunctionAndDerivative(Func f, Func 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");
}
}
}
diff --git a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
index 5a091f69..808cd76d 100644
--- a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
+++ b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
@@ -44,21 +44,21 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
// Roots at -2, 2
Func f1 = x => x * x - 4;
Func 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 f2 = x => (x - 3) * (x - 4);
Func 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 f1 = x => x * x * x - 2 * x + 2;
Func 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 f1 = x => 1/(x - 2) + 2;
Func 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 f2 = x => -1/(x - 2) + 2;
Func 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 f3 = x => 1/(x - 2) + x + 2;
Func 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 f1 = x => x * x + 4;
Func df1 = x => 2 * x;
- Assert.Throws(() => HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 50, 20));
+ Assert.Throws(() => HybridNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50, 20));
}
}
}