diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 9cde96e1..6b3a77b8 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -111,6 +111,7 @@
+
diff --git a/src/Numerics/RootFinding/Bisection.cs b/src/Numerics/RootFinding/Bisection.cs
index c0b907cc..729edeb8 100644
--- a/src/Numerics/RootFinding/Bisection.cs
+++ b/src/Numerics/RootFinding/Bisection.cs
@@ -33,6 +33,9 @@ using MathNet.Numerics.Properties;
namespace MathNet.Numerics.RootFinding
{
+ ///
+ /// Bisection root-finding algorithm.
+ ///
public static class Bisection
{
/// Find a solution of the equation f(x)=0.
diff --git a/src/Numerics/RootFinding/Brent.cs b/src/Numerics/RootFinding/Brent.cs
index 524e600f..00173dd8 100644
--- a/src/Numerics/RootFinding/Brent.cs
+++ b/src/Numerics/RootFinding/Brent.cs
@@ -32,6 +32,10 @@ using System;
namespace MathNet.Numerics.RootFinding
{
+ ///
+ /// Algorithm by by Brent, Van Wijngaarden, Dekker et al.
+ /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press
+ ///
public static class Brent
{
/// Find a solution of the equation f(x)=0.
@@ -41,10 +45,6 @@ namespace MathNet.Numerics.RootFinding
/// Desired accuracy. The root will be refined until the accuracy or the maximum number of iterations is reached.
/// Maximum number of iterations. Usually 100.
/// Returns the root with the specified accuracy.
- ///
- /// Algorithm by by Brent, Van Wijngaarden, Dekker et al.
- /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press
- ///
///
public static double FindRoot(Func f, double lowerBound, double upperBound, double accuracy, int maxIterations)
{
@@ -64,10 +64,6 @@ namespace MathNet.Numerics.RootFinding
/// Maximum number of iterations. Usually 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.
- ///
- /// Algorithm by by Brent, Van Wijngaarden, Dekker et al.
- /// Implementation inspired by Press, Teukolsky, Vetterling, and Flannery, "Numerical Recipes in C", 2nd edition, Cambridge University Press
- ///
public static bool TryFindRoot(Func f, double lowerBound, double upperBound, double accuracy, int maxIterations, out double root)
{
double fmin = f(lowerBound);
diff --git a/src/Numerics/RootFinding/NewtonRaphson.cs b/src/Numerics/RootFinding/NewtonRaphson.cs
new file mode 100644
index 00000000..98d13561
--- /dev/null
+++ b/src/Numerics/RootFinding/NewtonRaphson.cs
@@ -0,0 +1,94 @@
+//
+// 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
+{
+ ///
+ /// Pure Newton-Raphson 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 NewtonRaphson
+ {
+ /// Find a solution of the equation f(x)=0.
+ /// The function to find roots from.
+ /// The first derivative of the function to find roots from.
+ /// Initial guess of the root.
+ /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval.
+ /// The high 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.
+ /// Returns the root with the specified accuracy.
+ ///
+ public static double FindRoot(Func f, Func df, double initialGuess, double lowerBound, double upperBound, double accuracy, int maxIterations)
+ {
+ double root;
+ if (TryFindRoot(f, df, initialGuess, lowerBound, upperBound, accuracy, maxIterations, out root))
+ {
+ return root;
+ }
+ throw new NonConvergenceException("The algorithm failed or has exceeded the number of iterations allowed. Consider to use RobustNewtonRaphson instead.");
+ }
+
+ /// Find a solution of the equation f(x)=0.
+ /// The function to find roots from.
+ /// The first derivative of the function to find roots from.
+ /// Initial guess of the root.
+ /// The low value of the range where the root is supposed to be. Aborts if it leaves the interval.
+ /// The high 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, Func df, double initialGuess, double lowerBound, double upperBound, double accuracy, int maxIterations, out double root)
+ {
+ root = initialGuess;
+ for (int i = 0; i < maxIterations && root >= lowerBound && root <= upperBound; i++)
+ {
+ // Evaluation
+ double fx = f(root);
+ double dfx = df(root);
+
+ // Netwon-Raphson step
+ double step = fx/dfx;
+ root -= step;
+
+ if (Math.Abs(step) < accuracy && Math.Abs(fx) < accuracy)
+ {
+ return true;
+ }
+ }
+
+ return false;
+ }
+ }
+}
diff --git a/src/Numerics/RootFinding/RobustNewtonRaphson.cs b/src/Numerics/RootFinding/RobustNewtonRaphson.cs
index 5ee9a374..13424418 100644
--- a/src/Numerics/RootFinding/RobustNewtonRaphson.cs
+++ b/src/Numerics/RootFinding/RobustNewtonRaphson.cs
@@ -32,6 +32,10 @@ using System;
namespace MathNet.Numerics.RootFinding
{
+ ///
+ /// Robust Newton-Raphson root-finding algorithm that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.
+ ///
+ ///
public static class RobustNewtonRaphson
{
/// Find a solution of the equation f(x)=0.
@@ -43,7 +47,6 @@ namespace MathNet.Numerics.RootFinding
/// Maximum number of iterations. Example: 100.
/// How many parts an interval should be split into for zero crossing scanning in case of lacking bracketing. Example: 20.
/// Returns the root with the specified accuracy.
- /// Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.
///
public static double FindRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision)
{
@@ -65,7 +68,6 @@ namespace MathNet.Numerics.RootFinding
/// How many parts an interval should be split into for zero crossing scanning in case of lacking bracketing. Example: 20.
/// 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.
- /// Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.
public static bool TryFindRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
{
double fmin = f(lowerBound);
diff --git a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
index 155284af..4b841938 100644
--- a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
+++ b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
@@ -43,21 +43,21 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
// Roots at -2, 2
Func f1 = x => x * x - 4;
Func df1 = x => 2 * x;
- Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
- Assert.AreEqual(-2, RobustNewtonRaphson.FindRoot(f1, df1, -5, -1, 1e-14, 100, 20));
- Assert.AreEqual(2, RobustNewtonRaphson.FindRoot(f1, df1, 1, 4, 1e-14, 100, 20));
- Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100, 20)));
- Assert.AreEqual(-2, RobustNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100, 20));
- Assert.AreEqual(2, RobustNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100, 20));
+ Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, 0.5, -5, 6, 1e-14, 100)));
+ Assert.AreEqual(-2, NewtonRaphson.FindRoot(f1, df1, -3.0, -5, -1, 1e-14, 100));
+ Assert.AreEqual(2, NewtonRaphson.FindRoot(f1, df1, 2.5, 1, 4, 1e-14, 100));
+ Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 0.6, -5, 6, 1e-14, 100)));
+ Assert.AreEqual(-2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -3, -5, -1, 1e-14, 100));
+ Assert.AreEqual(2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 2.5, 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(RobustNewtonRaphson.FindRoot(f2, df2, -5, 5, 1e-14, 100, 20)));
- Assert.AreEqual(3, RobustNewtonRaphson.FindRoot(f2, df2, -5, 3.5, 1e-14, 100, 20));
- Assert.AreEqual(4, RobustNewtonRaphson.FindRoot(f2, df2, 3.2, 5, 1e-14, 100, 20));
- Assert.AreEqual(3, RobustNewtonRaphson.FindRoot(f2, df2, 2.1, 3.9, 0.001, 50, 20), 0.001);
- Assert.AreEqual(3, RobustNewtonRaphson.FindRoot(f2, df2, 2.1, 3.4, 0.001, 50, 20), 0.001);
+ Assert.AreEqual(0, f2(NewtonRaphson.FindRoot(f2, df2, 0.5, -5, 6, 1e-14, 100)));
+ Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, -0.75, -5, 3.5, 1e-14, 100));
+ Assert.AreEqual(4, NewtonRaphson.FindRoot(f2, df2, 4.1, 3.2, 5, 1e-14, 100));
+ Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 3, 2.1, 3.9, 0.001, 50), 0.001);
+ Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.75, 2.1, 3.4, 0.001, 50), 0.001);
}
[Test]
@@ -65,38 +65,8 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
{
Func f1 = x => x * x * x - 2 * x + 2;
Func df1 = x => 3 * x * x - 2;
- Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
- Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(f1, df1, -2, 4, 1e-14, 100, 20)));
- }
-
- [Test]
- public void Pole()
- {
- Func f1 = x => 1/(x - 2) + 2;
- Func df1 = x => -1/(x*x - 4*x + 4);
- Assert.AreEqual(1.5, RobustNewtonRaphson.FindRoot(f1, df1, 1, 2, 1e-14, 100, 20));
- Assert.AreEqual(1.5, RobustNewtonRaphson.FindRoot(f1, df1, 1, 6, 1e-14, 100, 20));
- Assert.AreEqual(1.5, FindRoots.OfFunctionAndDerivative(f1, df1, 1, 6));
-
- Func f2 = x => -1/(x - 2) + 2;
- Func df2 = x => 1/(x*x - 4*x + 4);
- Assert.AreEqual(2.5, RobustNewtonRaphson.FindRoot(f2, df2, 2, 3, 1e-14, 100, 20));
- Assert.AreEqual(2.5, RobustNewtonRaphson.FindRoot(f2, df2, -2, 3, 1e-14, 100, 20));
- Assert.AreEqual(2.5, FindRoots.OfFunctionAndDerivative(f2, df2, -2, 3));
-
- Func f3 = x => 1/(x - 2) + x + 2;
- Func df3 = x => -1/(x*x - 4*x + 4) + 1;
- Assert.AreEqual(-Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, -2, -1, 1e-14, 100, 20), 1e-14);
- Assert.AreEqual(Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, 1, 1.99, 1e-14, 100, 20));
- Assert.AreEqual(Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, -1.5, 1.99, 1e-14, 100, 20));
- Assert.AreEqual(Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, 1, 6, 1e-14, 100, 20));
-
- Func f4 = x => 1/(2 - x) - x + 6;
- Func df4 = x => 1/(x*x - 4*x + 4) - 1;
- Assert.AreEqual(4 + Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, 5, 6, 1e-14, 100, 20), 1e-14);
- Assert.AreEqual(4 - Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, 2.01, 3, 1e-14, 100, 20));
- Assert.AreEqual(4 - Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, 2.01, 5, 1e-14, 100, 20));
- Assert.AreEqual(4 - Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, -2, 4, 1e-14, 100, 20));
+ Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, 0.5, -5, 6, 1e-14, 100)));
+ //Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, 1, -2, 4, 1e-14, 100)));
}
[Test]
@@ -105,15 +75,15 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
// with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5
Func f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3);
Func df1 = x => Evaluate.Polynomial(x, 5, 8, 9);
- Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10, 100, 20), 1e-6);
- Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10, 100, 20), 1e-6);
+ Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, -1.5, -2, -1, 1e-10, 100), 1e-6);
+ Assert.AreEqual(-1.265328088928, NewtonRaphson.FindRoot(f1, df1, 0, -5, 5, 1e-10, 100), 1e-6);
// real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50
Func f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2);
Func df2 = x => Evaluate.Polynomial(x, -50, 8, 6);
- Assert.AreEqual(-6.1466562197069, RobustNewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10, 100, 20), 1e-6);
- Assert.AreEqual(0.12124737195841, RobustNewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10, 100, 20), 1e-6);
- Assert.AreEqual(4.0254088477485, RobustNewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10, 100, 20), 1e-6);
+ Assert.AreEqual(-6.1466562197069, NewtonRaphson.FindRoot(f2, df2, -6.5, -8, -5, 1e-10, 100), 1e-6);
+ Assert.AreEqual(0.12124737195841, NewtonRaphson.FindRoot(f2, df2, 0, -1, 1, 1e-10, 100), 1e-6);
+ Assert.AreEqual(4.0254088477485, NewtonRaphson.FindRoot(f2, df2, 4, 3, 5, 1e-10, 100), 1e-6);
}
[Test]
@@ -121,7 +91,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
{
Func f1 = x => x * x + 4;
Func df1 = x => 2 * x;
- Assert.Throws(() => RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50, 20));
+ Assert.Throws(() => NewtonRaphson.FindRoot(f1, df1, 0, -5, 5, 1e-14, 50));
}
}
}
diff --git a/src/UnitTests/RootFindingTests/RobustNewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/RobustNewtonRaphsonTest.cs
new file mode 100644
index 00000000..9f649ea1
--- /dev/null
+++ b/src/UnitTests/RootFindingTests/RobustNewtonRaphsonTest.cs
@@ -0,0 +1,127 @@
+//
+// 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]
+ public class RobustNewtonRaphsonTest
+ {
+ [Test]
+ public void MultipleRoots()
+ {
+ // Roots at -2, 2
+ Func f1 = x => x * x - 4;
+ Func df1 = x => 2 * x;
+ Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(-2, RobustNewtonRaphson.FindRoot(f1, df1, -5, -1, 1e-14, 100, 20));
+ Assert.AreEqual(2, RobustNewtonRaphson.FindRoot(f1, df1, 1, 4, 1e-14, 100, 20));
+ Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(-2, RobustNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100, 20));
+ Assert.AreEqual(2, RobustNewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100, 20));
+
+ // Roots at 3, 4
+ Func f2 = x => (x - 3) * (x - 4);
+ Func df2 = x => 2 * x - 7;
+ Assert.AreEqual(0, f2(RobustNewtonRaphson.FindRoot(f2, df2, -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(3, RobustNewtonRaphson.FindRoot(f2, df2, -5, 3.5, 1e-14, 100, 20));
+ Assert.AreEqual(4, RobustNewtonRaphson.FindRoot(f2, df2, 3.2, 5, 1e-14, 100, 20));
+ Assert.AreEqual(3, RobustNewtonRaphson.FindRoot(f2, df2, 2.1, 3.9, 0.001, 50, 20), 0.001);
+ Assert.AreEqual(3, RobustNewtonRaphson.FindRoot(f2, df2, 2.1, 3.4, 0.001, 50, 20), 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(RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(0, f1(RobustNewtonRaphson.FindRoot(f1, df1, -2, 4, 1e-14, 100, 20)));
+ }
+
+ [Test]
+ public void Pole()
+ {
+ Func f1 = x => 1/(x - 2) + 2;
+ Func df1 = x => -1/(x*x - 4*x + 4);
+ Assert.AreEqual(1.5, RobustNewtonRaphson.FindRoot(f1, df1, 1, 2, 1e-14, 100, 20));
+ Assert.AreEqual(1.5, RobustNewtonRaphson.FindRoot(f1, df1, 1, 6, 1e-14, 100, 20));
+ Assert.AreEqual(1.5, FindRoots.OfFunctionAndDerivative(f1, df1, 1, 6));
+
+ Func f2 = x => -1/(x - 2) + 2;
+ Func df2 = x => 1/(x*x - 4*x + 4);
+ Assert.AreEqual(2.5, RobustNewtonRaphson.FindRoot(f2, df2, 2, 3, 1e-14, 100, 20));
+ Assert.AreEqual(2.5, RobustNewtonRaphson.FindRoot(f2, df2, -2, 3, 1e-14, 100, 20));
+ Assert.AreEqual(2.5, FindRoots.OfFunctionAndDerivative(f2, df2, -2, 3));
+
+ Func f3 = x => 1/(x - 2) + x + 2;
+ Func df3 = x => -1/(x*x - 4*x + 4) + 1;
+ Assert.AreEqual(-Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, -2, -1, 1e-14, 100, 20), 1e-14);
+ Assert.AreEqual(Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, 1, 1.99, 1e-14, 100, 20));
+ Assert.AreEqual(Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, -1.5, 1.99, 1e-14, 100, 20));
+ Assert.AreEqual(Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f3, df3, 1, 6, 1e-14, 100, 20));
+
+ Func f4 = x => 1/(2 - x) - x + 6;
+ Func df4 = x => 1/(x*x - 4*x + 4) - 1;
+ Assert.AreEqual(4 + Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, 5, 6, 1e-14, 100, 20), 1e-14);
+ Assert.AreEqual(4 - Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, 2.01, 3, 1e-14, 100, 20));
+ Assert.AreEqual(4 - Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, 2.01, 5, 1e-14, 100, 20));
+ Assert.AreEqual(4 - Math.Sqrt(3), RobustNewtonRaphson.FindRoot(f4, df4, -2, 4, 1e-14, 100, 20));
+ }
+
+ [Test]
+ public void Cubic()
+ {
+ // with complex roots (looking for the real root only): 3x^3 + 4x^2 + 5x + 6, derivative 9x^2 + 8x + 5
+ Func f1 = x => Evaluate.Polynomial(x, 6, 5, 4, 3);
+ Func df1 = x => Evaluate.Polynomial(x, 5, 8, 9);
+ Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -2, -1, 1e-10, 100, 20), 1e-6);
+ Assert.AreEqual(-1.265328088928, RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-10, 100, 20), 1e-6);
+
+ // real roots only: 2x^3 + 4x^2 - 50x + 6, derivative 6x^2 + 8x - 50
+ Func f2 = x => Evaluate.Polynomial(x, 6, -50, 4, 2);
+ Func df2 = x => Evaluate.Polynomial(x, -50, 8, 6);
+ Assert.AreEqual(-6.1466562197069, RobustNewtonRaphson.FindRoot(f2, df2, -8, -5, 1e-10, 100, 20), 1e-6);
+ Assert.AreEqual(0.12124737195841, RobustNewtonRaphson.FindRoot(f2, df2, -1, 1, 1e-10, 100, 20), 1e-6);
+ Assert.AreEqual(4.0254088477485, RobustNewtonRaphson.FindRoot(f2, df2, 3, 5, 1e-10, 100, 20), 1e-6);
+ }
+
+ [Test]
+ public void NoRoot()
+ {
+ Func f1 = x => x * x + 4;
+ Func df1 = x => 2 * x;
+ Assert.Throws(() => RobustNewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50, 20));
+ }
+ }
+}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index bbe923c8..24e0b96b 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -773,6 +773,7 @@
+