Browse Source

RootFinding: refactoring, simplify

pull/121/head
Christoph Ruegg 13 years ago
parent
commit
04a5f30d99
  1. 3
      src/Numerics/Numerics.csproj
  2. 45
      src/Numerics/RootFinding/FindRoots.cs
  3. 60
      src/Numerics/RootFinding/RootFinder.cs
  4. 7
      src/UnitTests/RootFindingTests/BrentTest.cs
  5. 2
      src/UnitTests/UnitTests.csproj

3
src/Numerics/Numerics.csproj

@ -110,8 +110,7 @@
<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\FindRoots.cs" />
<Compile Include="RootFinding\RootFindingException.cs" />
<Compile Include="SpecialFunctions\Evaluate.cs" />
<Compile Include="SpecialFunctions\ModifiedStruve.cs" />

45
src/Numerics/RootFinding/BrentRootFinder.cs → src/Numerics/RootFinding/FindRoots.cs

@ -3,39 +3,32 @@ using MathNet.Numerics.Properties;
namespace MathNet.Numerics.RootFinding
{
public class BrentRootFinder : RootFinder
public static class FindRoots
{
public BrentRootFinder() : base()
/// <summary>Find a solution of the equation f(x)=0.</summary>
/// <param name="f">The function to find roots from.</param>
/// <param name="xmin">The low value of the range where the root is supposed to be.</param>
/// <param name="xmax">The high value of the range where the root is supposed to be.</param>
/// <param name="accuracy">Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.</param>
/// <param name="maxIterations">Maximum number of iterations. Usually 100.</param>
/// <returns>Returns the root with the specified accuracy.</returns>
/// <remarks>
/// 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
/// </remarks>
public static double BrentMethod(Func<double, double> 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 @@ namespace MathNet.Numerics.RootFinding
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 @@ namespace MathNet.Numerics.RootFinding
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));
}
/// <summary>Helper method useful for preventing rounding errors.</summary>

60
src/Numerics/RootFinding/RootFinder.cs

@ -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<double, double> _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<double, double> 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; }
}
/// <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();
}
}

7
src/UnitTests/RootFindingTests/BrentRootFindingTest.cs → src/UnitTests/RootFindingTests/BrentTest.cs

@ -4,14 +4,13 @@ using NUnit.Framework;
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);
}
}
}

2
src/UnitTests/UnitTests.csproj

@ -770,7 +770,7 @@
<Compile Include="Random\WH1982Tests.cs" />
<Compile Include="Random\WH2006Tests.cs" />
<Compile Include="Random\XorshiftTests.cs" />
<Compile Include="RootFindingTests\BrentRootFindingTest.cs" />
<Compile Include="RootFindingTests\BrentTest.cs" />
<Compile Include="SortingTests.cs" />
<Compile Include="SpecialFunctionsTests\ModifiedStruveTests.cs" />
<Compile Include="SpecialFunctionsTests\ModifiedBesselTests.cs" />

Loading…
Cancel
Save