From bd8742aec32035c5f10cc83b3facc6b8281d924f Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 19 Oct 2018 19:23:17 +0200 Subject: [PATCH] Polynomial: more integration into FindRoots --- MathNet.Numerics.sln.DotSettings | 1 + .../RootFindingTests/CubicTest.cs | 31 +++++++++++++ src/Numerics/FindRoots.cs | 22 +++++++--- src/Numerics/Polynomial.cs | 43 ++++++++++++------- 4 files changed, 75 insertions(+), 22 deletions(-) diff --git a/MathNet.Numerics.sln.DotSettings b/MathNet.Numerics.sln.DotSettings index 14b3d43d..dd611a16 100644 --- a/MathNet.Numerics.sln.DotSettings +++ b/MathNet.Numerics.sln.DotSettings @@ -81,6 +81,7 @@ OTHER DEALINGS IN THE SOFTWARE. True True True + True True diff --git a/src/Numerics.Tests/RootFindingTests/CubicTest.cs b/src/Numerics.Tests/RootFindingTests/CubicTest.cs index 43cdfc73..b8763dc0 100644 --- a/src/Numerics.Tests/RootFindingTests/CubicTest.cs +++ b/src/Numerics.Tests/RootFindingTests/CubicTest.cs @@ -97,6 +97,19 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.That(roots.Item3.Imaginary, Is.EqualTo(0).Within(1e-14)); } + [TestCase(6.0, -5.0, -2.0, 1.0, 3.0, -2.0, 1.0)] + public void ComplexRoots_TripleReal_AsPolynomial(double d, double c, double b, double a, double x1, double x2, double x3) + { + var roots = FindRoots.Polynomial(new[] { d, c, b, a }); + Assert.That(roots.Length, Is.EqualTo(3)); + Assert.That(roots[0].Real, Is.EqualTo(x1).Within(1e-14)); + Assert.That(roots[0].Imaginary, Is.EqualTo(0).Within(1e-14)); + Assert.That(roots[1].Real, Is.EqualTo(x2).Within(1e-14)); + Assert.That(roots[1].Imaginary, Is.EqualTo(0).Within(1e-14)); + Assert.That(roots[2].Real, Is.EqualTo(x3).Within(1e-14)); + Assert.That(roots[2].Imaginary, Is.EqualTo(0).Within(1e-14)); + } + [TestCase(-350.0, 162.0, -30.0, 2.0)] [TestCase(6.0, -5.0, -2.0, 1.0)] [TestCase(1.0, 5.0, 2.0, 1.0)] @@ -113,5 +126,23 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests Assert.That(Polynomial.Evaluate(roots.Item3, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); Assert.That(Polynomial.Evaluate(roots.Item3, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); } + + [TestCase(-350.0, 162.0, -30.0, 2.0)] + [TestCase(6.0, -5.0, -2.0, 1.0)] + [TestCase(1.0, 5.0, 2.0, 1.0)] + [TestCase(1.0, 5.0, 0.0, 1.0)] + [TestCase(1.0, 0.0, 2.0, 1.0)] + [TestCase(0.0, 0.0, 0.0, 2.0)] + public void ComplexRootsAreRoots_AsPolynomial(double d, double c, double b, double a) + { + var roots = FindRoots.Polynomial(new[] { d, c, b, a }); + Assert.That(roots.Length, Is.EqualTo(3)); + Assert.That(Polynomial.Evaluate(roots[0], d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); + Assert.That(Polynomial.Evaluate(roots[0], d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); + Assert.That(Polynomial.Evaluate(roots[1], d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); + Assert.That(Polynomial.Evaluate(roots[1], d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); + Assert.That(Polynomial.Evaluate(roots[2], d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); + Assert.That(Polynomial.Evaluate(roots[2], d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); + } } } diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs index 4d101d54..73008682 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -114,11 +114,21 @@ namespace MathNet.Numerics /// /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix /// - /// the values for the polynomial in ascending order e.G new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" - /// the roots of the polynomial - public static Complex[] Polynomial(double[] poly) + /// The coefficients of the polynomial in ascending order, e.g. new double[] {5, 0, 2} = "5 + 0 x^1 + 2 x^2" + /// The roots of the polynomial + public static Complex[] Polynomial(double[] coefficients) { - return new Polynomial(poly).Roots(); + return new Polynomial(coefficients).Roots(); + } + + /// + /// Find all roots of a polynomial by calculating the characteristic polynomial of the companion matrix + /// + /// The polynomial. + /// The roots of the polynomial + public static Complex[] Polynomial(Polynomial polynomial) + { + return polynomial.Roots(); } /// @@ -139,7 +149,7 @@ namespace MathNet.Numerics double location = 0.5*(intervalBegin + intervalEnd); double scale = 0.5*(intervalEnd - intervalBegin); - // evaluate first kind chebyshev nodes + // evaluate first kind chebychev nodes double angleFactor = Constants.Pi/(2*degree); var samples = new double[degree]; @@ -168,7 +178,7 @@ namespace MathNet.Numerics double location = 0.5*(intervalBegin + intervalEnd); double scale = 0.5*(intervalEnd - intervalBegin); - // evaluate second kind chebyshev nodes + // evaluate second kind chebychev nodes double angleFactor = Constants.Pi/(degree + 1); var samples = new double[degree]; diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index edceedfc..fae344bf 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -76,10 +76,21 @@ namespace MathNet.Numerics /// Create a constant polynomial. /// Example: 3.0 -> "p : x -> 3.0" /// - /// just the "x^0" part + /// The coefficient of the "x^0" monomial. public Polynomial(double coefficient) { - Coefficients = coefficient == 0.0 ? new double[0] : new[] { coefficient }; + if (coefficient == 0.0) + { +#if NET40 + Coefficients = new double[0]; +#else + Coefficients = Array.Empty(); +#endif + } + else + { + Coefficients = new[] { coefficient }; + } } /// @@ -101,6 +112,17 @@ namespace MathNet.Numerics { } + public static Polynomial Zero => new Polynomial(); + + /// + /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k + /// + public static Polynomial Fit(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) + { + var coefficients = Numerics.Fit.Polynomial(x, y, order, method); + return new Polynomial(coefficients); + } + static int EvaluateDegree(double[] coefficients) { for (int i = coefficients.Length - 1; i >= 0; i--) @@ -114,22 +136,11 @@ namespace MathNet.Numerics return -1; } - public static Polynomial Zero => new Polynomial(); - - /// - /// Least-Squares fitting the points (x,y) to a k-order polynomial y : x -> p0 + p1*x + p2*x^2 + ... + pk*x^k - /// - public static Polynomial Fit(double[] x, double[] y, int order, DirectRegressionMethod method = DirectRegressionMethod.QR) - { - var coefficients = Numerics.Fit.Polynomial(x, y, order, method); - return new Polynomial(coefficients); - } - #region Evaluation /// /// Evaluate a polynomial at point x. - /// Coefficients are ordered by power with power k at index k. + /// Coefficients are ordered ascending by power with power k at index k. /// Example: coefficients [3,-1,2] represent y=2x^2-x+3. /// /// The location where to evaluate the polynomial at. @@ -148,7 +159,7 @@ namespace MathNet.Numerics /// /// Evaluate a polynomial at point x. - /// Coefficients are ordered by power with power k at index k. + /// Coefficients are ordered ascending by power with power k at index k. /// Example: coefficients [3,-1,2] represent y=2x^2-x+3. /// /// The location where to evaluate the polynomial at. @@ -167,7 +178,7 @@ namespace MathNet.Numerics /// /// Evaluate a polynomial at point x. - /// Coefficients are ordered by power with power k at index k. + /// Coefficients are ordered ascending by power with power k at index k. /// Example: coefficients [3,-1,2] represent y=2x^2-x+3. /// /// The location where to evaluate the polynomial at.