Click here to Skip to main content
14,923,504 members
Articles / Programming Languages / C#
Article
Posted 18 Jul 2014

Stats

35.8K views
906 downloads
10 bookmarked

To Solve a Cubic Equation

Rate me:
Please Sign up or sign in to vote.
4.76/5 (8 votes)
18 Jul 2014CPOL3 min read
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 :)

The Wikipedia description starts with the qubic equation

Image 1

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

Image 2   and   Image 3

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

Image 4   and   Image 5

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:

Image 6   and   Image 7

For the complex angles we get:

Image 8   and   Image 9

And:

Image 10   and   Image 11

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:

Image 12

If we separate:

Image 13

And:

Image 14

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 :)

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

Image 15

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.

Image 16   and   Image 17

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:

Image 18

The angle becomes for u:

Image 19

The angle becomes for v:

Image 20

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:

Image 21

Image 22

Image 23

And as αv = 2π - α

Image 24

Image 25

Image 26

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

Image 27

Image 28

Image 29

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 2 *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 :)

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
Tester / Quality Assurance Annax Switzerland AG
Switzerland Switzerland
Computers are very straight... They always do exactly what we tell them to do... Only, much too often what we tell them to do is not really what we want them to do Smile | :)

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 with it. Besides doing software for HMI's on C# for business, I enjoy very much to implement interesting algorithms and analyse the mathematics they are based on in my leisure time Smile | :)

For more detailed descriptions and math visit me on my own page

www.mosismath.com

Comments and Discussions

 
QuestionA good approach Pin
Harvey Triana30-Jul-20 20:15
MemberHarvey Triana30-Jul-20 20:15 
Suggestionmake a easier program(I hope) to solve qubic equation Pin
Member 1460208024-Sep-19 19:44
MemberMember 1460208024-Sep-19 19:44 
QuestionBug Pin
marygold893-Jan-16 20:28
Membermarygold893-Jan-16 20:28 
AnswerRe: Bug Pin
Mosi_6215-Mar-16 9:52
professionalMosi_6215-Mar-16 9:52 
Questionhow to download the source code Pin
marygold893-Jan-16 6:35
Membermarygold893-Jan-16 6:35 
AnswerRe: how to download the source code Pin
Mosi_6214-Mar-16 9:16
professionalMosi_6214-Mar-16 9:16 
QuestionBug Pin
MartinBell15-Dec-15 6:34
MemberMartinBell15-Dec-15 6:34 
AnswerRe: Bug Pin
Mosi_6214-Mar-16 9:37
professionalMosi_6214-Mar-16 9:37 
GeneralRe: Bug Pin
Member 1285975619-Nov-16 10:53
MemberMember 1285975619-Nov-16 10:53 
QuestionDiscussion of the cubic roots Pin
YvesDaoust21-Jul-14 3:04
MemberYvesDaoust21-Jul-14 3:04 
AnswerRe: Discussion of the cubic roots Pin
Mosi_6222-Jul-14 8:35
professionalMosi_6222-Jul-14 8:35 
QuestionI'll take your word for it! Pin
.dan.g.20-Jul-14 20:15
professional.dan.g.20-Jul-14 20:15 
QuestionVote of 5 Pin
Kenneth Haugland18-Jul-14 18:04
professionalKenneth Haugland18-Jul-14 18:04 
AnswerRe: Vote of 5 Pin
Mosi_6219-Jul-14 5:40
professionalMosi_6219-Jul-14 5:40 

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

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