Solving the Cubic Equation
There are explanations online of how to solve a cubic equation that do a great job of introducing the principles of how to solve the cubic formula for certain types of cubic equations. However, I haven't seen one that convincingly explains how to find all the real roots algorithmically for all possible cubic equations. In this article we'll work cooperatively with a small program to discover how to solve the cubic equation for every real-value coefficient and non-zero leading coefficient.
The Depressed Cubic
Let's first start with a cubic polynomial.
The first task is to convert this form into a "depressed cubic". This is a transform of the cubic equation which brings its inflection point, i.e. where the 2nd derivative is zero, to lie on the y-axis. Why this is advantageous will be shown later.
Now that we know the x position of the original cubic equation's inflection point we can re-write the cubic equation of a new x term shifting the whole equation in the x-axis. This will change the values of our roots but since it's a linear transformation of them it will be easy to undo later.
And substituting it into the original cubic equation we get this.
At this point it would appear as if we've just made the equation even harder to solve. However the squared term works out such that it's always equal to zero which can be shown via substitution.
And we can divide out by the leading coefficient in order to push all coefficients into two terms.
Since our new coefficients are quite clunky and we won't need to look at their structure from here on out we can simply rename them. While it is strange that one is a derivative of the other this fact won't be used later.
So up to here we've gone through a number of steps to simply transform the x coordinate system to simplify the equation to this depressed cubic form. While it might seem like that was a lot of work to go from a cubic we can't solve to a slightly simpler cubic that we also can't solve that was the easy part.
Cardano's Formula
These next steps are based off of Gerolamo Cardano's solution based off of the solutions by Scipione del Ferro and Tartaglia. Why this substitution and constraint were chosen and why they work are not terribly straightforward and I don't understand it well myself so we'll take it for granted and move on.
From here we can substitute for the x coordinate and expand.
At this point the equation is as expanded as it can get. The next thing we'll want to do is to eliminate all the instances of one of the substituted variables. Since the first and second substituted variable are each taken to the 3rd power multiplying the whole equation by the cube of the first variable will make it possible to use the constraint to eliminate all instances of the second variable.
And now we're back to a rather small formula. It would appear that we've just made things worse though since we're now dealing with x to the 6th power. However, a small substitution turns this into a quadratic equation.
And from here we can use the quadratic equation to solve.
With this we're as deep into the algebraic solving as is needed. From here things go in reverse.
Writing the first set of substitution variables in terms of the most recent one gives us a cube root. And continuing to substitute into that equation gives us an equation with a cube root containing a square root. We'll want to keep track of which particular root we're receiving later so we can distinguish them because each remaining substitution variable can take on 6 values and together, naively, 36 values which can't all possibly be useful for describing roots.
Here we'll re-label some items so that it's shorter to write out later.
And removing the last of our substitution variables we get an expression for all possible roots though not all of them are valid at all times as we'll see. The root is still a function of the coefficients from earlier polynomial forms but that part is left implicit and the phase and sign part is explicit since it will need to be selected.
Detecting Correct Roots
Note that while the particular phase of the cube root doesn't matter if it is immediately raised to the 3rd power afterwards I'll be using a particular definition that your favorite programming environment may not use and may make the following results unusable. The definition we'll be using keeps the argument of the resulting complex number and for cube roots and between and for square roots. And if you were to open your browser's console on this page you'd be able to see this definition in action.
const u = Mathi.complex(1, 0.001);
Mathi.sqrt(u); // {a: 1.000000124999961, b: 0.0004999999375000274}
Mathi.argz(Mathi.sqrt(u)); // 0.0004999998333334333
const v = Mathi.complex(1, -0.001);
Mathi.sqrt(v); // {a: -1.000000124999961, b: 0.000499999937500203}
Mathi.argz(Mathi.sqrt(v)); // 3.1410926537564596
Whereas if you were to attempt to do the same in Python you might get a result which returns different bounds.
import cmath
u = 1+0.001j
u**(1/2) # (1.000000124999961+0.0004999999375000274j)
cmath.phase(u**(1/2)) # 0.0004999998333334333
v = 1-0.001j
v**(1/2) # (1.000000124999961-0.0004999999375000274j)
cmath.phase(v**(1/2)) # -0.0004999998333334333
Determining which combinations of roots are valid for which combinations of the depressed cubic's coefficients isn't obvious. So from here on we're going to show interactive graphs of depressed cubic functions. You should be able to vary the and coefficients by moving your mouse over the graph area or clicking on the graph area. The coefficient increases by going up and the coefficient increases by going right. You can see the resulting values for and below the graph area. And below those coordinates you'll see three tables quantifying some property which is only true where it's 1.000 and is false for numbers less than 1.000 which are shown in order to demonstrate the continuous and discontinuous approach to 1.000. The Complex Rootness table plugs the given root back into the depressed cubic equation and checks how close it actually is to 0. The Realness table checks how close to 0 the imaginary component of the root is. The Real Rootness table checks that the given root has both Complex Rootness and Realness.
The tables themselves are structured so that takes on values in the 6x6 tables as follows:
There are some helpful, but not necessary, things to note for taking a shortcut to understanding the next steps. Due to how and were defined earlier and need to be different values in order to satisfy Cardano's substitutions. Another useful observation is that and can be swapped and still yield the same value. Swapping and means that roots below the main (top left to bottom right) diagonal in the tables are necessarily duplicates.
The last thing of note before we look at the graphs is that there are two functions drawn in addition to the depressed cubic. The first is and the second is . is important because that's where the slope of the depressed cubic at the inflection point changes sign and is also where the domain of valid roots changes. is important because that's where changes between real and imaginary and is important because that's where the domain of valid roots changes and can also change from 1 real root to 3 real roots.
When interacting with the graph you should note certain behaviors in the Complex Rootness table. There are certain Complex Rootness entries which stay exactly 1.000 over a broad area. There are other entries in the table which approach 1.000 when approaching some 1-dimensional feature in one direction and also gradually fall off continuing through it in that direction. When these 1-dimensional features occur embedded inside a broad area with 3 working complex roots we can ignore them as their roots are special cases of the broader area.
When is greater than zero and not in the bottom left of the grid you should identify these complex roots with the first being real-valued.
When above the line where on the right-hand side of the graph you should identify these roots with the first being real-valued.
Going to its mirror image on the left side of the graph you should identify these roots again with the first being real-valued.
When below the line where you should identify these roots with all of them being real-valued.
With this we've solved the problem almost everywhere. But the "almost" we haven't checked is the seams between the areas we've solved for. Ideally these would smoothly continue the broader areas but that ends up not being the case. For the origin this is easy to check as every form of root gradually approaches 0 even if from different directions. For the seams themselves I'll give you more graphs constrained to these seams in order to check their behavior yourself.
On the side of this constrained graph you'll find a lot of possible complex roots due to the fact that when . And for that same reason you'll see the cubic radical with the added drop out. However, this ultimately yields the exact same degenerate equation for all 3 real values we would end up choosing. So starting again from the depressed cubic equation with the zero terms removed.
The actual values we should be choosing should be as follows with the first root being real-valued.
The underscores indicates that we don't actually care what that value is since it gets cancelled out anyways. So we're technically choosing 9 complex roots which degenerate to 3 distinct values. We could choose the other 9 complex roots but we already know they just swap added terms. The absolute value of is taken so that we can be sure of which roots are real-valued. Ideally we'd be able to show that the inputs to these functions are compatible with those of the areas it joins. However this can't be done with all the functions we've chosen but given both sets of equations it's simple to show that the equations for both areas equal these equations as we approach .
Following the same logic we can get the roots on the mirrored axis. Again the first root is real-valued.
With this we can be confident that these seams don't have issues when it comes to solving for the roots. But one of the remaining seams will cause some issues. The remaining seams are both along and so the function inputs which determine the sign are irrelevant and have been removed from the grid below the graph.
Despite this being a smaller grid any cell below the main diagonal is still a duplicate of anything above it. However, in this case those also hint at double-roots which would appear twice in a product-of-sums form of the polynomial. When simplifying the values for the roots the following equalities might be useful.
With this you should be able to identify these roots with all roots being real-valued.
And likewise for the mirrored case.
Here we run into a problem though. In the case the function for the roots of the areas on both sides of the half-axis became the function for the special case when that case's constraints were imposed. And that still occurs for the case but not the case. The equations are the same coming from the area above but not from the area below. If you go back to the first interactive graph you should be able to notice this discontinuity.
Removing Discontinuities
If you look back at the 1-dimensional root equations you might notice that when you assume a of the opposite sign and you flip the root (or ) value you get the other set of roots. In fact, if you look at the equation for the depressed cubic you can flip it in , flip it in then flip it in and get the same equation.
This means that if there's a situation where we're not confident in the stability of our computation we can opt to solve against this "other" depressed cubic. Flipping isn't particularly relevant to finding roots but flipping negates the sign of our roots and flipping will negate it in our equations but also lets us choose the root functions from the opposite side of the graph. So if we avoid solving directly for then we can remove the discontinuity.
In order to indicate that we're working in this "other" depressed cubic it will be notated with an overline. And since they're both affected by being flipped and get the same treatment.
This next graph doesn't have the typical grid that the previous ones have had.
Instead it chooses the root equations we've found earlier for the cases where
From this it should be possible (with some effort) to determine that you can find the complex roots at any point without a discontinuity and get the exact real roots.
And one last point to remember is that we need to add the