diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 6d25cd659..ca3e50a50 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -110,7 +110,8 @@
-
+
+
diff --git a/src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs b/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
similarity index 51%
rename from src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs
rename to src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
index f442ce27a..b066833fe 100644
--- a/src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs
+++ b/src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs
@@ -32,20 +32,40 @@
namespace MathNet.Numerics.RootFinding.Algorithms
{
- public static class NewtonRaphson
+ public static class HybridNewtonRaphson
{
/// Find a solution of the equation f(x)=0.
- /// Hybrid Newton-Raphson that when failing falls back to bisection.
+ /// 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 = 1e-8, int maxIterations = 100)
+ public static double FindSingleRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision)
+ {
+ double root;
+ if (TryFindSingleRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, subdivision, out root))
+ {
+ return root;
+ }
+ throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+ }
+
+ /// Find a solution of the equation f(x)=0.
+ /// Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.
+ public static bool TryFindSingleRoot(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
{
double fmin = f(lowerBound);
double fmax = f(upperBound);
- if (fmin == 0.0) return lowerBound;
- if (fmax == 0.0) return upperBound;
+ if (fmin == 0.0)
+ {
+ root = lowerBound;
+ return true;
+ }
+ if (fmax == 0.0)
+ {
+ root = upperBound;
+ return true;
+ }
- double root = 0.5*(lowerBound + upperBound);
+ root = 0.5*(lowerBound + upperBound);
double fx = f(root);
double lastStep = Math.Abs(upperBound - lowerBound);
for (int i = 0; i < maxIterations; i++)
@@ -56,9 +76,18 @@ public static double FindRoot(Func f, Func df, d
double step = fx/dfx;
root -= step;
- if (root < lowerBound || root > upperBound || Math.Abs(2*fx) > Math.Abs(lastStep*dfx))
+ bool overshoot = root > upperBound, undershoot = root < lowerBound;
+ if (overshoot || undershoot || Math.Abs(2*fx) > Math.Abs(lastStep*dfx))
{
- // Newton-Raphson step failed -> bisect instead
+ // Newton-Raphson step failed
+
+ // If same signs, try subdivision to scan for zero crossing intervals
+ if (Math.Sign(fmin) == Math.Sign(fmax) && TryScanForCrossingsWithRoots(f, df, lowerBound, upperBound, accuracy, maxIterations - i - 1, subdivision, out root))
+ {
+ return true;
+ }
+
+ // Bisection
root = 0.5*(upperBound + lowerBound);
fx = f(root);
lastStep = 0.5*Math.Abs(upperBound - lowerBound);
@@ -66,25 +95,35 @@ public static double FindRoot(Func f, Func df, d
{
lowerBound = root;
fmin = fx;
+ if (overshoot)
+ {
+ root = upperBound;
+ fx = fmax;
+ }
}
else
{
upperBound = root;
fmax = fx;
+ if (undershoot)
+ {
+ root = lowerBound;
+ fx = fmin;
+ }
}
continue;
}
-
+
if (Math.Abs(step) < accuracy)
{
- return root;
+ return true;
}
// Evaluation
fx = f(root);
lastStep = step;
- // Update Bounds
+ // Update bounds
if (Math.Sign(fx) != Math.Sign(fmin))
{
upperBound = root;
@@ -97,11 +136,26 @@ public static double FindRoot(Func f, Func df, d
}
else if (Math.Sign(fmin) != Math.Sign(fmax))
{
- return root;
+ return true;
}
}
- throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+ return false;
+ }
+
+ static bool TryScanForCrossingsWithRoots(Func f, Func df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
+ {
+ var zeroCrossings = ZeroCrossingBracketing.FindIntervalsWithin(f, lowerBound, upperBound, subdivision);
+ foreach (Tuple bounds in zeroCrossings)
+ {
+ if (TryFindSingleRoot(f, df, bounds.Item1, bounds.Item2, accuracy, maxIterations, subdivision, out root))
+ {
+ return true;
+ }
+ }
+
+ root = double.NaN;
+ return false;
}
}
}
diff --git a/src/Numerics/RootFinding/Algorithms/ZeroCrossingBracketing.cs b/src/Numerics/RootFinding/Algorithms/ZeroCrossingBracketing.cs
new file mode 100644
index 000000000..689155672
--- /dev/null
+++ b/src/Numerics/RootFinding/Algorithms/ZeroCrossingBracketing.cs
@@ -0,0 +1,44 @@
+using System;
+using System.Collections.Generic;
+
+namespace MathNet.Numerics.RootFinding.Algorithms
+{
+ public static class ZeroCrossingBracketing
+ {
+ public static IEnumerable> FindIntervalsWithin(Func f, double lowerBound, double upperBound, int parts)
+ {
+ // TODO: Consider binary-style search instead of linear scan
+
+ var fmin = f(lowerBound);
+ var fmax = f(upperBound);
+
+ if (Math.Sign(fmin) != Math.Sign(fmax))
+ {
+ yield return new Tuple(lowerBound, upperBound);
+ yield break;
+ }
+
+ double subdiv = (upperBound - lowerBound)/parts;
+ var smin = lowerBound;
+ int sign = Math.Sign(fmin);
+
+ for (int k = 0; k < parts; k++)
+ {
+ var smax = smin + subdiv;
+ var sfmax = f(smax);
+ if (double.IsInfinity(sfmax))
+ {
+ // expand interval to include pole
+ smin = smax;
+ continue;
+ }
+ if (Math.Sign(sfmax) != sign)
+ {
+ yield return new Tuple(smin, smax);
+ sign = Math.Sign(sfmax);
+ }
+ smin = smax;
+ }
+ }
+ }
+}
diff --git a/src/Numerics/RootFinding/FloatingPointRoots.cs b/src/Numerics/RootFinding/FloatingPointRoots.cs
index 5928a311d..4fa14d8c9 100644
--- a/src/Numerics/RootFinding/FloatingPointRoots.cs
+++ b/src/Numerics/RootFinding/FloatingPointRoots.cs
@@ -35,21 +35,22 @@ namespace MathNet.Numerics.RootFinding
{
public static class FloatingPointRoots
{
- public static double OfFunction(Func f, double lowerBound, double upperBound)
+ public static double OfFunction(Func f, double lowerBound, double upperBound, double accuracy = 1e-8)
{
- return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100);
+ return Brent.FindRoot(f, lowerBound, upperBound, accuracy, 100);
}
- public static double OfFunctionAndDerivative(Func f, Func df, double lowerBound, double upperBound)
+ public static double OfFunctionAndDerivative(Func f, Func df, double lowerBound, double upperBound, double accuracy = 1e-8)
{
- try
+ double root;
+ if (HybridNewtonRaphson.TryFindSingleRoot(f, df, lowerBound, upperBound, accuracy, 100, 20, out root))
{
- return NewtonRaphson.FindRoot(f, df, lowerBound, upperBound, 1e-8, 100);
- }
- catch (NonConvergenceException)
- {
- return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100);
+ return root;
}
+
+ return Brent.FindRoot(f, lowerBound, upperBound, accuracy, 100);
+
+ //throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
}
}
}
diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj
index 853986018..7adc1441d 100644
--- a/src/Portable/Portable.csproj
+++ b/src/Portable/Portable.csproj
@@ -990,8 +990,11 @@
RootFinding\Algorithms\Brent.cs
-
- RootFinding\Algorithms\NewtonRaphson.cs
+
+ RootFinding\Algorithms\HybridNewtonRaphson.cs
+
+
+ RootFinding\Algorithms\ZeroCrossingBracketing.cs
RootFinding\Bracketing.cs
diff --git a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
index 85aac18b2..5a091f699 100644
--- a/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
+++ b/src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs
@@ -29,6 +29,7 @@
//
using System;
+using MathNet.Numerics.RootFinding;
using MathNet.Numerics.RootFinding.Algorithms;
using NUnit.Framework;
@@ -43,21 +44,21 @@ 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));
+ Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(-2, HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, -1, 1e-14, 100, 20));
+ Assert.AreEqual(2, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 4, 1e-14, 100, 20));
+ Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(-2, HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100, 20));
+ Assert.AreEqual(2, HybridNewtonRaphson.FindSingleRoot(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(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);
+ Assert.AreEqual(0, f2(HybridNewtonRaphson.FindSingleRoot(f2, df2, -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, -5, 3.5, 1e-14, 100, 20));
+ Assert.AreEqual(4, HybridNewtonRaphson.FindSingleRoot(f2, df2, 3.2, 5, 1e-14, 100, 20));
+ Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2.1, 3.9, 0.001, 50, 20), 0.001);
+ Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2.1, 3.4, 0.001, 50, 20), 0.001);
}
[Test]
@@ -65,8 +66,31 @@ 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)));
+ Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+ Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(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, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 2, 1e-14, 100, 20));
+ Assert.AreEqual(1.5, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 6, 1e-14, 100, 20));
+ Assert.AreEqual(1.5, FloatingPointRoots.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, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2, 3, 1e-14, 100, 20));
+ Assert.AreEqual(2.5, HybridNewtonRaphson.FindSingleRoot(f2, df2, -2, 3, 1e-14, 100, 20));
+ Assert.AreEqual(2.5, FloatingPointRoots.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), HybridNewtonRaphson.FindSingleRoot(f3, df3, -2, -1, 1e-14, 100, 20), 1e-14);
+ Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, 1, 1.99, 1e-14, 100, 20));
+ Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, -1.5, 1.99, 1e-14, 100, 20));
+ Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, 1, 6, 1e-14, 100, 20));
}
[Test]
@@ -74,7 +98,7 @@ 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));
+ Assert.Throws(() => HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 50, 20));
}
}
}