Having trouble in finding roots of a cubic equation with C#-CodePudding

Recently my friends and I decided to make a multi-functional calculator to automatically find the roots of quadratic equations. It wasn't that challenging after all, so we decided to go onto the next level, to make a calculator for a cubic equation.

(ax^3 bx^2 cx d)

However, we stumbled across some trivial problems, and no matter how hard we tried, the results are still the same. we are beginners in terms of coding, so we're not sure if we are actually making some stupid mistakes here but at least we want to learn something from others.

Basically, we have tried a lot of combinations of different cubic equations, or even re-coding the whole thing. The problem is that the results we yield are always wrong, but only for the real parts of the second root and third root.

For a better understanding, we have tried 9x^3 8x^2 7x 6, as an example.
The correct answer, according to a cubic equation calculator website, is

(The website)

First root = -0.87285
Second root = -0.00802 0.87391 i
Third root = -0.00802 - 0.87391 i

However, our result to this equation is :

First root = -0.87285
Second root = -0.2963 0.87391 i
Third root = -0.2963 -0.87391 i

It is apparently noticeable that only parts of the second and third roots are wrong.
We have tried finding similar threads to help, but those are a little bit too difficult for us to understand or those problems are not the same as ours.
We look forward to finding the solution and the causes to this problem.

We have separated the formulas for finding the roots of a cubic equation into 5 parts. (rootp1-rootp5)
The formulas are coded according to the formulas that can be found in the Chinese version of Wikipedia's page of Cubic equation.
(The formulas)
We have also rounded the real parts of the roots to 5 digits.

(Some code lines might be redundant, but as I have mentioned, we are new to coding)

Code (C#) :

using System;
using System.Numerics;

namespace ComplexNumbers
{
    public class ComplexNo
    {
        public static void Main()
        {
            Console.WriteLine("Cubic Equation Calculator :");
            Console.Write("Insert first coefficient : ");
            double a = Convert.ToDouble(Console.ReadLine());
            Console.Write("Insert second coefficient : ");
            double b = Convert.ToDouble(Console.ReadLine());
            Console.Write("Insert third coefficient : ");
            double c = Convert.ToDouble(Console.ReadLine());
            Console.Write("Insert constant : ");
            double d = Convert.ToDouble(Console.ReadLine());

                Console.WriteLine(" ");
                Console.WriteLine("Solve for "   a   "x"   "^3 "   b   "x^2 "   c   "x "   d   " :");
                double rootp1 = -(b / (3 * a));
                double rootp2 = (b * c / (6 * Math.Pow(a, 2))) - (Math.Pow(b, 3) / (27 * Math.Pow(a, 3))) - (d / (2 * a));
                double rootp3 = (c / (3 * a)) - (Math.Pow(b, 2) / (9 * Math.Pow(a, 2)));
                Complex root1 = rootp1   Math.Cbrt(rootp2   Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3)))  
                                 Math.Cbrt(rootp2 - Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3)));
                Complex rootp4 = new Complex(-1 / 2, Math.Sqrt(3) / 2);
                Complex rootp5 = new Complex(-1 / 2, -(Math.Sqrt(3) / 2));
                Complex root2 = rootp1   (rootp4 * Math.Cbrt(rootp2   Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3))))  
                                 (rootp5 * Math.Cbrt(rootp2 - Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3))));
                Complex root3 = rootp1   (rootp5 * Math.Cbrt(rootp2   Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3))))  
                                 (rootp4 * Math.Cbrt(rootp2 - Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3))));

                Console.WriteLine(" ");
                Console.WriteLine("Results :");
                Console.WriteLine("First Root :");
                string root1rp = Convert.ToString(Math.Round(root1.Real, 5));
                string root1ip = Convert.ToString(Math.Round(root1.Imaginary, 5));
                Console.WriteLine(root1rp   "   "   root1ip   "i");
                Console.WriteLine("Second Root :");
                string root2rp = Convert.ToString(Math.Round(root2.Real, 5));
                string root2ip = Convert.ToString(Math.Round(root2.Imaginary, 5));
                Console.WriteLine(root2rp   "   "   root2ip   "i");
                Console.WriteLine("Third Root :");
                string root3rp = Convert.ToString(Math.Round(root3.Real, 5));
                string root3ip = Convert.ToString(Math.Round(root3.Imaginary, 5));
                Console.WriteLine(root3rp   "   "   root3ip   "i");

                Console.ReadLine();
        }
    }
}

(Sorry for making this thread so long and my bad grammar)

CodePudding user response：

The problem is with this line Complex rootp4 = new Complex(-1 / 2, Math.Sqrt(3) / 2);. -1/2 uses integer division and evaluates to 0.

This code will work.

        Console.WriteLine(" ");
        Console.WriteLine("Solve for "   a   "x"   "^3 "   b   "x^2 "   c   "x "   d   " :");
        double minusBover3a = -(b / (3 * a));
        double rootp2 = (b * c / (6 * Math.Pow(a, 2))) - (Math.Pow(b, 3) / (27 * Math.Pow(a, 3))) - (d / (2 * a));
        double rootp3 = (c / (3 * a)) - (Math.Pow(b, 2) / (9 * Math.Pow(a, 2)));

        double bigCubeRootPlus = Math.Cbrt(rootp2   Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3)));
        double bigCubeRootMinus = Math.Cbrt(rootp2 - Math.Sqrt(Math.Pow(rootp2, 2)   Math.Pow(rootp3, 3)));
        // ** THIS IS THE PROBLEM.  "-1/2" uses integer division, so this complex has 0 for real part
        Complex complexPlus = new Complex(-1.0 / 2, Math.Sqrt(3) / 2);
        Complex complexMinus = new Complex(-1.0 / 2, -(Math.Sqrt(3) / 2));

        Complex root1 = minusBover3a   bigCubeRootPlus   bigCubeRootMinus;
        Complex root2 = minusBover3a   complexPlus * bigCubeRootPlus   complexMinus * bigCubeRootMinus;
        Complex root3 = minusBover3a   complexMinus * bigCubeRootPlus   complexPlus * bigCubeRootMinus;

CodePudding user response：

I admit I haven't tried your code, but you may have a casting issue in the math in your original calculations. Try this:

double rootp2 = ((double)b * (double)c / (6D * Math.Pow(a, 2D))) - (Math.Pow(b, 3D) / (27D * Math.Pow(a, 3D))) - ((double)d / (2D * (double)a));

If that makes a difference, you'd have to propagate a similar change (cast all variables as (double) and inline constants as double with D) through the other calculations.

CodePudding user response：

The following is an alternative to your code:

    //  from https://www.daniweb.com/programming/software-development/
    //  code/454493/solving-the-cubic-equation-using-the-complex-struct

    //  algorithm described in
    //  https://en.wikipedia.org/wiki/Cubic_equation#General_cubic_formula
    const int NRoots = 3;

    double SquareRootof3 = Math.Sqrt(3);
    // the 3 cubic roots of 1
    var CubicUnity = new List<Complex>(NRoots)
                { new Complex(1, 0), 
                  new Complex(-0.5, -SquareRootof3 / 2.0), 
                  new Complex(-0.5, SquareRootof3 / 2.0) };
    // intermediate calculations
    double DELTA = 18 * a * b * c * d 
                 - 4 * b * b * b * d 
                   b * b * c * c 
                 - 4 * a * c * c * c 
                 - 27 * a * a * d * d;
    double DELTA0 = b * b - 3 * a * c;
    double DELTA1 = 2 * b * b * b 
                  - 9 * a * b * c 
                    27 * a * a * d;
    Complex DELTA2 = -27 * a * a * DELTA;
    Complex C = Complex.Pow((DELTA1   Complex.Pow(DELTA2, 0.5)) / 2, 1 / 3.0);

    for (int i = 0; i < NRoots; i  )
    {
        Complex M = CubicUnity[i] * C;
        Complex Root = -1.0 / (3 * a) * (b   M   DELTA0 / M);
        Console.WriteLine();
        Console.WriteLine($"Root {i 1}:");
        Console.WriteLine($"Real      {Root.Real:0.#####}");
        Console.WriteLine($"Imaginary {Root.Imaginary:0.#####}i");
    }