diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 941675c9..0d7d3fad 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -210,6 +210,7 @@ + diff --git a/src/Numerics/RootFinding/Secant.cs b/src/Numerics/RootFinding/Secant.cs new file mode 100644 index 00000000..69e755aa --- /dev/null +++ b/src/Numerics/RootFinding/Secant.cs @@ -0,0 +1,108 @@ +// +// 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.Properties; + +namespace MathNet.Numerics.RootFinding +{ + /// + /// Pure Secant root-finding algorithm without any recovery measures in cases it behaves badly. + /// The algorithm aborts immediately if the root leaves the bound interval. + /// + /// + public static class Secant + { + /// Find a solution of the equation f(x)=0. + /// The function to find roots from. + /// The first guess of the root within the bounds specified. + /// The second guess of the root within the bounds specified. + /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval. Default MinValue. + /// The high value of the range where the root is supposed to be. Aborts if it leaves the interval. Default MaxValue. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Default 1e-8. + /// Maximum number of iterations. Default 100. + /// Returns the root with the specified accuracy. + /// + public static double FindRoot(Func f, double guess, double secondGuess, double lowerBound = double.MinValue, double upperBound = double.MaxValue, double accuracy = 1e-8, int maxIterations = 100) + { + double root; + if (TryFindRoot(f, guess, secondGuess, lowerBound, upperBound, accuracy, maxIterations, out root)) + { + return root; + } + + throw new NonConvergenceException(Resources.RootFindingFailed); + } + + /// Find a solution of the equation f(x)=0. + /// The function to find roots from. + /// The first guess of the root within the bounds specified. + /// The second guess of the root within the bounds specified. + /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval. + /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval. + /// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached. Example: 1e-14. + /// Maximum number of iterations. Example: 100. + /// The root that was found, if any. Undefined if the function returns false. + /// True if a root with the specified accuracy was found, else false + public static bool TryFindRoot(Func f, double guess, double secondGuess, double lowerBound, double upperBound, double accuracy, int maxIterations, out double root) + { + root = secondGuess; + + // Either guess is outside of bounds + if (guess <= lowerBound || guess >= upperBound || secondGuess <= lowerBound || secondGuess >= upperBound) + { + return false; + } + + // Evaluation + double fguess = f(guess); + double froot = f(root); + + for (int i = 0; i <= maxIterations && root >= lowerBound && root <= upperBound; i++) + { + // Secant step + double step = froot * (root - guess) / (froot - fguess); + + guess = root; + fguess = froot; + + root -= step; + froot = f(root); + + if (Math.Abs(step) < accuracy && Math.Abs(froot) < accuracy) + { + return true; + } + } + + return false; + } + } +} diff --git a/src/UnitTests/RootFindingTests/SecantTest.cs b/src/UnitTests/RootFindingTests/SecantTest.cs new file mode 100644 index 00000000..ddce1efd --- /dev/null +++ b/src/UnitTests/RootFindingTests/SecantTest.cs @@ -0,0 +1,90 @@ +// +// 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, Category("RootFinding")] + public class SecantTest + { + [Test] + public void MultipleRoots() + { + // Roots at -2, 2 + Func f1 = x => x * x - 4; + Assert.AreEqual(0, f1(Secant.FindRoot(f1, 0, 6, -10, 10, 1e-14))); + Assert.AreEqual(-2, Secant.FindRoot(f1, -3, 0, -10, 10, 1e-14)); + Assert.AreEqual(2, Secant.FindRoot(f1, 1, 4, -10, 10, 1e-14)); + Assert.AreEqual(0, f1(Secant.FindRoot(x => f1(x), -6, 0, -10, 10, 1e-14))); + Assert.AreEqual(-2, Secant.FindRoot(x => f1(x), -5, -1, -10, 10, 1e-14)); + Assert.AreEqual(2, Secant.FindRoot(x => f1(x), 1, 4, -10, 10, 1e-14)); + + // Roots at 3, 4 + Func f2 = x => (x - 3) * (x - 4); + Assert.AreEqual(0, f2(Secant.FindRoot(f2, 3.5, 6, 0, 10, 1e-14))); + Assert.AreEqual(3, Secant.FindRoot(f2, -5, 3.5, -10, 10, 1e-14)); + Assert.AreEqual(4, Secant.FindRoot(f2, 3.2, 4.2, 0, 10, 1e-14)); + Assert.AreEqual(3, Secant.FindRoot(f2, 2.1, 3.1, 0, 10, 0.001, 50), 0.001); + Assert.AreEqual(3, Secant.FindRoot(f2, 2.1, 3.4, 0, 10, 0.001, 50), 0.001); + } + + [Test] + public void LocalMinima() + { + Func f1 = x => x * x * x - 2 * x + 2; + Assert.AreEqual(0, f1(Secant.FindRoot(f1, -5, 6, -10, 10, 1e-14))); + } + + [Test] + public void Cubic() + { + // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6 + Func f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3); + Assert.AreEqual(-1.265328088928, Secant.FindRoot(f1, -2, -1, -10, 10, 1e-10), 1e-6); + Assert.AreEqual(-1.265328088928, Secant.FindRoot(f1, -5, 6, -10, 10, 1e-10), 1e-6); + + // real roots only: 2x^3 + 4x^2 - 50x + 6 + Func f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2); + Assert.AreEqual(-6.1466562197069, Secant.FindRoot(f2, -8, -5, -10, 10, 1e-10), 1e-6); + Assert.AreEqual(0.12124737195841, Secant.FindRoot(f2, -1, 1, -10, 10, 1e-10), 1e-6); + Assert.AreEqual(4.0254088477485, Secant.FindRoot(f2, 3, 5, 0, 10, 1e-10), 1e-6); + } + + [Test] + public void NoRoot() + { + Func f1 = x => x * x + 4; + Assert.That(() => Secant.FindRoot(f1, -5, 5, -10, 10, 1e-14, 50), Throws.TypeOf()); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 713d8740..6552488c 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -367,6 +367,7 @@ +