Browse Source

RootFinding: extract bracketing to separate class

v2
Christoph Ruegg 14 years ago
parent
commit
2f9913921c
  1. 1
      src/Numerics/Numerics.csproj
  2. 48
      src/Numerics/RootFinding/Bracketing.cs
  3. 12
      src/Numerics/RootFinding/BrentRootFinder.cs
  4. 59
      src/Numerics/RootFinding/RootFinder.cs

1
src/Numerics/Numerics.csproj

@ -109,6 +109,7 @@
<Compile Include="Financial\AbsoluteRiskMeasures.cs" />
<Compile Include="LinearAlgebra\Generic\Matrix.BCL.cs" />
<Compile Include="LinearAlgebra\Generic\Vector.BCL.cs" />
<Compile Include="RootFinding\Bracketing.cs" />
<Compile Include="RootFinding\BrentRootFinder.cs" />
<Compile Include="RootFinding\RootFinder.cs" />
<Compile Include="RootFinding\RootFindingException.cs" />

48
src/Numerics/RootFinding/Bracketing.cs

@ -0,0 +1,48 @@
using System;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.RootFinding
{
public static class Bracketing
{
/// <summary>Detect a range containing at least one root.</summary>
/// <param name="f">The function to detect roots from.</param>
/// <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>
/// <param name="maxIterations">Maximum number of iterations. Usually 50.</param>
/// <returns>True if the bracketing operation succeeded, false otherwise.</returns>
/// <remarks>This iterative methods stops when two values with opposite signs are found.</remarks>
public static bool SearchOutward(Func<double, double> 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;i<maxIterations; i++)
{
if (Math.Sign(fmin) != Math.Sign(fmax))
{
return true;
}
if (Math.Abs(fmin) < Math.Abs(fmax))
{
xmin += factor * (xmin - xmax);
fmin = f(xmin);
}
else
{
xmax += factor * (xmax - xmin);
fmax = f(xmax);
}
}
return false;
}
}
}

12
src/Numerics/RootFinding/BrentRootFinder.cs

@ -111,5 +111,17 @@ namespace MathNet.Numerics.RootFinding
// The algorithm has exceeded the number of iterations allowed
throw new RootFindingException(Resources.AccuracyNotReached, i, XMin, XMax, Math.Abs(xMid));
}
/// <summary>Helper method useful for preventing rounding errors.</summary>
/// <returns>a*sign(b)</returns>
static double Sign(double a, double b)
{
return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
}
static bool Close(double d1, double d2)
{
return Math.Abs(d1 - d2) <= double.Epsilon;
}
}
}

59
src/Numerics/RootFinding/RootFinder.cs

@ -1,5 +1,4 @@
using System;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.RootFinding
{
@ -13,7 +12,6 @@ namespace MathNet.Numerics.RootFinding
double _xmin = double.MinValue;
double _xmax = double.MaxValue;
Func<double, double> _func;
private double _bracketingFactor = 1.6;
public RootFinder() : this(DefaultMaxIterations, DefaultAccuracy)
{
@ -36,16 +34,6 @@ namespace MathNet.Numerics.RootFinding
public double Accuracy { get; set; }
public double BracketingFactor
{
get { return _bracketingFactor; }
set
{
if (value <= 0.0) throw new ArgumentOutOfRangeException();
_bracketingFactor = value;
}
}
public int Iterations
{
set
@ -56,41 +44,6 @@ namespace MathNet.Numerics.RootFinding
protected get { return _maxNumIters; }
}
/// <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 RootFindingException(string.Format(Resources.ArgumentOutOfRangeGreater,"xmax","xmin"), 0, 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 RootFindingException(Resources.RootNotFound, i, 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>
@ -103,17 +56,5 @@ namespace MathNet.Numerics.RootFinding
}
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;
}
}
}

Loading…
Cancel
Save