diff --git a/src/Numerics/RootFinding/Bisection.cs b/src/Numerics/RootFinding/Bisection.cs index 729edeb8..4ab9c73a 100644 --- a/src/Numerics/RootFinding/Bisection.cs +++ b/src/Numerics/RootFinding/Bisection.cs @@ -29,71 +29,109 @@ // using System; -using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding { /// - /// Bisection root-finding algorithm. + /// Bisection root-finding algorithm without any recovery measures in case of lacking bracketing. /// public static class Bisection { /// Find a solution of the equation f(x)=0. + /// The function to find roots from. + /// Guess for the low value of the range where the root is supposed to be. Will be expanded if needed. + /// Guess for the high value of the range where the root is supposed to be. Will be expanded if needed. + /// 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. + /// Factor at which to expand the bounds, if needed. Default 1.6. + /// Maximum number of expand iterations. Default 100. + /// Returns the root with the specified accuracy. /// - public static double FindRootExpand(Func f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, double expandFactor = 1.6, int maxExpandIteratons = 100) + public static double FindRootExpand(Func f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, int maxIterations = 100, double expandFactor = 1.6, int maxExpandIteratons = 100) { ZeroCrossingBracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons); - return FindRoot(f, guessLowerBound, guessUpperBound, accuracy); + return FindRoot(f, guessLowerBound, guessUpperBound, accuracy, maxIterations); } /// 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. 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 = 1e-8) + 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)) + { + return root; + } + throw new NonConvergenceException("The algorithm has failed or 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. + 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); + // already there? if (Math.Abs(fmin) < accuracy) { - return lowerBound; + root = lowerBound; + return true; } if (Math.Abs(fmax) < accuracy) { - return upperBound; + root = upperBound; + return true; } + root = 0.5*(lowerBound + upperBound); + + // bad bracketing? if (Math.Sign(fmin) == Math.Sign(fmax)) { - throw new NonConvergenceException(Resources.RootMustBeBracketedByBounds); + return false; } - while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision) + for (int i = 0; i <= maxIterations; i++) { - double midpoint = 0.5*(upperBound + lowerBound); - if ((midpoint == lowerBound) || (midpoint == upperBound)) + if (Math.Abs(fmax - fmin) < 0.5*accuracy && upperBound.AlmostEqual(lowerBound)) { - return midpoint; + return true; } - double midval = f(midpoint); + double midval = f(root); if (Math.Sign(midval) == Math.Sign(fmin)) { - lowerBound = midpoint; + lowerBound = root; fmin = midval; } else if (Math.Sign(midval) == Math.Sign(fmax)) { - upperBound = midpoint; + upperBound = root; fmax = midval; } else { - return midpoint; + return true; } + + root = 0.5*(lowerBound + upperBound); } - return 0.5*(lowerBound + upperBound); + return false; } } } diff --git a/src/UnitTests/RootFindingTests/BisectionTest.cs b/src/UnitTests/RootFindingTests/BisectionTest.cs index aced22cc..c60cdb5c 100644 --- a/src/UnitTests/RootFindingTests/BisectionTest.cs +++ b/src/UnitTests/RootFindingTests/BisectionTest.cs @@ -78,7 +78,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests public void Oneeq1() { Func f1 = z => 8 * Math.Pow((4 - z) * z, 2) / (Math.Pow(6 - 3 * z, 2) * (2 - z)) - 0.186; - double x = Bisection.FindRoot(f1, 0.1, 0.9); + double x = Bisection.FindRoot(f1, 0.1, 0.9, accuracy: 1e-9, maxIterations: 80); Assert.AreEqual(0.277759543089215, x, 1e-9); Assert.AreEqual(0, f1(x), 1e-16); }