diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index f3f830c63..2eee6c0bb 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -110,8 +110,7 @@ - - + diff --git a/src/Numerics/RootFinding/BrentRootFinder.cs b/src/Numerics/RootFinding/FindRoots.cs similarity index 70% rename from src/Numerics/RootFinding/BrentRootFinder.cs rename to src/Numerics/RootFinding/FindRoots.cs index 6f5b9096b..adf93d844 100644 --- a/src/Numerics/RootFinding/BrentRootFinder.cs +++ b/src/Numerics/RootFinding/FindRoots.cs @@ -3,39 +3,32 @@ namespace MathNet.Numerics.RootFinding { - public class BrentRootFinder : RootFinder + public static class FindRoots { - public BrentRootFinder() : base() + /// 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. + /// Returns the root with the specified accuracy. + /// + /// 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 double BrentMethod(Func f, double xmin, double xmax, double accuracy = 1e-8, int maxIterations = 100) { - } - public BrentRootFinder(int numIters, double accuracy) : base(numIters, accuracy) - { - } - - protected override double Find() - { - /* The implementation of the algorithm was inspired by - Press, Teukolsky, Vetterling, and Flannery, - "Numerical Recipes in C", 2nd edition, Cambridge - University Press - */ - double min1, min2; double p, q, r, s, xAcc1, xMid = 0; double d = 0.0, e = 0.0; // set up - double xmin = XMin; - double fxmin = Func(XMin); - double xmax = XMax; - double fxmax = Func(XMax); - + double fxmin = f(xmin); + double fxmax = f(xmax); double root = xmax; double froot = fxmax; - // solve - int i = 0; - for (; i <= Iterations; i++) + for (int i = 0; i <= maxIterations; i++) { if (Math.Sign(froot) == Math.Sign(fxmax)) { @@ -54,7 +47,7 @@ University Press fxmax = fxmin; } // Convergence check - xAcc1 = 2.0 * DoubleAccuracy * Math.Abs(root) + 0.5 * Accuracy; + xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy; xMid = (xmax - root) / 2.0; if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0)) { @@ -105,11 +98,11 @@ University Press root += d; else root += Sign(xAcc1, xMid); - froot = Func(root); + froot = f(root); } // The algorithm has exceeded the number of iterations allowed - throw new RootFindingException(Resources.AccuracyNotReached, i, XMin, XMax, Math.Abs(xMid)); + throw new RootFindingException(Resources.AccuracyNotReached, maxIterations, xmin, xmax, Math.Abs(xMid)); } /// Helper method useful for preventing rounding errors. diff --git a/src/Numerics/RootFinding/RootFinder.cs b/src/Numerics/RootFinding/RootFinder.cs deleted file mode 100644 index 71c3e8f06..000000000 --- a/src/Numerics/RootFinding/RootFinder.cs +++ /dev/null @@ -1,60 +0,0 @@ -using System; - -namespace MathNet.Numerics.RootFinding -{ - public abstract class RootFinder - { - protected const double DoubleAccuracy = 9.99200722162641E-16; - private const int DefaultMaxIterations = 30; - private const double DefaultAccuracy = 1e-8; - - int _maxNumIters; - double _xmin = double.MinValue; - double _xmax = double.MaxValue; - Func _func; - - public RootFinder() : this(DefaultMaxIterations, DefaultAccuracy) - { - } - - public RootFinder(int numIters, double accuracy) - { - _maxNumIters = numIters; - Accuracy = accuracy; - } - - protected double XMin { get { return _xmin; } } - protected double XMax { get { return _xmax; } } - - public Func Func - { - get { return _func; } - set { _func = value; } - } - - public double Accuracy { get; set; } - - public int Iterations - { - set - { - if (value <= 0) throw new ArgumentOutOfRangeException(); - _maxNumIters = value; - } - protected get { return _maxNumIters; } - } - - /// Prototype algorithm for solving the equation f(x)=0. - /// 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. - /// Returns the root with the specified accuracy. - public virtual double Solve(double x1, double x2) - { - _xmin = x1; - _xmax = x2; - return Find(); - } - - protected abstract double Find(); - } -} diff --git a/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs similarity index 50% rename from src/UnitTests/RootFindingTests/BrentRootFindingTest.cs rename to src/UnitTests/RootFindingTests/BrentTest.cs index 4a264a91b..2956caadc 100644 --- a/src/UnitTests/RootFindingTests/BrentRootFindingTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -4,14 +4,13 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests { [TestFixture] - public class BrentRootFindingTest + public class BrentTest { [Test] public void MultipleRoots() { - var solver = new BrentRootFinder(100, 1e-14) {Func = x => x*x - 4}; - double root = solver.Solve(-5, 5); - Assert.AreEqual(0, solver.Func(root)); + double root = FindRoots.BrentMethod(x => x*x - 4, -5, 5, 1e-14, 100); + Assert.AreEqual(0, root*root - 4); } } } diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 916dc10fc..7c771f552 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -770,7 +770,7 @@ - +