From 1b347b415536dba05ebe295621ba897d54ee5e98 Mon Sep 17 00:00:00 2001 From: taschna Date: Thu, 11 Jul 2013 12:51:15 +0200 Subject: [PATCH] Add test for Brent and Bisection. Fix termination condition of Bisection.FindRoot. Add test case Oneeq2a (http://www.polymath-software.com/library/nle/Oneeq2a.htm) for Brent and Bisection. Fix termination condition of Bisection.FindRoot (new test case results in infinite loop otherwise). --- src/Numerics/RootFinding/Bisection.cs | 5 ++++ .../RootFindingTests/BisectionTest.cs | 25 +++++++++++++++++++ src/UnitTests/RootFindingTests/BrentTest.cs | 25 +++++++++++++++++++ 3 files changed, 55 insertions(+) diff --git a/src/Numerics/RootFinding/Bisection.cs b/src/Numerics/RootFinding/Bisection.cs index 739b26e6..c0b907cc 100644 --- a/src/Numerics/RootFinding/Bisection.cs +++ b/src/Numerics/RootFinding/Bisection.cs @@ -67,6 +67,11 @@ namespace MathNet.Numerics.RootFinding while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision) { double midpoint = 0.5*(upperBound + lowerBound); + if ((midpoint == lowerBound) || (midpoint == upperBound)) + { + return midpoint; + } + double midval = f(midpoint); if (Math.Sign(midval) == Math.Sign(fmin)) diff --git a/src/UnitTests/RootFindingTests/BisectionTest.cs b/src/UnitTests/RootFindingTests/BisectionTest.cs index 8e038ad8..aced22cc 100644 --- a/src/UnitTests/RootFindingTests/BisectionTest.cs +++ b/src/UnitTests/RootFindingTests/BisectionTest.cs @@ -82,5 +82,30 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.AreEqual(0.277759543089215, x, 1e-9); Assert.AreEqual(0, f1(x), 1e-16); } + + [Test] + public void Oneeq2a() + { + Func f1 = T => + { + const double x1 = 0.332112; + const double x2 = 1 - x1; + const double G12 = 0.07858889; + const double G21 = 0.30175355; + double P2 = Math.Pow(10, 6.87776 - 1171.53 / (224.366 + T)); + double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + T)); + double t1 = x1 + x2 * G12; + double t2 = x2 + x1 * G21; + double gamma2 = Math.Exp(-Math.Log(t2) - (x1 * (G12 * t2 - G21 * t1)) / (t1 * t2)); + double gamma1 = Math.Exp(-Math.Log(t1) + (x2 * (G12 * t2 - G21 * t1)) / (t1 * t2)); + double k1 = gamma1 * P1 / 760; + double k2 = gamma2 * P2 / 760; + return 1 - k1 * x1 - k2 * x2; + }; + + double x = Bisection.FindRoot(f1, 0, 100); + Assert.AreEqual(58.7000023925600, x, 1e-5); + Assert.AreEqual(0, f1(x), 1e-14); + } } } diff --git a/src/UnitTests/RootFindingTests/BrentTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs index b4e8a70f..9c06c2f4 100644 --- a/src/UnitTests/RootFindingTests/BrentTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -95,5 +95,30 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.AreEqual(0.277759543089215, x, 1e-9); Assert.AreEqual(0, f1(x), 1e-16); } + + [Test] + public void Oneeq2a() + { + Func f1 = T => + { + const double x1 = 0.332112; + const double x2 = 1 - x1; + const double G12 = 0.07858889; + const double G21 = 0.30175355; + double P2 = Math.Pow(10, 6.87776 - 1171.53 / (224.366 + T)); + double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + T)); + double t1 = x1 + x2 * G12; + double t2 = x2 + x1 * G21; + double gamma2 = Math.Exp(-Math.Log(t2) - (x1 * (G12 * t2 - G21 * t1)) / (t1 * t2)); + double gamma1 = Math.Exp(-Math.Log(t1) + (x2 * (G12 * t2 - G21 * t1)) / (t1 * t2)); + double k1 = gamma1 * P1 / 760; + double k2 = gamma2 * P2 / 760; + return 1 - k1 * x1 - k2 * x2; + }; + + double x = Brent.FindRoot(f1, 0, 100, 1e-14, 100); + Assert.AreEqual(58.7000023925600, x, 1e-5); + Assert.AreEqual(0, f1(x), 1e-14); + } } }