diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs
index 183d6651..bc7a083d 100644
--- a/src/Numerics/FindRoots.cs
+++ b/src/Numerics/FindRoots.cs
@@ -32,10 +32,22 @@ using System;
using MathNet.Numerics.Properties;
using MathNet.Numerics.RootFinding;
+#if NOSYSNUMERICS
+ using Complex = Numerics.Complex;
+#else
+ using Complex = System.Numerics.Complex;
+#endif
+
namespace MathNet.Numerics
{
public static class FindRoots
{
+ /// Find a solution of the equation f(x)=0.
+ /// The function to find roots from.
+ /// The low value of the range where the root is supposed to be.
+ /// The high value of the range where the root is supposed to be.
+ /// 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.
public static double OfFunction(Func f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
{
double root;
@@ -53,6 +65,13 @@ namespace MathNet.Numerics
throw new NonConvergenceException(Resources.RootFindingFailed);
}
+ /// 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.
+ /// The low value of the range where the root is supposed to be.
+ /// The high value of the range where the root is supposed to be.
+ /// 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.
public static double OfFunctionDerivative(Func f, Func df, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
{
double root;
@@ -69,5 +88,24 @@ namespace MathNet.Numerics
throw new NonConvergenceException(Resources.RootFindingFailed);
}
+
+ ///
+ /// Find both complex roots of the quadratic equation c + b*x + a*x^2 = 0.
+ /// Note the special coefficient order ascending by exponent (consistent with polynomials).
+ ///
+ public static Tuple Quadratic(double c, double b, double a)
+ {
+ if (b == 0d)
+ {
+ var t = new Complex(-c/a, 0d).SquareRoot();
+ return new Tuple(t, -t);
+ }
+
+ var q = b > 0d
+ ? -0.5*(b + new Complex(b*b - 4*a*c, 0d).SquareRoot())
+ : -0.5*(b - new Complex(b*b - 4*a*c, 0d).SquareRoot());
+
+ return new Tuple(q/a, c/q);
+ }
}
}
diff --git a/src/UnitTests/RootFindingTests/FindRootsTest.cs b/src/UnitTests/RootFindingTests/FindRootsTest.cs
index 60ea98f2..ddb22a8f 100644
--- a/src/UnitTests/RootFindingTests/FindRootsTest.cs
+++ b/src/UnitTests/RootFindingTests/FindRootsTest.cs
@@ -29,7 +29,7 @@
//
using System;
-using MathNet.Numerics.RootFinding;
+using System.Numerics;
using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.RootFindingTests
@@ -154,5 +154,52 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests
Assert.AreEqual(0.842524411168525, x, 1e-2);
Assert.AreEqual(0, f1(x), 1e-13);
}
+
+ private void AssertComplexEqual(Complex expected, Complex actual, double delta)
+ {
+ Assert.AreEqual(expected.Real, actual.Real, delta);
+ Assert.AreEqual(expected.Imaginary, actual.Imaginary, delta);
+ }
+
+ [TestCase(2d, 3d, 1d)]
+ [TestCase(-2d, 3d, 2d)]
+ [TestCase(2d, -3d, -2d)]
+ [TestCase(3d, 0d, 2d)]
+ public void QuadraticExpanded(double u, double v, double t)
+ {
+ // t*(x+u)*(x+v) = t*u*v + t*(u+v)*x + t*x^2
+ double c = t*u*v, b = t*(u + v), a = t;
+ var x = FindRoots.Quadratic(c, b, a);
+ Complex x1 = x.Item1, x2 = x.Item2;
+
+ // Expected Roots: -u, -v
+ Complex r1 = new Complex(-u, 0d), r2 = new Complex(-v, 0d);
+ Assert.IsTrue(x1.AlmostEqualWithError(r1, 1e-14) && x2.AlmostEqualWithError(r2, 1e-14)
+ || x1.AlmostEqualWithError(r2, 1e-14) && x2.AlmostEqualWithError(r1, 1e-14));
+
+ // Verify they really are roots
+ AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14);
+ AssertComplexEqual(Complex.Zero, c + b*x2 + a*x2*x2, 1e-14);
+ }
+
+ [TestCase(1d, 1d, 1d, -0.5, -0.866025403784439, -0.5, 0.866025403784439)]
+ [TestCase(3d, 4d, 5d, -0.4, -0.663324958071080, -0.4, 0.663324958071080)]
+ [TestCase(-4d, 2d, 1.5d, -2.43050087404306, 0d, 1.09716754070973, 0d)]
+ [TestCase(4d, 0d, 1.5d, 0d, -1.63299316185545, 0d, 1.63299316185545)]
+ public void QuadraticExpected(double c, double b, double a, double x1R, double x1I, double x2R, double x2I)
+ {
+ // c + b*x + a*x^2
+ var x = FindRoots.Quadratic(c, b, a);
+ Complex x1 = x.Item1, x2 = x.Item2;
+
+ // Expected roots?
+ Complex r1 = new Complex(x1R, x1I), r2 = new Complex(x2R, x2I);
+ Assert.IsTrue(x1.AlmostEqualWithError(r1, 1e-14) && x2.AlmostEqualWithError(r2, 1e-14)
+ || x1.AlmostEqualWithError(r2, 1e-14) && x2.AlmostEqualWithError(r1, 1e-14));
+
+ // Verify they really are roots
+ AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14);
+ AssertComplexEqual(Complex.Zero, c + b*x2 + a*x2*x2, 1e-14);
+ }
}
}