Click here to Skip to main content
Click here to Skip to main content

Tagged as

To Solve a Cubic Equation

, 18 Jul 2014 CPOL
Rate this:
Please Sign up or sign in to vote.
Gerolamo Cardano published a method to solve a cubic equation in 1545. There is a description of this method on Wikipedia. But it is not too detailed and on the German Wikipedia. In fact, the last part is missing and without this part, one cannot implement it into an algorithm. So I thought I could

Cardano's Method

Gerolamo Cardano published a method to solve a cubic equation in 1545. There is a description of this method on Wikipedia. But it is not too detailed and on the German Wikipedia. In fact, the last part is missing and without this part, one cannot implement it into an algorithm. So I thought I could try to pick up there where the Wikipedia description ends Smile | :)

The Wikipedia description starts with the qubic equation

And even though some details are missing, the Wikipedia description is OK until the part:

  and  

Now things become interesting and most descriptions get vague.

There are the three possibilities: D >= 0 , D = 0 and D < 0.

D > 0

If D > 0 the root of D is real and we get one real solution and two complex solutions for u as well as for v:

The real solutions are

  and  

The two complex solutions are of the same length as u0 and v0, each one has the complex angle 120° for one and 240° for the other solution. That means with:

  and  

For the complex angles we get:

  and  

And:

  and  

So we get three possible solutions for u and three possible solutions for v and if we combine them, we get nine possible solutions. That’s pretty much it, but not all of them are valid. Furthermore, above we had the clause u*v = - p/3 and this is not the case for all nine combinations of u and v. It’s valid only for the combinations uo * vo, u1 * v2 and u2 * v1. Let’s have look at u1 * v2. That’s:

If we separate:

And:

That only works for ε1 * ε2. If we combine ε1 * ε1, there remains an imaginary part and that can’t give - p/3. This is a part that is most often missing in the description of Cardano’s formula:

So, given that D > 0 there are the three solutions. But only the first one is real.

y1 = u0 + v0

y2 = u1 + v2

y3 = u2 + v1

And on the way to this the substitution:

x = y – a/3 = u + v – a/3

has been made. This should not be forgotten as we are looking for "x" finally Smile | :)

If we put this into code it looks like this:

if (d > 0)
{
    u = Xroot(-q / 2.0 + Math.Sqrt(d), 3.0);
    v = Xroot(-q / 2.0 - Math.Sqrt(d), 3.0);
    x1.real = u + v - a / 3.0;
    x2.real = -(u + v) / 2.0 - a / 3.0;
    x2.imag = Math.Sqrt(3.0) / 2.0 * (u - v);
    x3.real = x2.real;
    x3.imag = -x2.imag;
}

D = 0

If D = 0, we have the same three solutions, but as D = 0, u0 = v0 and

And y3 = y2. That means we get only two solutions.

In code that’s:

if (d == 0)
{
    u = Xroot(-q / 2.0, 3.0);
    v = Xroot(-q / 2.0, 3.0);
    x1.real = u + v - a / 3.0;
    x2.real = (u + v) / 2.0 - a / 3.0;
}

D < 0

The case D < 0 is the so called "Casus irreducibilis." Now the root of D becomes complex and the cubic roots for u and v become complex too.

  and  

To solve this complex root we interpret the complex number as number of the length r and complex angle α.

With the substitution from further above:

4((q/2)2 + (p/3)3) = 4D

We get the same length for u and v:

The angle becomes for u:

The angle becomes for v:

To get the cubic root out of this now we have to calculate the cubic root of r and basically divide α by 3 that’s it. At least for the one first solution. The cubic root of a complex number has not just one but three solutions. The angle α we get once by multiplying α/3 by 3, twice by multiplying (2π + α)/3 by 3 and thirdly by multiplying (4π + α)/3 by 3. That makes three solutions for u and three solutions for v. They look like:

And as αv = 2π - α

That would give nine possible solutions again. But as we are interested in real solutions, we can reduce that to three real solutions. For the three combinations

u0 + v2 u1 + v1 u2 + v0

the imaginary part becomes 0 and the solutions are

This is the most complicate part:

if (d < 0)
{
    r = Math.Sqrt(-p * p * p / 27.0);
    alpha = Math.Atan(Math.Sqrt(-d) / -q * 2.0);
    if (q > 0)                         // if q > 0 the angle becomes PI - alpha
        alpha = 2.0 * Math.PI - alpha;

    x1.real = Xroot(r, 3.0) * (Math.Cos((6.0 * Math.PI - alpha) / 3.0)
	    + Math.Cos(alpha / 3.0)) - a / 3.0;
    x2.real = Xroot(r, 3.0) * (Math.Cos((2.0 * Math.PI + alpha) / 3.0)
		+ Math.Cos((4.0 * Math.PI - alpha) / 3.0)) - a / 3.0;
    x3.real = Xroot(r, 3.0) * (Math.Cos((4.0 * Math.PI + alpha) / 3.0)
		+ Math.Cos((2.0 * Math.PI - alpha) / 3.0)) - a / 3.0;
}

But that’s all it takes Smile | :)

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Mosi_62
Software Developer (Senior) Aebi Maschinenfabrik, Burgdorf
Switzerland Switzerland
Writing Software is one of the most creative tings one can do. I have been doing this for more than ten years now and still having a lot of fun whit it. Besides doing software for farming Tractors for business, I enjoy very much to implement interesting algorithms and analyse the mathematics they are based on in my leisure time Smile | :)

Comments and Discussions

 
QuestionDiscussion of the cubic roots [modified] PinmemberYvesDaoust21-Jul-14 3:04 
AnswerRe: Discussion of the cubic roots PinprofessionalMosi_6222-Jul-14 8:35 
QuestionI'll take your word for it! Pinmvp.dan.g.20-Jul-14 20:15 
QuestionVote of 5 PinprofessionalKenneth Haugland18-Jul-14 18:04 
AnswerRe: Vote of 5 PinprofessionalMosi_6219-Jul-14 5:40 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web01 | 2.8.141022.2 | Last Updated 18 Jul 2014
Article Copyright 2014 by Mosi_62
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid