14,923,504 members
Articles / Programming Languages / C#
Article
Posted 18 Jul 2014

35.8K views
10 bookmarked

# To Solve a Cubic Equation

Rate me:
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

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

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

## Share

 Tester / Quality Assurance Annax Switzerland AG 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

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

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

www.mosismath.com

 First Prev Next
 A good approach Harvey Triana30-Jul-20 20:15 Harvey Triana 30-Jul-20 20:15
 make a easier program(I hope) to solve qubic equation Member 1460208024-Sep-19 19:44 Member 14602080 24-Sep-19 19:44
 Bug marygold893-Jan-16 20:28 marygold89 3-Jan-16 20:28
 Re: Bug Mosi_6215-Mar-16 9:52 Mosi_62 15-Mar-16 9:52
 Re: how to download the source code Mosi_6214-Mar-16 9:16 Mosi_62 14-Mar-16 9:16
 Bug MartinBell15-Dec-15 6:34 MartinBell 15-Dec-15 6:34
 Re: Bug Mosi_6214-Mar-16 9:37 Mosi_62 14-Mar-16 9:37
 Re: Bug Member 1285975619-Nov-16 10:53 Member 12859756 19-Nov-16 10:53
 Discussion of the cubic roots YvesDaoust21-Jul-14 3:04 YvesDaoust 21-Jul-14 3:04
 Re: Discussion of the cubic roots Mosi_6222-Jul-14 8:35 Mosi_62 22-Jul-14 8:35
 I'll take your word for it! .dan.g.20-Jul-14 20:15 .dan.g. 20-Jul-14 20:15
 Vote of 5 Kenneth Haugland18-Jul-14 18:04 Kenneth Haugland 18-Jul-14 18:04
 Re: Vote of 5 Mosi_6219-Jul-14 5:40 Mosi_62 19-Jul-14 5:40
 Last Visit: 31-Dec-99 18:00     Last Update: 14-Jun-21 23:46 Refresh 1