uo9xpb74
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

added BrentRootFinder

v2
Candy Chiu 14 years ago
committed by Christoph Ruegg
parent
bae2c9614f
commit
51e529e07e
3 changed files with 313 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 114
      src/Numerics/RootFinders/BrentRootFinder.cs
  2. 175
      src/Numerics/RootFinders/RootFinder.cs
  3. 24
      src/Numerics/RootFinders/RootFinderTest.cs

114
src/Numerics/RootFinders/BrentRootFinder.cs
View File

@ -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));
}
}
}

175
src/Numerics/RootFinders/RootFinder.cs
View File

@ -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<double, double> _func;
//protected Func<double, double> m_Of;
int _maxNumIters;
double _xmin = double.MinValue;
double _xmax = double.MaxValue;
double _accuracy;
Func<double, double> _func;
private double bracketingFactor = 1.6;
/// <summary>Constructor.</summary>
/// <param name="f">A continuous function.</param>
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<double, double> 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
/// <summary>Detect a range containing at least one root.</summary>
/// <param name="xmin">Lower value of the range.</param>
/// <param name="xmax">Upper value of the range</param>
/// <param name="factor">The growing factor of research. Usually 1.6.</param>
/// <returns>True if the bracketing operation succeeded, else otherwise.</returns>
/// <remarks>This iterative methods stops when two values with opposite signs are found.</remarks>
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);
}
/// <summary>Prototype algorithm for solving the equation f(x)=0.</summary>
/// <param name="x1">The low value of the range where the root is supposed to be.</param>
/// <param name="x2">The high value of the range where the root is supposed to be.</param>
/// <returns>Returns the root with the specified accuracy.</returns>
public virtual double Solve(double x1, double x2)
{
_xmin = x1;
_xmax = x2;
return Find();
}
protected abstract double Find();
/// <summary>Helper method useful for preventing rounding errors.</summary>
/// <returns>a*sign(b)</returns>
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;
}
}
}

24
src/Numerics/RootFinders/RootFinderTest.cs
View File

@ -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<double, double> f = (x) => { return x * x - 4; };
_solver.Func = f;
double root = _solver.Solve(-5, 5);
Assert.AreEqual(0, f(root));
}
}
}
Write Preview
Loading…
Cancel
Save