I am currently trying to solve a quartic equation using Ferrari's method from Wikipedia. I want to retrieve only the real roots, discarding the imaginary one. My implementation does not return the good value for real roots. I can't find the mistakes in the formula.
My CubicEquation works as expected, and my bi-quadratic solving either. Now, I am only missing the Ferrari's method to be done but I can't make it work!
Here's my class:
public class QuarticFunction {
private final double a;
private final double b;
private final double c;
private final double d;
private final double e;
public QuarticFunction(double a, double b, double c, double d, double e) {
this.a = a;
this.b = b;
this.c = c;
this.d = d;
this.e = e;
}
public final double[] findRealRoots() {
if (a == 0) {
return new CubicEquation(b, c, d, e).findRealRoots();
}
if (isBiquadratic()) { //This part works as expected
return solveUsingBiquadraticMethod();
}
return solveUsingFerrariMethodWikipedia();
}
private double[] solveUsingFerrariMethodWikipedia() {
// http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution
// ERROR: Wrong numbers when Two Real Roots + Two Complex Roots
// ERROR: Wrong numbers when Four Real Roots
QuarticFunction depressedQuartic = toDepressed();
if (depressedQuartic.isBiquadratic()) {
return depressedQuartic.solveUsingBiquadraticMethod();
}
double y = findFerraryY(depressedQuartic);
double originalRootConversionPart = -b / (4 * a);
double firstPart = Math.sqrt(depressedQuartic.c + 2 * y);
double positiveSecondPart = Math.sqrt(-(3 * depressedQuartic.c + 2 * y + 2 * depressedQuartic.d / Math.sqrt(depressedQuartic.c + 2 * y)));
double negativeSecondPart = Math.sqrt(-(3 * depressedQuartic.c + 2 * y - 2 * depressedQuartic.d / Math.sqrt(depressedQuartic.c + 2 * y)));
double x1 = originalRootConversionPart + (firstPart + positiveSecondPart) / 2;
double x2 = originalRootConversionPart + (-firstPart + negativeSecondPart) / 2;
double x3 = originalRootConversionPart + (firstPart - positiveSecondPart) / 2;
double x4 = originalRootConversionPart + (-firstPart - negativeSecondPart) / 2;
Set<Double> realRoots = findAllRealRoots(x1, x2, x3, x4);
return toDoubleArray(realRoots);
}
private double findFerraryY(QuarticFunction depressedQuartic) {
double a3 = 1;
double a2 = 5 / 2 * depressedQuartic.c;
double a1 = 2 * Math.pow(depressedQuartic.c, 2) - depressedQuartic.e;
double a0 = Math.pow(depressedQuartic.c, 3) / 2 - depressedQuartic.c * depressedQuartic.e / 2
- Math.pow(depressedQuartic.d, 2) / 8;
//CubicEquation works as expected! No need to worry! It returns either 1 or 3 roots.
CubicEquation cubicEquation = new CubicEquation(a3, a2, a1, a0);
double[] roots = cubicEquation.findRealRoots();
for (double y : roots) {
if (depressedQuartic.c + 2 * y != 0) {
return y;
}
}
throw new IllegalStateException("Ferrari method should have at least one y");
}
public final boolean isBiquadratic() {
return Double.compare(b, 0) == 0 && Double.compare(d, 0) == 0;
}
private double[] solveUsingBiquadraticMethod() {
//It works as expected!
QuadraticLine quadraticEquation = new QuadraticLine(a, c, e);
if (!quadraticEquation.hasRoots()) {
return new double[] {};
}
double[] quadraticRoots = quadraticEquation.findRoots();
Set<Double> roots = new HashSet<>();
for (double quadraticRoot : quadraticRoots) {
if (quadraticRoot > 0) {
roots.add(Math.sqrt(quadraticRoot));
roots.add(-Math.sqrt(quadraticRoot));
} else if (quadraticRoot == 0.00) {
roots.add(0.00);
}
}
return toDoubleArray(roots);
}
public QuarticFunction toDepressed() {
// http://en.wikipedia.org/wiki/Quartic_function#Converting_to_a_depressed_quartic
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 * d * Math.pow(a, 2)) / (8 * Math.pow(a, 3));
double r = (-3 * Math.pow(b, 4) + 256 * e * Math.pow(a, 3) - 64 * d * b * Math.pow(a, 2) + 16 * c * a
* Math.pow(b, 2)) / (256 * Math.pow(a, 4));
return new QuarticFunction(1, 0, p, q, r);
}
private Set<Double> findAllRealRoots(double... roots) {
Set<Double> realRoots = new HashSet<>();
for (double root : roots) {
if (!Double.isNaN(root)) {
realRoots.add(root);
}
}
return realRoots;
}
private double[] toDoubleArray(Collection<Double> values) {
double[] doubleArray = new double[values.size()];
int i = 0;
for (double value : values) {
doubleArray[i] = value;
++i;
}
return doubleArray;
}
}
I have tried a simpler implementation found here, but I had a new problem. Now, half of the roots are good, but the other half is wrong. Here is my tentative:
private double[] solveUsingFerrariMethodTheorem() {
// https://proofwiki.org/wiki/Ferrari's_Method
CubicEquation cubicEquation = getCubicEquationForFerrariMethodTheorem();
double[] cubicRoots = cubicEquation.findRealRoots();
double y1 = findFirstNonZero(cubicRoots);
double inRootOfP = Math.pow(b, 2) / Math.pow(a, 2) - 4 * c / a + 4 * y1;
double p1 = b / a + Math.sqrt(inRootOfP);
double p2 = b / a - Math.sqrt(inRootOfP);
double inRootOfQ = Math.pow(y1, 2) - 4 * e / a;
double q1 = y1 - Math.sqrt(inRootOfQ);
double q2 = y1 + Math.sqrt(inRootOfQ);
double x1 = (-p1 + Math.sqrt(Math.pow(p1, 2) - 8 * q1)) / 4;
double x2 = (-p2 - Math.sqrt(Math.pow(p2, 2) - 8 * q1)) / 4;
double x3 = (-p1 + Math.sqrt(Math.pow(p1, 2) - 8 * q2)) / 4;
double x4 = (-p2 - Math.sqrt(Math.pow(p2, 2) - 8 * q2)) / 4;
Set<Double> realRoots = findAllRealRoots(x1, x2, x3, x4);
return toDoubleArray(realRoots);
}
private CubicEquation getCubicEquationForFerrariMethodTheorem() {
double cubicA = 1;
double cubicB = -c / a;
double cubicC = b * d / Math.pow(a, 2) - 4 * e / a;
double cubicD = 4 * c * e / Math.pow(a, 2) - Math.pow(b, 2) * e / Math.pow(a, 3) - Math.pow(d, 2)
/ Math.pow(a, 2);
return new CubicEquation(cubicA, cubicB, cubicC, cubicD);
}
private double findFirstNonZero(double[] values) {
for (double value : values) {
if (Double.compare(value, 0) != 0) {
return value;
}
}
throw new IllegalArgumentException(values + " does not contain any non-zero value");
}
I don't know what I am missing. I have spent hours debugging trying to see the errors, but now I'm completely lost (plus some headaches). What am I doing wrong? I don't care which formulas to use, since I only want one of them to work.
I had two mistakes.
Here is my resulting implementation