From 886de41952b4b6a14fb3a79c320ce311f4ccb75b Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 25 Jul 2013 13:23:22 +0200 Subject: [PATCH] RootFinding: quadratic case, inline docs --- src/Numerics/FindRoots.cs | 38 ++++++++++++++ .../RootFindingTests/FindRootsTest.cs | 49 ++++++++++++++++++- 2 files changed, 86 insertions(+), 1 deletion(-) 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); + } } }