From b0ed692ab2a5b4cff369e477cad6637eb19ffe6e Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 22 Feb 2014 19:10:23 +0100 Subject: [PATCH] RootFinding: new implementation for complex cubic roots, tests --- src/Numerics/FindRoots.cs | 4 +- src/Numerics/RootFinding/Cubic.cs | 81 ++++++++++++++++----- src/UnitTests/RootFindingTests/CubicTest.cs | 41 +++++++++-- 3 files changed, 99 insertions(+), 27 deletions(-) diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs index 5860c987..5dfeec0f 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2014 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -112,7 +112,7 @@ namespace MathNet.Numerics /// public static Tuple Cubic(double d, double c, double b, double a) { - return RootFinding.Cubic.Roots(b/a, c/a, d/a); + return RootFinding.Cubic.Roots(d, c, b, a); } /// diff --git a/src/Numerics/RootFinding/Cubic.cs b/src/Numerics/RootFinding/Cubic.cs index f6d81479..6b206a31 100644 --- a/src/Numerics/RootFinding/Cubic.cs +++ b/src/Numerics/RootFinding/Cubic.cs @@ -1,4 +1,34 @@ -using System; +// +// 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-2014 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; #if !NOSYSNUMERICS using Complex = System.Numerics.Complex; @@ -34,7 +64,11 @@ namespace MathNet.Numerics.RootFinding return Math.Pow(Math.Abs(n), 1d / 3d) * Math.Sign(n); } - public static Tuple RealRoots(double a2, double a1, double a0) + /// + /// Find all real-valued roots of the cubic equation a0 + a1*x + a2*x^2 + x^3 = 0. + /// Note the special coefficient order ascending by exponent (consistent with polynomials). + /// + public static Tuple RealRoots(double a0, double a1, double a2) { double Q, R; QR(a2, a1, a0, out Q, out R); @@ -68,27 +102,38 @@ namespace MathNet.Numerics.RootFinding return new Tuple(x1, x2, x3); } - public static Tuple Roots(double a2, double a1, double a0) + /// + /// Find all three complex roots of the cubic equation d + c*x + b*x^2 + a*x^3 = 0. + /// Note the special coefficient order ascending by exponent (consistent with polynomials). + /// + public static Tuple Roots(double d, double c, double b, double a) { - // use eqn (54)-(56) - double Q, R; - QR(a2, a1, a0, out Q, out R); + double A = b*b - 3*a*c; + double B = 2*b*b*b - 9*a*b*c + 27*a*a*d; + double s = -1/(3*a); - var D = Q * Q * Q + R * R; + double D = (B*B - 4*A*A*A)/(-27*a*a); + if (D == 0d) + { + if (A == 0d) + { + var u = new Complex(s*b, 0d); + return new Tuple(u, u, u); + } - var rootD = Complex.Sqrt(D); - var S = Complex.Pow(R + rootD, 1d / 3d); - var T = Complex.Pow(R - rootD, 1d / 3d); - var shift = -a2 / 3d; - var sharedI = 0.5 * Complex.ImaginaryOne * Constants.Sqrt3 * (S - T); + var v = new Complex((9*a*d - b*c)/(2*A), 0d); + var w = new Complex((4*a*b*c - 9*a*a*d - b*b*b)/(a*A), 0d); + return new Tuple(v, v, w); + } - var x1 = shift + (S + T); - var x2 = shift - 0.5 * (S + T); - var x3 = x2; - x2 += sharedI; - x3 -= sharedI; + var C = (A == 0) + ? new Complex(B, 0d).CubicRoots() + : ((B + Complex.Sqrt(B*B - 4*A*A*A))/2).CubicRoots(); - return new Tuple(x1, x2, x3); + return new Tuple( + s*(b + C.Item1 + A/C.Item1), + s*(b + C.Item2 + A/C.Item2), + s*(b + C.Item3 + A/C.Item3)); } } } diff --git a/src/UnitTests/RootFindingTests/CubicTest.cs b/src/UnitTests/RootFindingTests/CubicTest.cs index c8d71a3f..b0e66abf 100644 --- a/src/UnitTests/RootFindingTests/CubicTest.cs +++ b/src/UnitTests/RootFindingTests/CubicTest.cs @@ -1,5 +1,4 @@ -using System; -using MathNet.Numerics.RootFinding; +using MathNet.Numerics.RootFinding; using NUnit.Framework; namespace MathNet.Numerics.UnitTests.RootFindingTests @@ -13,8 +12,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests var a0 = 1d; var a1 = 5d; var a2 = 2d; - - var roots = Cubic.RealRoots(a2, a1, a0); + var roots = Cubic.RealRoots(a0, a1, a2); var root = roots.Item1; var funcValue = root * (root * (root + a2) + a1) + a0; Assert.That(funcValue, Is.EqualTo(0).Within(1e-14)); @@ -28,7 +26,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests var a0 = -2d; var a1 = 5d; var a2 = -4d; - var roots = Cubic.RealRoots(a2, a1, a0); + var roots = Cubic.RealRoots(a0, a1, a2); Assert.That(roots.Item1, Is.EqualTo(2).Within(1e-14)); Assert.That(roots.Item2, Is.EqualTo(1).Within(1e-14)); Assert.AreEqual(double.NaN, roots.Item3); @@ -40,7 +38,7 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests var a0 = -4d; var a1 = 8d; var a2 = -5d; - var roots = Cubic.RealRoots(a2, a1, a0); + var roots = Cubic.RealRoots(a0, a1, a2); Assert.That(roots.Item1, Is.EqualTo(1).Within(1e-14)); Assert.That(roots.Item2, Is.EqualTo(2).Within(1e-14)); Assert.AreEqual(double.NaN, roots.Item3); @@ -52,10 +50,39 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests var a0 = 6d; var a1 = -5d; var a2 = -2d; - var roots = Cubic.RealRoots(a2, a1, a0); + var roots = Cubic.RealRoots(a0, a1, a2); Assert.That(roots.Item1, Is.EqualTo(3).Within(1e-14)); Assert.That(roots.Item2, Is.EqualTo(-2).Within(1e-14)); Assert.That(roots.Item3, Is.EqualTo(1).Within(1e-14)); } + + [TestCase(6.0, -5.0, -2.0, 1.0, -2.0, 3.0, 1.0)] + public void ComplexRoots_TripleReal(double d, double c, double b, double a, double x1, double x2, double x3) + { + var roots = FindRoots.Cubic(d, c, b, a); + Assert.That(roots.Item1.Real, Is.EqualTo(x1).Within(1e-14)); + Assert.That(roots.Item1.Imaginary, Is.EqualTo(0).Within(1e-14)); + Assert.That(roots.Item2.Real, Is.EqualTo(x2).Within(1e-14)); + Assert.That(roots.Item2.Imaginary, Is.EqualTo(0).Within(1e-14)); + Assert.That(roots.Item3.Real, Is.EqualTo(x3).Within(1e-14)); + Assert.That(roots.Item3.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)] + [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(double d, double c, double b, double a) + { + var roots = FindRoots.Cubic(d, c, b, a); + Assert.That(Evaluate.Polynomial(roots.Item1, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); + Assert.That(Evaluate.Polynomial(roots.Item1, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); + Assert.That(Evaluate.Polynomial(roots.Item2, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); + Assert.That(Evaluate.Polynomial(roots.Item2, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); + Assert.That(Evaluate.Polynomial(roots.Item3, d, c, b, a).Real, Is.EqualTo(0).Within(1e-12)); + Assert.That(Evaluate.Polynomial(roots.Item3, d, c, b, a).Imaginary, Is.EqualTo(0).Within(1e-12)); + } } }