diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b839b90c2..6d25cd659 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -110,6 +110,7 @@ + diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 0a336eb01..2a09d9bb3 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -753,6 +753,15 @@ public static string ProposalDistributionNoUpperBound { } } + /// + /// Looks up a localized string similar to The lower and upper bounds must bracket a single root.. + /// + public static string RootMustBeBracketedByBounds { + get { + return ResourceManager.GetString("RootMustBeBracketedByBounds", resourceCulture); + } + } + /// /// Looks up a localized string similar to The algorithm ended without root in the range.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index e0645a91d..a402e54b9 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -381,4 +381,7 @@ The algorithm ended without root in the range. + + The lower and upper bounds must bracket a single root. + \ No newline at end of file diff --git a/src/Numerics/Properties/Resources1.Designer.cs b/src/Numerics/Properties/Resources1.Designer.cs index 13e9f0d28..2a09d9bb3 100644 --- a/src/Numerics/Properties/Resources1.Designer.cs +++ b/src/Numerics/Properties/Resources1.Designer.cs @@ -60,6 +60,15 @@ internal Resources() { } } + /// + /// Looks up a localized string similar to The accuracy couldn't be reached with the specified number of iterations.. + /// + public static string AccuracyNotReached { + get { + return ResourceManager.GetString("AccuracyNotReached", resourceCulture); + } + } + /// /// Looks up a localized string similar to The array arguments must have the same length.. /// @@ -744,6 +753,24 @@ public static string ProposalDistributionNoUpperBound { } } + /// + /// Looks up a localized string similar to The lower and upper bounds must bracket a single root.. + /// + public static string RootMustBeBracketedByBounds { + get { + return ResourceManager.GetString("RootMustBeBracketedByBounds", resourceCulture); + } + } + + /// + /// Looks up a localized string similar to The algorithm ended without root in the range.. + /// + public static string RootNotFound { + get { + return ResourceManager.GetString("RootNotFound", resourceCulture); + } + } + /// /// Looks up a localized string similar to The number of rows must greater than or equal to the number of columns.. /// diff --git a/src/Numerics/RootFinding/Algorithms/Bisection.cs b/src/Numerics/RootFinding/Algorithms/Bisection.cs index ca5ca3e40..09d1772f6 100644 --- a/src/Numerics/RootFinding/Algorithms/Bisection.cs +++ b/src/Numerics/RootFinding/Algorithms/Bisection.cs @@ -29,41 +29,39 @@ // using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.RootFinding.Algorithms { public static class Bisection { - public static double FindRootExpand(Func f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-5, double expandFactor = 1.6, int maxExpandIteratons = 100) + /// Find a solution of the equation f(x)=0. + /// + public static double FindRootExpand(Func f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, double expandFactor = 1.6, int maxExpandIteratons = 100) { Bracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons); return FindRoot(f, guessLowerBound, guessUpperBound, accuracy); } /// Find a solution of the equation f(x)=0. - public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy = 1e-5) + /// + public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy = 1e-8) { double fmin = f(lowerBound); double fmax = f(upperBound); - if (fmin == 0.0) - return lowerBound; - if (fmax == 0.0) - return upperBound; - - ValidateEvaluation(fmin, lowerBound); - ValidateEvaluation(fmax, upperBound); + if (fmin == 0.0) return lowerBound; + if (fmax == 0.0) return upperBound; if (Math.Sign(fmin) == Math.Sign(fmax)) { - throw new NonConvergenceException("Bounds do not necessarily span a root."); + throw new NonConvergenceException(Resources.RootMustBeBracketedByBounds); } while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision) { double midpoint = 0.5*(upperBound + lowerBound); double midval = f(midpoint); - ValidateEvaluation(midval, midpoint); if (Math.Sign(midval) == Math.Sign(fmin)) { @@ -83,13 +81,5 @@ public static double FindRoot(Func f, double lowerBound, double return 0.5*(lowerBound + upperBound); } - - static void ValidateEvaluation(double output, double input) - { - if (Double.IsInfinity(output) || Double.IsInfinity(output)) - { - throw new Exception(String.Format("Objective function returned non-finite result: f({0}) = {1}", input, output)); - } - } } } diff --git a/src/Numerics/RootFinding/Algorithms/Brent.cs b/src/Numerics/RootFinding/Algorithms/Brent.cs index df68e819f..75f789d8d 100644 --- a/src/Numerics/RootFinding/Algorithms/Brent.cs +++ b/src/Numerics/RootFinding/Algorithms/Brent.cs @@ -139,8 +139,7 @@ public static double FindRoot(Func f, double lowerBound, double froot = f(root); } - // The algorithm has exceeded the number of iterations allowed - throw new NonConvergenceException(); + throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed"); } /// Helper method useful for preventing rounding errors. diff --git a/src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs b/src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs new file mode 100644 index 000000000..f442ce27a --- /dev/null +++ b/src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs @@ -0,0 +1,107 @@ +// +// 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; + +namespace MathNet.Numerics.RootFinding.Algorithms +{ + public static class NewtonRaphson + { + /// Find a solution of the equation f(x)=0. + /// Hybrid Newton-Raphson that when failing falls back to bisection. + /// + public static double FindRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100) + { + double fmin = f(lowerBound); + double fmax = f(upperBound); + + if (fmin == 0.0) return lowerBound; + if (fmax == 0.0) return upperBound; + + double root = 0.5*(lowerBound + upperBound); + double fx = f(root); + double lastStep = Math.Abs(upperBound - lowerBound); + for (int i = 0; i < maxIterations; i++) + { + double dfx = df(root); + + // Netwon-Raphson step + double step = fx/dfx; + root -= step; + + if (root < lowerBound || root > upperBound || Math.Abs(2*fx) > Math.Abs(lastStep*dfx)) + { + // Newton-Raphson step failed -> bisect instead + root = 0.5*(upperBound + lowerBound); + fx = f(root); + lastStep = 0.5*Math.Abs(upperBound - lowerBound); + if (Math.Sign(fx) == Math.Sign(fmin)) + { + lowerBound = root; + fmin = fx; + } + else + { + upperBound = root; + fmax = fx; + } + continue; + } + + if (Math.Abs(step) < accuracy) + { + return root; + } + + // Evaluation + fx = f(root); + lastStep = step; + + // Update Bounds + if (Math.Sign(fx) != Math.Sign(fmin)) + { + upperBound = root; + fmax = fx; + } + else if (Math.Sign(fx) != Math.Sign(fmax)) + { + lowerBound = root; + fmin = fx; + } + else if (Math.Sign(fmin) != Math.Sign(fmax)) + { + return root; + } + } + + throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed"); + } + } +} diff --git a/src/Numerics/RootFinding/FloatingPointRoots.cs b/src/Numerics/RootFinding/FloatingPointRoots.cs index 7e67492b0..5928a311d 100644 --- a/src/Numerics/RootFinding/FloatingPointRoots.cs +++ b/src/Numerics/RootFinding/FloatingPointRoots.cs @@ -39,5 +39,17 @@ public static double OfFunction(Func f, double lowerBound, doubl { return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100); } + + public static double OfFunctionAndDerivative(Func f, Func df, double lowerBound, double upperBound) + { + try + { + return NewtonRaphson.FindRoot(f, df, lowerBound, upperBound, 1e-8, 100); + } + catch (NonConvergenceException) + { + return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100); + } + } } } diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj index c3daafe56..853986018 100644 --- a/src/Portable/Portable.csproj +++ b/src/Portable/Portable.csproj @@ -990,6 +990,9 @@ RootFinding\Algorithms\Brent.cs + + RootFinding\Algorithms\NewtonRaphson.cs + RootFinding\Bracketing.cs diff --git a/src/UnitTests/RootFindingTests/BisectionTest.cs b/src/UnitTests/RootFindingTests/BisectionTest.cs index e86ff19ca..99650e3ca 100644 --- a/src/UnitTests/RootFindingTests/BisectionTest.cs +++ b/src/UnitTests/RootFindingTests/BisectionTest.cs @@ -46,6 +46,9 @@ public void MultipleRoots() Assert.AreEqual(0, f1(Bisection.FindRootExpand(f1, 3, 4, 1e-14))); Assert.AreEqual(-2, Bisection.FindRoot(f1, -5, -1, 1e-14)); Assert.AreEqual(2, Bisection.FindRoot(f1, 1, 4, 1e-14)); + Assert.AreEqual(0, f1(Bisection.FindRoot(x => -f1(x), 0, 5, 1e-14))); + Assert.AreEqual(-2, Bisection.FindRoot(x => -f1(x), -5, -1, 1e-14)); + Assert.AreEqual(2, Bisection.FindRoot(x => -f1(x), 1, 4, 1e-14)); // Roots at 3, 4 Func f2 = x => (x - 3) * (x - 4); @@ -56,6 +59,14 @@ public void MultipleRoots() Assert.AreEqual(3, Bisection.FindRoot(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(Bisection.FindRoot(f1, -5, 5, 1e-14)), 1e-14); + Assert.AreEqual(0, f1(Bisection.FindRoot(f1, -2, 4, 1e-14)), 1e-14); + } + [Test] public void NoRoot() { diff --git a/src/UnitTests/RootFindingTests/BrentTest.cs b/src/UnitTests/RootFindingTests/BrentTest.cs index 6e94d9856..6dea542f5 100644 --- a/src/UnitTests/RootFindingTests/BrentTest.cs +++ b/src/UnitTests/RootFindingTests/BrentTest.cs @@ -45,6 +45,9 @@ public void MultipleRoots() Assert.AreEqual(0, f1(Brent.FindRoot(f1, -5, 5, 1e-14, 100))); Assert.AreEqual(-2, Brent.FindRoot(f1, -5, -1, 1e-14, 100)); Assert.AreEqual(2, Brent.FindRoot(f1, 1, 4, 1e-14, 100)); + Assert.AreEqual(0, f1(Brent.FindRoot(x => -f1(x), -5, 5, 1e-14, 100))); + Assert.AreEqual(-2, Brent.FindRoot(x => -f1(x), -5, -1, 1e-14, 100)); + Assert.AreEqual(2, Brent.FindRoot(x => -f1(x), 1, 4, 1e-14, 100)); // Roots at 3, 4 Func f2 = x => (x - 3)*(x - 4); @@ -55,6 +58,14 @@ public void MultipleRoots() Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001); } + [Test] + public void LocalMinima() + { + Func f1 = x => x * x * x - 2 * x + 2; + Assert.AreEqual(0, f1(Brent.FindRoot(f1, -5, 5, 1e-14, 100)), 1e-14); + Assert.AreEqual(0, f1(Brent.FindRoot(f1, -2, 4, 1e-14, 100)), 1e-14); + } + [Test] public void NoRoot() { diff --git a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs new file mode 100644 index 000000000..85aac18b2 --- /dev/null +++ b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs @@ -0,0 +1,80 @@ +// +// 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.Algorithms; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.RootFindingTests +{ + [TestFixture] + public class NewtonRaphsonTest + { + [Test] + public void MultipleRoots() + { + // Roots at -2, 2 + Func f1 = x => x * x - 4; + Func df1 = x => 2 * x; + Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100))); + Assert.AreEqual(-2, NewtonRaphson.FindRoot(f1, df1, -5, -1, 1e-14, 100)); + Assert.AreEqual(2, NewtonRaphson.FindRoot(f1, df1, 1, 4, 1e-14, 100)); + Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100))); + Assert.AreEqual(-2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100)); + Assert.AreEqual(2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100)); + + // Roots at 3, 4 + Func f2 = x => (x - 3) * (x - 4); + Func df2 = x => 2 * x - 7; + Assert.AreEqual(0, f2(NewtonRaphson.FindRoot(f2, df2, -5, 5, 1e-14, 100))); + Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, -5, 3.5, 1e-14, 100)); + Assert.AreEqual(4, NewtonRaphson.FindRoot(f2, df2, 3.2, 5, 1e-14, 100)); + Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.1, 3.9, 0.001, 50), 0.001); + Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.1, 3.4, 0.001, 50), 0.001); + } + + [Test] + public void LocalMinima() + { + Func f1 = x => x * x * x - 2 * x + 2; + Func df1 = x => 3 * x * x - 2; + Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100))); + Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -2, 4, 1e-14, 100))); + } + + [Test] + public void NoRoot() + { + Func f1 = x => x * x + 4; + Func df1 = x => 2 * x; + Assert.Throws(() => NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50)); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index fcd89f996..685264833 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -772,6 +772,7 @@ +