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

RootFinding: move bisection to root-finding namespace, update portable project

v2
Christoph Ruegg 14 years ago
parent
729ba52bf3
commit
e523675a43
7 changed files with 186 additions and 133 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      src/Numerics/Numerics.csproj
  2. 101
      src/Numerics/Optimization/BisectionRootFinder.cs
  3. 126
      src/Numerics/RootFinding/Algorithms/Bisection.cs
  4. 3
      src/Portable/Portable.csproj
  5. 30
      src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs
  6. 55
      src/UnitTests/RootFindingTests/BisectionTest.cs
  7. 2
      src/UnitTests/UnitTests.csproj

2
src/Numerics/Numerics.csproj
View File

@ -113,7 +113,7 @@
<Compile Include="RootFinding\Bracketing.cs" />
<Compile Include="RootFinding\Algorithms\Brent.cs" />
<Compile Include="RootFinding\FloatingPointRoots.cs" />
<Compile Include="Optimization\BisectionRootFinder.cs" />
<Compile Include="RootFinding\Algorithms\Bisection.cs" />
<Compile Include="SpecialFunctions\Evaluate.cs" />
<Compile Include="SpecialFunctions\ModifiedStruve.cs" />
<Compile Include="SpecialFunctions\ModifiedBessel.cs" />

101
src/Numerics/Optimization/BisectionRootFinder.cs
View File

@ -1,101 +0,0 @@
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace MathNet.Numerics.Optimization
{
public class BisectionRootFinder
{
public double ObjectiveTolerance { get; set; }
public double XTolerance { get; set; }
public double LowerExpansionFactor { get; set; }
public double UpperExpansionFactor { get; set; }
public int MaxExpansionSteps { get; set; }
public BisectionRootFinder(double objective_tolerance=1e-5, double x_tolerance=1e-5, double lower_expansion_factor=-1.0, double upper_expansion_factor=-1.0, int max_expansion_steps=10)
{
this.ObjectiveTolerance = objective_tolerance;
this.XTolerance = x_tolerance;
this.LowerExpansionFactor = lower_expansion_factor;
this.UpperExpansionFactor = upper_expansion_factor;
this.MaxExpansionSteps = max_expansion_steps;
}
public double FindRoot(Func<double, double> objective_function, double lower_bound, double upper_bound)
{
double lower_val = objective_function(lower_bound);
double upper_val = objective_function(upper_bound);
if (lower_val == 0.0)
return lower_bound;
if (upper_val == 0.0)
return upper_bound;
this.ValidateEvaluation(lower_val, lower_bound);
this.ValidateEvaluation(upper_val, upper_bound);
if (Math.Sign(lower_val) == Math.Sign(upper_val) && this.LowerExpansionFactor <= 1.0 && this.UpperExpansionFactor <= 1.0)
throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case.");
int expansion_steps = 0;
while (Math.Sign(lower_val) == Math.Sign(upper_val) && expansion_steps < this.MaxExpansionSteps)
{
double midpoint = 0.5 * (upper_bound + lower_bound);
double range = upper_bound - lower_bound;
if (this.UpperExpansionFactor <= 0.0 || (this.LowerExpansionFactor > 0.0 && Math.Abs(lower_val) < Math.Abs(upper_val)) )
{
lower_bound = upper_bound - this.LowerExpansionFactor * range;
lower_val = objective_function(lower_bound);
this.ValidateEvaluation(lower_val, lower_bound);
}
else
{
upper_bound = lower_bound + this.UpperExpansionFactor * range;
upper_val = objective_function(upper_bound);
this.ValidateEvaluation(upper_val, upper_bound);
}
expansion_steps += 1;
}
if (expansion_steps == this.MaxExpansionSteps)
throw new MaximumIterationsException("Could not bound root in maximum expansion iterations.");
while (Math.Abs(upper_val - lower_val) > 0.5 * this.ObjectiveTolerance || Math.Abs(upper_bound - lower_bound) > 0.5 * this.XTolerance)
{
double midpoint = 0.5 * (upper_bound + lower_bound);
double midval = objective_function(midpoint);
this.ValidateEvaluation(midval, midpoint);
if (Math.Sign(midval) == Math.Sign(lower_val))
{
lower_bound = midpoint;
lower_val = midval;
}
else if (Math.Sign(midval) == Math.Sign(upper_val))
{
upper_bound = midpoint;
upper_val = midval;
}
else
{
return midpoint;
}
}
return 0.5 * (lower_bound + upper_bound);
}
private void ValidateEvaluation(double output, double input)
{
if (!IsFinite(output))
throw new Exception(String.Format("Objective function returned non-finite result: f({0}) = {1}", input, output));
}
private static bool IsFinite(double x)
{
return !(Double.IsInfinity(x) || Double.IsNaN(x));
}
}
}

126
src/Numerics/RootFinding/Algorithms/Bisection.cs
View File

@ -0,0 +1,126 @@
// <copyright file="Bisection.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
namespace MathNet.Numerics.RootFinding.Algorithms
{
public class Bisection
{
public Bisection(double objective_tolerance = 1e-5, double x_tolerance = 1e-5, double lower_expansion_factor = -1.0, double upper_expansion_factor = -1.0, int max_expansion_steps = 10)
{
ObjectiveTolerance = objective_tolerance;
XTolerance = x_tolerance;
LowerExpansionFactor = lower_expansion_factor;
UpperExpansionFactor = upper_expansion_factor;
MaxExpansionSteps = max_expansion_steps;
}
public double ObjectiveTolerance { get; set; }
public double XTolerance { get; set; }
public double LowerExpansionFactor { get; set; }
public double UpperExpansionFactor { get; set; }
public int MaxExpansionSteps { get; set; }
public double FindRoot(Func<double, double> objective_function, double lower_bound, double upper_bound)
{
double lower_val = objective_function(lower_bound);
double upper_val = objective_function(upper_bound);
if (lower_val == 0.0)
return lower_bound;
if (upper_val == 0.0)
return upper_bound;
ValidateEvaluation(lower_val, lower_bound);
ValidateEvaluation(upper_val, upper_bound);
if (Math.Sign(lower_val) == Math.Sign(upper_val) && LowerExpansionFactor <= 1.0 && UpperExpansionFactor <= 1.0)
throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case.");
int expansion_steps = 0;
while (Math.Sign(lower_val) == Math.Sign(upper_val) && expansion_steps < MaxExpansionSteps)
{
double range = upper_bound - lower_bound;
if (UpperExpansionFactor <= 0.0 || (LowerExpansionFactor > 0.0 && Math.Abs(lower_val) < Math.Abs(upper_val)))
{
lower_bound = upper_bound - LowerExpansionFactor*range;
lower_val = objective_function(lower_bound);
ValidateEvaluation(lower_val, lower_bound);
}
else
{
upper_bound = lower_bound + UpperExpansionFactor*range;
upper_val = objective_function(upper_bound);
ValidateEvaluation(upper_val, upper_bound);
}
expansion_steps += 1;
}
if (expansion_steps == MaxExpansionSteps)
throw new NonConvergenceException();
while (Math.Abs(upper_val - lower_val) > 0.5*ObjectiveTolerance || Math.Abs(upper_bound - lower_bound) > 0.5*XTolerance)
{
double midpoint = 0.5*(upper_bound + lower_bound);
double midval = objective_function(midpoint);
ValidateEvaluation(midval, midpoint);
if (Math.Sign(midval) == Math.Sign(lower_val))
{
lower_bound = midpoint;
lower_val = midval;
}
else if (Math.Sign(midval) == Math.Sign(upper_val))
{
upper_bound = midpoint;
upper_val = midval;
}
else
{
return midpoint;
}
}
return 0.5*(lower_bound + upper_bound);
}
void ValidateEvaluation(double output, double input)
{
if (!IsFinite(output))
throw new Exception(String.Format("Objective function returned non-finite result: f({0}) = {1}", input, output));
}
static bool IsFinite(double x)
{
return !(Double.IsInfinity(x) || Double.IsNaN(x));
}
}
}

3
src/Portable/Portable.csproj
View File

@ -984,6 +984,9 @@
<Compile Include="..\Numerics\Random\Xorshift.cs">
<Link>Random\Xorshift.cs</Link>
</Compile>
<Compile Include="..\Numerics\RootFinding\Algorithms\Bisection.cs">
<Link>RootFinding\Algorithms\Bisection.cs</Link>
</Compile>
<Compile Include="..\Numerics\RootFinding\Algorithms\Brent.cs">
<Link>RootFinding\Algorithms\Brent.cs</Link>
</Compile>

30
src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs
View File

@ -1,30 +0,0 @@
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using NUnit.Framework;
using MathNet.Numerics.Optimization;
namespace MathNet.Numerics.UnitTests.OptimizationTests
{
[TestFixture]
class TestBisectionRootFinder
{
[Test]
public void FindRoot_Works()
{
var algorithm = new BisectionRootFinder(0.001, 0.001);
var f1 = new Func<double, double>((x) => (x - 3) * (x - 4));
var r1 = algorithm.FindRoot(f1, 2.1, 3.9);
Assert.That(Math.Abs(f1(r1)), Is.LessThan(0.001));
Assert.That(Math.Abs(r1 - 3.0), Is.LessThan(0.001));
var f2 = new Func<double, double>((x) => (x - 3) * (x - 4));
var r2 = algorithm.FindRoot(f1, 2.1, 3.4);
Assert.That(Math.Abs(f2(r2)), Is.LessThan(0.001));
Assert.That(Math.Abs(r2 - 3.0), Is.LessThan(0.001));
}
}
}

55
src/UnitTests/RootFindingTests/BisectionTest.cs
View File

@ -0,0 +1,55 @@
// <copyright file="BisectionTest.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
using MathNet.Numerics.RootFinding.Algorithms;
using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.RootFindingTests
{
[TestFixture]
internal class BisectionTest
{
[Test]
public void FindRoot_Works()
{
var algorithm = new Bisection(0.001, 0.001);
var f1 = new Func<double, double>((x) => (x - 3)*(x - 4));
double r1 = algorithm.FindRoot(f1, 2.1, 3.9);
Assert.That(Math.Abs(f1(r1)), Is.LessThan(0.001));
Assert.That(Math.Abs(r1 - 3.0), Is.LessThan(0.001));
var f2 = new Func<double, double>((x) => (x - 3)*(x - 4));
double r2 = algorithm.FindRoot(f1, 2.1, 3.4);
Assert.That(Math.Abs(f2(r2)), Is.LessThan(0.001));
Assert.That(Math.Abs(r2 - 3.0), Is.LessThan(0.001));
}
}
}

2
src/UnitTests/UnitTests.csproj
View File

@ -751,7 +751,7 @@
<Compile Include="NumberTheoryTests\GcdRelatedTest.cs" />
<Compile Include="NumberTheoryTests\GcdRelatedTestBigInteger.cs" />
<Compile Include="NumberTheoryTests\IntegerTheoryTest.cs" />
<Compile Include="OptimizationTests\TestBisectionRootFinder.cs" />
<Compile Include="RootFindingTests\BisectionTest.cs" />
<Compile Include="PermutationTest.cs" />
<Compile Include="PrecisionTest.cs" />
<Compile Include="Properties\AssemblyInfo.cs" />

Write Preview
Loading…
Cancel
Save