diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index c3d61851..5b16823e 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -113,6 +113,7 @@ + diff --git a/src/Numerics/Optimization/BisectionRootFinder.cs b/src/Numerics/Optimization/BisectionRootFinder.cs new file mode 100644 index 00000000..f19444f1 --- /dev/null +++ b/src/Numerics/Optimization/BisectionRootFinder.cs @@ -0,0 +1,78 @@ +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 UpperExapansionFactor { 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.UpperExapansionFactor = upper_expansion_factor; + this.MaxExpansionSteps = max_expansion_steps; + } + + public double FindRoot(Func 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.UpperExapansionFactor <= 1.0) + throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case."); + + 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)); + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs b/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs new file mode 100644 index 00000000..561e9965 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs @@ -0,0 +1,30 @@ +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((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((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)); + } + + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 7c771f55..83be2101 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -751,6 +751,7 @@ + @@ -1087,6 +1088,7 @@ +