diff --git a/src/Numerics/RootFinders/BrentRootFinder.cs b/src/Numerics/RootFinders/BrentRootFinder.cs new file mode 100644 index 00000000..06c926c6 --- /dev/null +++ b/src/Numerics/RootFinders/BrentRootFinder.cs @@ -0,0 +1,114 @@ +using System; + +namespace MathNet.Numerics.RootFinders +{ + public class BrentRootFinder : RootFinder + { + public BrentRootFinder() : base() + { + } + 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 root = xmax; + double froot = fxmax; + + // solve + int i = 0; + for (; i <= Iterations; i++) + { + if (Math.Sign(froot) == Math.Sign(fxmax)) + { + // Rename xMin_, root_, xMax_ and adjust bounds + xmax = xmin; + fxmax = fxmin; + e = d = root - xmin; + } + if (Math.Abs(fxmax) < Math.Abs(froot)) + { + xmin = root; + root = xmax; + xmax = xmin; + fxmin = froot; + froot = fxmax; + fxmax = fxmin; + } + // Convergence check + xAcc1 = 2.0 * DOUBLE_ACCURACY * Math.Abs(root) + 0.5 * Accuracy; + xMid = (xmax - root) / 2.0; + if (Math.Abs(xMid) <= xAcc1 || Close(froot, 0.0)) + { + return root; + } + if (Math.Abs(e) >= xAcc1 && + Math.Abs(fxmin) > Math.Abs(froot)) + { + + // Attempt inverse quadratic interpolation + s = froot / fxmin; + if (Close(xmin, xmax)) + { + p = 2.0 * xMid * s; + q = 1.0 - s; + } + else + { + q = fxmin / fxmax; + r = froot / fxmax; + p = s * (2.0 * xMid * q * (q - r) - (root - xmin) * (r - 1.0)); + q = (q - 1.0) * (r - 1.0) * (s - 1.0); + } + if (p > 0.0) q = -q; // Check whether in bounds + p = Math.Abs(p); + min1 = 3.0 * xMid * q - Math.Abs(xAcc1 * q); + min2 = Math.Abs(e * q); + if (2.0 * p < Math.Min(min1, min2)) + { + e = d; // Accept interpolation + d = p / q; + } + else + { + d = xMid; // Interpolation failed, use bisection + e = d; + } + } + else + { + // Bounds decreasing too slowly, use bisection + d = xMid; + e = d; + } + xmin = root; + fxmin = froot; + if (Math.Abs(d) > xAcc1) + root += d; + else + root += Sign(xAcc1, xMid); + froot = Func(root); + } + + // The algorithm has exceeded the number of iterations allowed + throw new RootFinderException(ACCURACY_NOT_REACHED, i, new Range(XMin, XMax), Math.Abs(xMid)); + } + } +} diff --git a/src/Numerics/RootFinders/RootFinder.cs b/src/Numerics/RootFinders/RootFinder.cs new file mode 100644 index 00000000..4ee24498 --- /dev/null +++ b/src/Numerics/RootFinders/RootFinder.cs @@ -0,0 +1,175 @@ +using System; + +namespace MathNet.Numerics.RootFinders +{ + public struct Range + { + double Min, Max; + + public Range(double min, double max) + { + Min = min; Max = max; + } + } + public class RootFinderException : Exception + { + private int m_Iteration; + private Range m_Range; + private double m_Accuracy; + + public RootFinderException(string message, int iteration, Range range, double accuracy) + : base(message) + { + m_Iteration = iteration; + m_Range = range; + m_Accuracy = accuracy; + } + + public int Iteration + { + get { return m_Iteration; } + set { m_Iteration = value; } + } + + public Range Range + { + get { return m_Range; } + set { m_Range = value; } + } + + public double Accuracy { set; get; } + } + public abstract class RootFinder + { + + protected const string INVALID_RANGE="Invalid range while finding root"; + protected const string ACCURACY_NOT_REACHED = "The accuracy couldn't be reached with the specified number of iterations"; + protected const string ROOT_NOT_FOUND = "The algorithm ended without root in the range"; + protected const string ROOT_NOT_BRACKETED = "The algorithm could not start because the root seemed not to be bracketed"; + protected const string INVALID_ALGORITHM = "This algorithm is not able to solve this equation"; + protected const double DOUBLE_ACCURACY = 9.99200722162641E-16; + private const int DEFAULT_MAX_ITERATIONS = 30; + private const double DEFAULT_ACCURACY = 1e-8; + + //protected int _maxNumIters; + //protected double _xmin = double.MinValue; + //protected double _xmax = double.MaxValue; + //protected double _accuracy; + //protected Func _func; + //protected Func m_Of; + int _maxNumIters; + double _xmin = double.MinValue; + double _xmax = double.MaxValue; + double _accuracy; + Func _func; + private double bracketingFactor = 1.6; + + + + /// Constructor. + /// A continuous function. + public RootFinder() : this(DEFAULT_MAX_ITERATIONS, DEFAULT_ACCURACY) + { + } + + public RootFinder(int numIters, double accuracy) + { + _maxNumIters = numIters; + _accuracy = accuracy; + } + + #region Properties + protected double XMin { get { return _xmin; } } + protected double XMax { get { return _xmax; } } + + public Func Func + { + get { return _func; } + set { _func = value; } + } + public double BracketingFactor + { + get { return bracketingFactor; } + set + { + if (value <= 0.0) throw new ArgumentOutOfRangeException(); + bracketingFactor = value; + } + } + public int Iterations + { + set + { + if (value <= 0) throw new ArgumentOutOfRangeException(); + _maxNumIters = value; + } + protected get { return _maxNumIters; } + } + public double Accuracy + { + get { return _accuracy; } + set { _accuracy = value; } + } + #endregion Properties + + /// Detect a range containing at least one root. + /// Lower value of the range. + /// Upper value of the range + /// The growing factor of research. Usually 1.6. + /// True if the bracketing operation succeeded, else otherwise. + /// This iterative methods stops when two values with opposite signs are found. + public bool SearchBracketsOutward(ref double xmin, ref double xmax, double factor) + { + if (xmin >= xmax) + { + throw new RootFinderException(INVALID_RANGE, 0, new Range(xmin, xmax), 0.0); + } + + double fmin = _func(xmin); + double fmax = _func(xmax); + + int i = 0; + while (i++ < _maxNumIters) + { + if (Math.Sign(fmin) != Math.Sign(fmax)) return true; + if (Math.Abs(fmin) < Math.Abs(fmax)) + { + xmin += factor * (xmin - xmax); + fmin = _func(xmin); + } + else + { + xmax += factor * (xmax - xmin); + fmax = _func(xmax); + } + } + + throw new RootFinderException(ROOT_NOT_FOUND, i, new Range(fmin, fmax), 0.0); + } + + /// 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(); + + /// Helper method useful for preventing rounding errors. + /// a*sign(b) + protected static double Sign(double a, double b) + { + return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); + } + + protected static bool Close(double d1, double d2) + { + return Math.Abs(d1 - d2) <= double.Epsilon; + } + } +} diff --git a/src/Numerics/RootFinders/RootFinderTest.cs b/src/Numerics/RootFinders/RootFinderTest.cs new file mode 100644 index 00000000..02f731de --- /dev/null +++ b/src/Numerics/RootFinders/RootFinderTest.cs @@ -0,0 +1,24 @@ +namespace MathNet.Numerics.RootFinders +{ + using System; + using System.Collections.Generic; + using MbUnit.Framework; + using Gallio.Framework; + using SuntrustPortfolio.Numerics; + + [TestFixture] + public class RootFinderTest + { + BrentRootFinder _solver = new BrentRootFinder(100, 1e-14); + + [Test] + public void MultipleRoots() + { + Func f = (x) => { return x * x - 4; }; + _solver.Func = f; + double root = _solver.Solve(-5, 5); + + Assert.AreEqual(0, f(root)); + } + } +}