This thread is a continuation of my last asked question.
My friends and I are trying to make a multi-functional calculator that can automatically calculate
the roots of quadratic, cubic, and even quartic equations.
However, we have run into some problems, previously on cubic equations, now on quartic equations.
Unlike the last thread, we will provide images including the answers calculated by our program and the answers calculated by WolframAlpha instead of putting them into text here.
The quartic equation will be denoted as in
ax^4 bx^3 cx^2 dx e
Our method of finding the roots of a quartic equation. (Found on another thread from another website)
a=1 b=2 c=-41 d=-42 e=360
a=1 b=-0.25 c=-0.85 d=1.45 e=-4.35
a=1 b=2 c=3 d=4 e=5
As seen in the images, the discriminant calculated by the program is always correct. However, the roots of the equation can sometimes cause errors (NaN), and the nature of roots can also cause errors as well but so as often as the roots do.
With all the information provided above, we wish to know the answers to the following questions :
- What exactly is causing these errors and problems,
- What amendment should I make to my code so that the program can function properly and as intended,
- Whether our method of finding the roots of a quartic equation is correct, or is it even possible to find the roots of all quartic equations by a general formula,
- Last but not least, is there a way to improve, or condense the code lines in our program but still function as we intended it to.
As aforementioned in the previous thread, we have just started learning C#.
For your convenience, we have put notes to some code lines.
Codes (C#) :
using System;
using System.Numerics;
namespace ComplexNumbers
{
public class ComplexNo
{
public static void Main()
{
// prompting for coefficients and constant
Console.WriteLine("Quartic 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 fourth coefficient : ");
double d = Convert.ToDouble(Console.ReadLine());
Console.Write("Insert constant : ");
double e = Convert.ToDouble(Console.ReadLine());
// discriminant
double discriminant = 256 * Math.Pow(a, 3) * Math.Pow(e, 3) - (192 * Math.Pow(a, 2) * b * d * Math.Pow(e, 2)) -
(128 * Math.Pow(a, 2) * Math.Pow(c, 2) * Math.Pow(e, 2)) (144 * Math.Pow(a, 2) * c * Math.Pow(d, 2) * e) -
(27 * Math.Pow(a, 2) * Math.Pow(d, 4)) (144 * a * Math.Pow(b, 2) * c * Math.Pow(e, 2)) -
(6 * a * Math.Pow(b, 2) * Math.Pow(d, 2) * e) - (80 * a * b * Math.Pow(c, 2) * d * e)
(18 * a * b * c * Math.Pow(d, 3)) (16 * a * Math.Pow(c, 4) * e) - (4 * a * Math.Pow(c, 3) * Math.Pow(d, 2)) -
(27 * Math.Pow(b, 4) * Math.Pow(e, 2)) (18 * Math.Pow(b, 3) * c * d * e) -
(4 * Math.Pow(b, 3) * Math.Pow(d, 3)) - (4 * Math.Pow(b, 2) * Math.Pow(c, 3) * e)
(Math.Pow(b, 2) * Math.Pow(c, 2) * Math.Pow(d, 2));
// different parts of the roots presented in p1-p6
double p1 = (2 * Math.Pow(c, 3)) - (9 * b * c * d) (27 * a * Math.Pow(d, 2)) (27 * Math.Pow(b, 2) * e) -
(72 * a * c * e);
double p2 = p1 Math.Sqrt(-4.0 * Math.Pow(Math.Pow(c, 2) - (3 * b * d) (12 * a * e), 3) Math.Pow(p1, 2));
double p3 = ((Math.Pow(c, 2) - (3 * b * d) (12 * a * e)) / (3 * a * Math.Cbrt(p2 / 2))) (Math.Cbrt(p2 / 2) / (3 * a));
double p4 = Math.Sqrt((Math.Pow(b, 2) / (4 * Math.Pow(a, 2))) - (2 * c / (3 * a)) p3);
double p5 = (Math.Pow(b, 2) / (2 * Math.Pow(a, 2))) - (4 * c / (3 * a)) - p3;
double p6 = ((-(Math.Pow(b, 3) / Math.Pow(a, 3))) (4 * b * c / Math.Pow(a, 2)) - (8 * d / a)) / (4 * p4);
// calculating roots with p1-p6
Complex x1 = (-(b / (4 * a))) - (p4 / 2) - (Math.Sqrt(p5 - p6) / 2);
Complex x2 = (-(b / (4 * a))) - (p4 / 2) (Math.Sqrt(p5 - p6) / 2);
Complex x3 = (-(b / (4 * a))) (p4 / 2) - (Math.Sqrt(p5 - p6) / 2);
Complex x4 = (-(b / (4 * a))) (p4 / 2) (Math.Sqrt(p5 - p6) / 2);
// polynomials used for testing the nature of roots
double P = (8 * a * c) - (3 * Math.Pow(b, 2));
double R = Math.Pow(b, 3) (8 * d * Math.Pow(a, 2)) - (4 * a * b * c);
double deltazero = Math.Pow(c, 2) - (3 * b * d) (12 * a * e);
double D = (64 * Math.Pow(a, 3) * e) - (16 * Math.Pow(a, 2) * Math.Pow(c, 2))
(16 * a * Math.Pow(b, 2) * c) - (16 * Math.Pow(a, 2) * b * d) - (3 * Math.Pow(b, 4));
// results
Console.WriteLine(" ");
Console.WriteLine("Solve for " a "x^4 " " " b "x^3 " " " c "x^2 " " " d "x " " " e " :");
Console.WriteLine(" ");
Console.WriteLine("Results :");
Console.WriteLine("First Root :");
string x1rp = Convert.ToString(Math.Round(x1.Real, 5)); // rounding up both real and unreal parts
string x1ip = Convert.ToString(Math.Round(x1.Imaginary, 5));
if (x1.Imaginary == 0)
{
Console.WriteLine(x1rp);
}
else if (x1.Imaginary < 0) // testing if unreal part is negative
{
Console.WriteLine(x1rp " - " Math.Round(Complex.Abs(x1.Imaginary)) "i");
}
else
{
Console.WriteLine(x1rp " " x1ip "i");
}
string x2rp = Convert.ToString(Math.Round(x2.Real, 5)); // rounding up both real and unreal parts
string x2ip = Convert.ToString(Math.Round(x2.Imaginary, 5));
if (x2.Imaginary == 0)
{
Console.WriteLine(x2rp);
}
else if (x2.Imaginary < 0) // testing if unreal part is negative
{
Console.WriteLine(x2rp " - " Math.Round(Complex.Abs(x2.Imaginary)) "i");
}
else
{
Console.WriteLine(x2rp " " x2ip "i");
}
string x3rp = Convert.ToString(Math.Round(x3.Real, 5)); // rounding up both real and unreal parts
string x3ip = Convert.ToString(Math.Round(x3.Imaginary, 5));
if (x3.Imaginary == 0)
{
Console.WriteLine(x3rp);
}
else if (x3.Imaginary < 0) // testing if unreal part is negative
{
Console.WriteLine(x3rp " - " Math.Round(Complex.Abs(x3.Imaginary)) "i");
}
else
{
Console.WriteLine(x3rp " " x3ip "i");
}
string x4rp = Convert.ToString(Math.Round(x4.Real, 5)); // rounding up both real and unreal parts
string x4ip = Convert.ToString(Math.Round(x4.Imaginary, 5));
if (x4.Imaginary == 0)
{
Console.WriteLine(x4rp);
}
else if (x4.Imaginary < 0) // testing if unreal part is negative
{
Console.WriteLine(x4rp " - " Math.Round(Complex.Abs(x4.Imaginary)) "i");
}
else
{
Console.WriteLine(x4rp " " x4ip "i");
}
// discriminant tests
if (discriminant == 0) // if delta = 0
{
if (P < 0 && D < 0 && deltazero != 0)
{
string result4 = "The equation yields a pair of double roots and two real simple roots.";
Console.WriteLine(result4);
}
else if (D > 0 || (P > 0 && (D != 0 || R != 0)))
{
string result5 = "The equation yields a pair of double roots and two complex conjugate roots.";
Console.WriteLine(result5);
}
else if (deltazero == 0 && D != 0)
{
string result6 = "The equation yields a triple root and a simple root.";
Console.WriteLine(result6);
}
else if (D == 0)
{
if (P < 0)
{
string result7 = "The equation yields two pairs of real double roots.";
Console.WriteLine(result7);
}
else if (P > 0 && R == 0)
{
string result8 = "The equation yields two pairs of double complex conjugate roots.";
Console.WriteLine(result8);
}
else if (deltazero == 0)
{
string result9pt2 = Convert.ToString(-(b / (4 * a)));
string result9 = "All of the roots of this equation are equal to -b/4a";
Console.WriteLine(result9 "(" result9pt2 ")");
}
else
{
Console.WriteLine("ERROR");
}
}
else
{
Console.WriteLine("ERROR");
}
}
else if (discriminant > 0) // if delta > 0
{
if (P < 0 && D < 0)
{
string result2 = "The equation yields four distinct real roots.";
Console.WriteLine(result2);
}
else if (P > 0 && D > 0)
{
string result3 = "The equation yields two pairs of unreal complex conjugate roots.";
Console.WriteLine(result3);
}
else
{
Console.WriteLine("ERROR");
}
}
else if (discriminant < 0) // if delta < 0
{
string result1 = "The equation yields two distinct real roots and two complex conjugate non-real roots.";
Console.WriteLine(result1);
}
// discriminant of equation
Console.WriteLine(" ");
Console.WriteLine("Discriminant :");
Console.WriteLine(discriminant);
Console.ReadLine();
}
}
}
==========================================================
Edit #1 :
I tried another method to solve the roots.
However, the imaginary part is still wrong for the last 2 roots.
If anyone is able to spot my mistake, I would be grateful.
// different parts of the roots presented in p1-p6
double p = ((8 * a * c) - (3 * Math.Pow(b, 2))) / (8 * Math.Pow(a, 2));
double q = (Math.Pow(b, 3) - (4 * a * b * c) (8 * Math.Pow(a, 2) * d)) / (8 * Math.Pow(a, 3));
double D0 = Math.Pow(c, 2) - (3 * b * d) (12 * a * e);
double D1 = (2 * Math.Pow(c, 3)) - (9 * b * c * d) (27 * Math.Pow(b, 2) * e) (27 * a * Math.Pow(d, 2)) - (72 * a * c * e);
Complex Q = Complex.Pow((D1 Complex.Sqrt(Complex.Pow(D1, 2) - (4 * Complex.Pow(D0, 3)))) / 2.0, 1.0 / 3.0);
Complex S = 1.0 / 2.0 * Complex.Sqrt((-2.0 / 3.0 * p) (1.0 / (3 * a) * (Q (D0 / Q))));
// calculating roots with p1-p6
Complex x1 = (-(b / 4 * a)) - S (1.0 / 2.0 * Complex.Sqrt((-4.0 * Complex.Pow(S, 2)) - (2 * p) (q / S)));
Complex x2 = (-(b / 4 * a)) - S - (1.0 / 2.0 * Complex.Sqrt((-4.0 * Complex.Pow(S, 2)) - (2 * p) (q / S)));
Complex x3 = (-(b / 4 * a)) S (1.0 / 2.0 * Complex.Sqrt((-4.0 * Complex.Pow(S, 2)) - (2 * p) (q / S)));
Complex x4 = (-(b / 4 * a)) S - (1.0 / 2.0 * Complex.Sqrt((-4.0 * Complex.Pow(S, 2)) - (2 * p) (q / S)));
==========================================================
Edit #2
Problem solved.
The program is now able to solve quartic equations. The cause of the problems mentioned above is actually the minor typos I made in the code lines.
CodePudding user response：
I got an algorithm from here that works. It seems simpler than yours.
double D0 = c * c - 3 * b * d 12 * a * e;
double D1 = 2 * c * c * c - 9 * b * c * d 27 * b * b * e 27 * a * d * d - 72 * a * c * e;
double p = (8 * a * c - 3 * b * b) / (8 * a * a);
double q = (b * b * b - 4 * a * b * c 8 * a * a * d) / (8 * a * a * a);
Complex Q = Complex.Pow((D1 Complex.Sqrt(D1 * D1 - 4 * D0 * D0 * D0)) / 2, 1.0 / 3.0);
Complex S = Complex.Sqrt(-2 * p / 3 (Q D0 / Q) / (3 * a)) / 2;
Complex u = Complex.Sqrt(-4 * S * S - 2 * p q / S) / 2;
Complex v = Complex.Sqrt(-4 * S * S - 2 * p - q / S) / 2;
Complex x1 = -b / (4 * a) - S u;
Complex x2 = -b / (4 * a) - S - u;
Complex x3 = -b / (4 * a) S v;
Complex x4 = -b / (4 * a) S - v;
Also, according to this page, you can use this instead of the very long discriminant formula just to get its sign:
double discriminantSign = c * c - 3 * b * d 12 * a * e;
I also noticed 4 typos, you forgot to mention the 5 decimal when you apply the Math.Round function, so your results are rounded to an integer.
EDIT: If you still have issues in you imaginary part, it's because you didn't correct the typos:
Console.WriteLine(x1rp " - " Math.Round(Complex.Abs(x1.Imaginary)) "i");
Into:
Console.WriteLine(x1rp " - " Math.Round(Complex.Abs(x1.Imaginary), 5) "i");
You take 5 decimals everywhere except in 4 lines as above. So it's rounding the integer. Replace them all.