From a12684649eff18b834fb67f49e1e713458a221dc Mon Sep 17 00:00:00 2001 From: taschna Date: Wed, 17 Jul 2013 09:31:37 +0200 Subject: [PATCH] RootFinding: Add fallback to Bisection in FindRoots.OfFunction. Add fallback to Bisection in FindRoots.OfFunction. Add test cases for FindRoots.OfFunction. --- src/Numerics/FindRoots.cs | 13 +- src/UnitTests/RootFindingTests/BrentTest.cs | 2 + .../RootFindingTests/FindRootsTest.cs | 158 ++++++++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 4 files changed, 173 insertions(+), 1 deletion(-) create mode 100644 src/UnitTests/RootFindingTests/FindRootsTest.cs diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs index d6387547..deef6948 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -37,7 +37,18 @@ namespace MathNet.Numerics { public static double OfFunction(Func f, double lowerBound, double upperBound, double accuracy = 1e-8) { - return Brent.FindRoot(f, lowerBound, upperBound, accuracy, 100); + double root; + if (Brent.TryFindRoot(f, lowerBound, upperBound, accuracy, 100, out root)) + { + return root; + } + + if (Bisection.TryFindRoot(f, lowerBound, upperBound, accuracy, 100, out root)) + { + return root; + } + + throw new NonConvergenceException("The algorithm has failed or exceeded the number of iterations allowed"); } public static double OfFunctionAndDerivative(Func f, Func df, double lowerBound, double upperBound, double accuracy = 1e-8) diff --git a/src/UnitTests/RootFindingTests/BrentTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs index 451840a6..b57df784 100644 --- a/src/UnitTests/RootFindingTests/BrentTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -90,6 +90,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void Oneeq1() { + // Test case from http://www.polymath-software.com/library/nle/Oneeq1.htm Func f1 = z => 8 * Math.Pow((4 - z) * z, 2) / (Math.Pow(6 - 3 * z, 2) * (2 - z)) - 0.186; double x = Brent.FindRoot(f1, 0.1, 0.9, 1e-14, 100); Assert.AreEqual(0.277759543089215, x, 1e-9); @@ -99,6 +100,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests [Test] public void Oneeq2a() { + // Test case from http://www.polymath-software.com/library/nle/Oneeq2a.htm Func f1 = T => { const double x1 = 0.332112; diff --git a/src/UnitTests/RootFindingTests/FindRootsTest.cs b/src/UnitTests/RootFindingTests/FindRootsTest.cs new file mode 100644 index 00000000..60ea98f2 --- /dev/null +++ b/src/UnitTests/RootFindingTests/FindRootsTest.cs @@ -0,0 +1,158 @@ +// +// 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. +// + +using System; +using MathNet.Numerics.RootFinding; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.RootFindingTests +{ + [TestFixture] + internal class FindRootsTest + { + [Test] + public void MultipleRoots() + { + // Roots at -2, 2 + Func f1 = x => x * x - 4; + Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, 0, 5, 1e-14))); + Assert.AreEqual(-2, FindRoots.OfFunction(f1, -5, -1, 1e-14)); + Assert.AreEqual(2, FindRoots.OfFunction(f1, 1, 4, 1e-14)); + Assert.AreEqual(0, f1(FindRoots.OfFunction(x => -f1(x), 0, 5, 1e-14))); + Assert.AreEqual(-2, FindRoots.OfFunction(x => -f1(x), -5, -1, 1e-14)); + Assert.AreEqual(2, FindRoots.OfFunction(x => -f1(x), 1, 4, 1e-14)); + + // Roots at 3, 4 + Func f2 = x => (x - 3) * (x - 4); + Assert.AreEqual(0, f2(FindRoots.OfFunction(f2, 3.5, 5, 1e-14)), 1e-14); + Assert.AreEqual(3, FindRoots.OfFunction(f2, -5, 3.5, 1e-14)); + Assert.AreEqual(4, FindRoots.OfFunction(f2, 3.2, 5, 1e-14)); + Assert.AreEqual(3, FindRoots.OfFunction(f2, 2.1, 3.9, 0.001), 0.001); + Assert.AreEqual(3, FindRoots.OfFunction(f2, 2.1, 3.4, 0.001), 0.001); + } + + [Test] + public void LocalMinima() + { + Func f1 = x => x * x * x - 2 * x + 2; + Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, -5, 5, 1e-14)), 1e-14); + Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, -2, 4, 1e-14)), 1e-14); + } + + [Test] + public void NoRoot() + { + Func f1 = x => x * x + 4; + Assert.Throws(() => FindRoots.OfFunction(f1, -5, 5, 1e-14)); + } + + [Test] + public void Oneeq3() + { + // Test case from http://www.polymath-software.com/library/nle/Oneeq3.htm + Func f1 = T => Math.Exp(21000 / T) / (T * T) - 1.11e11; + double x = FindRoots.OfFunction(f1, 550, 560, 1e-2); + Assert.AreEqual(551.773822885233, x, 1e-5); + Assert.AreEqual(0, f1(x), 1e-2); + } + + [Test] + public void Oneeq5() + { + // Test case from http://www.polymath-software.com/library/nle/Oneeq5.htm + Func f1 = TR => + { + const double ALPHA = 7.256; + const double BETA = 2.298E-3; + const double GAMMA = 0.283E-6; + const double DH = -57798.0; + return DH + TR * (ALPHA + TR * (BETA / 2 + TR * GAMMA / 3)) - 298.0 * (ALPHA + 298.0 * (BETA / 2 + 298.0 * GAMMA / 3)); + }; + + double x = FindRoots.OfFunction(f1, 3000, 5000, 1e-10); + Assert.AreEqual(4305.30991366774, x, 1e-5); + Assert.AreEqual(0, f1(x), 1e-10); + } + + [Test] + public void Oneeq6a() + { + // Test case from http://www.polymath-software.com/library/nle/Oneeq6a.htm + Func f1 = V => + { + const double R = 0.08205; + const double T = 273.0; + const double B0 = 0.05587; + const double A0 = 2.2769; + const double C = 128300.0; + const double A = 0.01855; + const double B = -0.01587; + const double P = 100.0; + const double Beta = R * T * B0 - A0 - R * C / (T * T); + const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); + const double Delta = R * B0 * B * C / (T * T); + + return R * T / V + Beta / (V * V) + Gama / (V * V * V) + Delta / (V * V * V * V) - P; + }; + + double x = FindRoots.OfFunction(f1, 0.1, 1); + Assert.AreEqual(0.174749531708621, x, 1e-5); + Assert.AreEqual(0, f1(x), 1e-13); + } + + [Test] + public void Oneeq7() + { + // Test case from http://www.polymath-software.com/library/nle/Oneeq7.htm + Func f1 = x => x / (1 - x) - 5 * Math.Log(0.4 * (1 - x) / (0.4 - 0.5 * x)) + 4.45977; + double r = FindRoots.OfFunction(f1, 0, 0.79, 1e-2); + Assert.AreEqual(0.757396293891, r, 1e-5); + Assert.AreEqual(0, f1(r), 1e-13); + } + + [Test] + public void Oneeq8() + { + // Test case from http://www.polymath-software.com/library/nle/Oneeq8.htm + Func f1 = v => + { + const double a = 240; + const double b = 40; + const double c = 200; + + return a * v * v + b * Math.Pow(v, 7 / 4) - c; + }; + + double x = FindRoots.OfFunction(f1, 0.01, 1); + Assert.AreEqual(0.842524411168525, x, 1e-2); + Assert.AreEqual(0, f1(x), 1e-13); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 24e0b96b..b2adc9ac 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -773,6 +773,7 @@ +