diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 9249bd1a..2eee6c0b 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -109,6 +109,9 @@
+
+
+
diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs
index 13e9f0d2..0a336eb0 100644
--- a/src/Numerics/Properties/Resources.Designer.cs
+++ b/src/Numerics/Properties/Resources.Designer.cs
@@ -60,6 +60,15 @@ namespace MathNet.Numerics.Properties {
}
}
+ ///
+ /// Looks up a localized string similar to The accuracy couldn't be reached with the specified number of iterations..
+ ///
+ public static string AccuracyNotReached {
+ get {
+ return ResourceManager.GetString("AccuracyNotReached", resourceCulture);
+ }
+ }
+
///
/// Looks up a localized string similar to The array arguments must have the same length..
///
@@ -744,6 +753,15 @@ namespace MathNet.Numerics.Properties {
}
}
+ ///
+ /// Looks up a localized string similar to The algorithm ended without root in the range..
+ ///
+ public static string RootNotFound {
+ get {
+ return ResourceManager.GetString("RootNotFound", resourceCulture);
+ }
+ }
+
///
/// Looks up a localized string similar to The number of rows must greater than or equal to the number of columns..
///
diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx
index 2d7c5c10..e0645a91 100644
--- a/src/Numerics/Properties/Resources.resx
+++ b/src/Numerics/Properties/Resources.resx
@@ -369,10 +369,16 @@
We only support sparse matrix with less than int.MaxValue elements.
-
- Sample points should be sorted in strictly ascending order
+
+ Sample points should be sorted in strictly ascending order
-
- All sample points should be unique.
+
+ All sample points should be unique.
+
+
+ The accuracy couldn't be reached with the specified number of iterations.
+
+
+ The algorithm ended without root in the range.
\ No newline at end of file
diff --git a/src/Numerics/RootFinding/Bracketing.cs b/src/Numerics/RootFinding/Bracketing.cs
new file mode 100644
index 00000000..0cc1f603
--- /dev/null
+++ b/src/Numerics/RootFinding/Bracketing.cs
@@ -0,0 +1,48 @@
+using System;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.RootFinding
+{
+ public static class Bracketing
+ {
+ /// Detect a range containing at least one root.
+ /// The function to detect roots from.
+ /// 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)
+ {
+ if (xmin >= xmax)
+ {
+ throw new ArgumentOutOfRangeException("xmax", string.Format(Resources.ArgumentOutOfRangeGreater, "xmax", "xmin"));
+ }
+
+ double fmin = f(xmin);
+ double fmax = f(xmax);
+
+ for(int i=0;iFind 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)
+ {
+ double xMid = 0;
+ double d = 0.0, e = 0.0;
+
+ double fxmin = f(xmin);
+ double fxmax = f(xmax);
+ double root = xmax;
+ double froot = fxmax;
+
+ for (int i = 0; i <= maxIterations; i++)
+ {
+ // adjust bounds
+ if (Math.Sign(froot) == Math.Sign(fxmax))
+ {
+ 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
+ double xAcc1 = 2.0 * Precision.DoubleMachinePrecision * Math.Abs(root) + 0.5 * accuracy;
+ xMid = (xmax - root) / 2.0;
+ if (Math.Abs(xMid) <= xAcc1 || froot.AlmostEqualWithAbsoluteError(0, froot, accuracy))
+ {
+ return root;
+ }
+
+ if (Math.Abs(e) >= xAcc1 && Math.Abs(fxmin) > Math.Abs(froot))
+ {
+ // Attempt inverse quadratic interpolation
+ double s = froot / fxmin;
+ double p;
+ double q;
+ if (xmin.AlmostEqual(xmax))
+ {
+ p = 2.0 * xMid * s;
+ q = 1.0 - s;
+ }
+ else
+ {
+ q = fxmin / fxmax;
+ double 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)
+ {
+ // Check whether in bounds
+ q = -q;
+ }
+ p = Math.Abs(p);
+ if (2.0 * p < Math.Min(3.0 * xMid * q - Math.Abs(xAcc1 * q), Math.Abs(e * q)))
+ {
+ // Accept interpolation
+ e = d;
+ d = p / q;
+ }
+ else
+ {
+ // Interpolation failed, use bisection
+ d = xMid;
+ 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 = f(root);
+ }
+
+ // The algorithm has exceeded the number of iterations allowed
+ throw new RootFindingException(Resources.AccuracyNotReached, maxIterations, xmin, xmax, Math.Abs(xMid));
+ }
+
+ /// Helper method useful for preventing rounding errors.
+ /// a*sign(b)
+ static double Sign(double a, double b)
+ {
+ return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
+ }
+ }
+}
diff --git a/src/Numerics/RootFinding/RootFindingException.cs b/src/Numerics/RootFinding/RootFindingException.cs
new file mode 100644
index 00000000..3aaac181
--- /dev/null
+++ b/src/Numerics/RootFinding/RootFindingException.cs
@@ -0,0 +1,21 @@
+using System;
+
+namespace MathNet.Numerics.RootFinding
+{
+ public class RootFindingException : Exception
+ {
+ public RootFindingException(string message, int iteration, double rangeMin, double rangeMax, double accuracy)
+ : base(message)
+ {
+ Iteration = iteration;
+ RangeMin = rangeMin;
+ RangeMax = rangeMax;
+ Accuracy = accuracy;
+ }
+
+ public int Iteration { get; set; }
+ public double RangeMin { get; set; }
+ public double RangeMax { get; set; }
+ public double Accuracy { set; get; }
+ }
+}
\ No newline at end of file
diff --git a/src/UnitTests/RootFindingTests/BrentTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs
new file mode 100644
index 00000000..2956caad
--- /dev/null
+++ b/src/UnitTests/RootFindingTests/BrentTest.cs
@@ -0,0 +1,16 @@
+using MathNet.Numerics.RootFinding;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.RootFindingTests
+{
+ [TestFixture]
+ public class BrentTest
+ {
+ [Test]
+ public void MultipleRoots()
+ {
+ 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 1307db01..7c771f55 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -770,6 +770,7 @@
+