diff --git a/src/Numerics/RootFinding/Algorithms/Bisection.cs b/src/Numerics/RootFinding/Algorithms/Bisection.cs index 6941f1b2..ca5ca3e4 100644 --- a/src/Numerics/RootFinding/Algorithms/Bisection.cs +++ b/src/Numerics/RootFinding/Algorithms/Bisection.cs @@ -34,8 +34,14 @@ namespace MathNet.Numerics.RootFinding.Algorithms { public static class Bisection { + public static double FindRootExpand(Func f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-5, double expandFactor = 1.6, int maxExpandIteratons = 100) + { + Bracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons); + return FindRoot(f, guessLowerBound, guessUpperBound, accuracy); + } + /// Find a solution of the equation f(x)=0. - public static double FindRoot(Func f, double lowerBound, double upperBound, double fTolerance = 1e-5, double xTolerance = 1e-5, double lowerExpansionFactor = -1.0, double upperExpansionFactor = -1.0, int maxExpansionSteps = 10) + public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy = 1e-5) { double fmin = f(lowerBound); double fmax = f(upperBound); @@ -48,32 +54,12 @@ namespace MathNet.Numerics.RootFinding.Algorithms ValidateEvaluation(fmin, lowerBound); ValidateEvaluation(fmax, upperBound); - if (Math.Sign(fmin) == Math.Sign(fmax) && lowerExpansionFactor <= 1.0 && upperExpansionFactor <= 1.0) - throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case."); - - int expansionSteps = 0; - while (Math.Sign(fmin) == Math.Sign(fmax) && expansionSteps < maxExpansionSteps) + if (Math.Sign(fmin) == Math.Sign(fmax)) { - double range = upperBound - lowerBound; - if (upperExpansionFactor <= 0.0 || (lowerExpansionFactor > 0.0 && Math.Abs(fmin) < Math.Abs(fmax))) - { - lowerBound = upperBound - lowerExpansionFactor * range; - fmin = f(lowerBound); - ValidateEvaluation(fmin, lowerBound); - } - else - { - upperBound = lowerBound + upperExpansionFactor * range; - fmax = f(upperBound); - ValidateEvaluation(fmax, upperBound); - } - expansionSteps += 1; + throw new NonConvergenceException("Bounds do not necessarily span a root."); } - if (expansionSteps == maxExpansionSteps) - throw new NonConvergenceException(); - - while (Math.Abs(fmax - fmin) > 0.5 * fTolerance || Math.Abs(upperBound - lowerBound) > 0.5 * xTolerance) + while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision) { double midpoint = 0.5*(upperBound + lowerBound); double midval = f(midpoint); diff --git a/src/Numerics/RootFinding/Algorithms/Brent.cs b/src/Numerics/RootFinding/Algorithms/Brent.cs index 5ef3aa29..df68e819 100644 --- a/src/Numerics/RootFinding/Algorithms/Brent.cs +++ b/src/Numerics/RootFinding/Algorithms/Brent.cs @@ -77,7 +77,7 @@ namespace MathNet.Numerics.RootFinding.Algorithms // convergence check double xAcc1 = 2.0*Precision.DoubleMachinePrecision*Math.Abs(root) + 0.5*accuracy; double xMid = (upperBound - root)/2.0; - if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualWithAbsoluteError(0, froot, accuracy)) + if (Math.Abs(xMid) <= xAcc1 && froot.AlmostEqualWithAbsoluteError(0, froot, accuracy)) { return root; } diff --git a/src/Numerics/RootFinding/Bracketing.cs b/src/Numerics/RootFinding/Bracketing.cs index f90fb8bb..019c4291 100644 --- a/src/Numerics/RootFinding/Bracketing.cs +++ b/src/Numerics/RootFinding/Bracketing.cs @@ -37,21 +37,21 @@ namespace MathNet.Numerics.RootFinding { /// Detect a range containing at least one root. /// The function to detect roots from. - /// Lower value of the range. - /// Upper value of the range + /// Lower value of the range. + /// Upper value of the range /// The growing factor of research. Usually 1.6. /// Maximum number of iterations. Usually 50. /// True if the bracketing operation succeeded, false otherwise. /// This iterative methods stops when two values with opposite signs are found. - public static bool SearchOutward(Func f, ref double xmin, ref double xmax, double factor = 1.6, int maxIterations = 50) + public static bool Expand(Func f, ref double lowerBound, ref double upperBound, double factor = 1.6, int maxIterations = 50) { - if (xmin >= xmax) + if (lowerBound >= upperBound) { - throw new ArgumentOutOfRangeException("xmax", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin")); + throw new ArgumentOutOfRangeException("upperBound", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin")); } - double fmin = f(xmin); - double fmax = f(xmax); + double fmin = f(lowerBound); + double fmax = f(upperBound); for (int i = 0; i < maxIterations; i++) { @@ -62,13 +62,13 @@ namespace MathNet.Numerics.RootFinding if (Math.Abs(fmin) < Math.Abs(fmax)) { - xmin += factor*(xmin - xmax); - fmin = f(xmin); + lowerBound += factor*(lowerBound - upperBound); + fmin = f(lowerBound); } else { - xmax += factor*(xmax - xmin); - fmax = f(xmax); + upperBound += factor*(upperBound - lowerBound); + fmax = f(upperBound); } } diff --git a/src/UnitTests/RootFindingTests/BisectionTest.cs b/src/UnitTests/RootFindingTests/BisectionTest.cs index c99452e9..e86ff19c 100644 --- a/src/UnitTests/RootFindingTests/BisectionTest.cs +++ b/src/UnitTests/RootFindingTests/BisectionTest.cs @@ -40,15 +40,27 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void MultipleRoots() { - var f1 = new Func(x => (x - 3)*(x - 4)); - double r1 = Bisection.FindRoot(f1, 2.1, 3.9, 0.001, 0.001); - Assert.That(Math.Abs(f1(r1)), Is.LessThan(0.001)); - Assert.That(Math.Abs(r1 - 3.0), Is.LessThan(0.001)); + // Roots at -2, 2 + Func f1 = x => x * x - 4; + Assert.AreEqual(0, f1(Bisection.FindRoot(f1, 0, 5, 1e-14))); + Assert.AreEqual(0, f1(Bisection.FindRootExpand(f1, 3, 4, 1e-14))); + Assert.AreEqual(-2, Bisection.FindRoot(f1, -5, -1, 1e-14)); + Assert.AreEqual(2, Bisection.FindRoot(f1, 1, 4, 1e-14)); - var f2 = new Func(x => (x - 3)*(x - 4)); - double r2 = Bisection.FindRoot(f1, 2.1, 3.4, 0.001, 0.001); - Assert.That(Math.Abs(f2(r2)), Is.LessThan(0.001)); - Assert.That(Math.Abs(r2 - 3.0), Is.LessThan(0.001)); + // Roots at 3, 4 + Func f2 = x => (x - 3) * (x - 4); + Assert.AreEqual(0, f2(Bisection.FindRoot(f2, 3.5, 5, 1e-14)), 1e-14); + Assert.AreEqual(3, Bisection.FindRoot(f2, -5, 3.5, 1e-14)); + Assert.AreEqual(4, Bisection.FindRoot(f2, 3.2, 5, 1e-14)); + Assert.AreEqual(3, Bisection.FindRoot(f2, 2.1, 3.9, 0.001), 0.001); + Assert.AreEqual(3, Bisection.FindRoot(f2, 2.1, 3.4, 0.001), 0.001); + } + + [Test] + public void NoRoot() + { + Func f1 = x => x * x + 4; + Assert.Throws(() => Bisection.FindRoot(f1, -5, 5, 1e-14)); } } } diff --git a/src/UnitTests/RootFindingTests/BrentTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs index d150a9b1..6e94d985 100644 --- a/src/UnitTests/RootFindingTests/BrentTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -51,6 +51,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.AreEqual(0, f2(Brent.FindRoot(f2, -5, 5, 1e-14, 100))); Assert.AreEqual(3, Brent.FindRoot(f2, -5, 3.5, 1e-14, 100)); Assert.AreEqual(4, Brent.FindRoot(f2, 3.2, 5, 1e-14, 100)); + Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.9, 0.001, 50), 0.001); + Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001); + } + + [Test] + public void NoRoot() + { + Func f1 = x => x * x + 4; + Assert.Throws(() => Brent.FindRoot(f1, -5, 5, 1e-14, 50)); } } }