diff --git a/src/Numerics/RootFinding/Brent.cs b/src/Numerics/RootFinding/Brent.cs index 00173dd8..291c33e2 100644 --- a/src/Numerics/RootFinding/Brent.cs +++ b/src/Numerics/RootFinding/Brent.cs @@ -42,11 +42,11 @@ namespace MathNet.Numerics.RootFinding /// 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. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8. + /// Maximum number of iterations. Default 100. /// Returns the root with the specified accuracy. /// - public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy, int maxIterations) + public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100) { double root; if (TryFindRoot(f, lowerBound, upperBound, accuracy, maxIterations, out root)) diff --git a/src/Numerics/RootFinding/NewtonRaphson.cs b/src/Numerics/RootFinding/NewtonRaphson.cs index 98d13561..5fa8d3b9 100644 --- a/src/Numerics/RootFinding/NewtonRaphson.cs +++ b/src/Numerics/RootFinding/NewtonRaphson.cs @@ -42,14 +42,33 @@ namespace MathNet.Numerics.RootFinding /// 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. - /// Initial guess of the root. /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval. /// The high value of the range where the root is supposed to be. Aborts if it leaves the interval. - /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Example: 1e-14. - /// Maximum number of iterations. Example: 100. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8. + /// Maximum number of iterations. Default 100. + /// Returns the root with the specified accuracy. + /// + public static double FindRoot(Func f, Func 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."); + } + + /// 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. + /// Initial guess of the root. + /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval. Default MinValue. + /// The high value of the range where the root is supposed to be. Aborts if it leaves the interval. Default MaxValue. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8. + /// Maximum number of iterations. Default 100. /// Returns the root with the specified accuracy. /// - public static double FindRoot(Func f, Func df, double initialGuess, double lowerBound, double upperBound, double accuracy, int maxIterations) + public static double FindRootNearGuess(Func f, Func 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)) diff --git a/src/Numerics/RootFinding/RobustNewtonRaphson.cs b/src/Numerics/RootFinding/RobustNewtonRaphson.cs index 13424418..18520334 100644 --- a/src/Numerics/RootFinding/RobustNewtonRaphson.cs +++ b/src/Numerics/RootFinding/RobustNewtonRaphson.cs @@ -43,12 +43,12 @@ namespace MathNet.Numerics.RootFinding /// 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. Example: 1e-14. - /// Maximum number of iterations. Example: 100. - /// How many parts an interval should be split into for zero crossing scanning in case of lacking bracketing. Example: 20. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8. + /// Maximum number of iterations. Default 100. + /// How many parts an interval should be split into for zero crossing scanning in case of lacking bracketing. Default 20. /// Returns the root with the specified accuracy. /// - public static double FindRoot(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 = 1e-8, int maxIterations = 100, int subdivision = 20) { double root; if (TryFindRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, subdivision, out root)) diff --git a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs index 4b841938..6024cae9 100644 --- a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs +++ b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs @@ -43,21 +43,44 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests // Roots at -2, 2 Func f1 = x => x * x - 4; Func 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 f2 = x => (x - 3) * (x - 4); Func 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 f1 = x => x * x - 4; + Func 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 f2 = x => (x - 3) * (x - 4); + Func 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 f1 = x => x * x * x - 2 * x + 2; Func 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 f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3); Func 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 f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2); Func 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 f1 = x => x * x + 4; Func df1 = x => 2 * x; - Assert.Throws(() => NewtonRaphson.FindRoot(f1, df1, 0, -5, 5, 1e-14, 50)); + Assert.Throws(() => NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50)); } } }