Tilted Forum Project Discussion Community  

Go Back   Tilted Forum Project Discussion Community > The Academy > Tilted Knowledge and How-To

LinkBack Thread Tools
Old 12-25-2007, 09:26 AM   #1 (permalink)
has a plan
Hain's Avatar
Location: middle of Whywouldanyonebethere
Algorithm for exact roots of Cubic Function

Short of the long of it:
I need a function that gives the exact roots of a cubic function. I was using this to get exact values, however some tests have shown obvious failures in the algorithm. Given:
f(x) = a x<sup>3</sup> + b x<sup>2</sup> + c x + d = a (x - x<sub>1</sub>) (x - x<sub>2</sub>) (x - x<sub>3</sub>) ~ the algorithm is:
q = ( (3*a*c) - (b^2) ) / (9 * (a^2) );
r = ( (9*a*b*c) - (27*(a^2)*d ) - 2*(b^3)) / (54*(a^3));

s = (r + sqrt(q^3 + r^2))^(1/3);
t = (r - sqrt(q^3 + r^2))^(1/3);

x1 = s + t - (b / (3*a));
x2 = -(1/2)*(s+t) - (b / (3*a)) + (sqrt(3)/2)*(s-t)*j;
x3 = -(1/2)*(s+t) - (b / (3*a)) - (sqrt(3)/2)*(s-t)*j;

Long of it:
I am modeling data that requires all roots of a cubic function, even imaginary if they occur. I liked the one outlined in the link provided however as it doesn't work with all numbers I throw at it.

I throw a = 1, b = -2, c = 1, d = 0, I get answers like this:
x1 = 0.5000 - 0.2887i
x2 = 0.5000 - 0.2887i
x3 = 1.0000 + 0.5774i
The answers should be 0, 1, 1. If you remember your algebra II class, you know that every cubic function with real coefficients will always have one real root. It has worked for almost all other functions I throw at it... Here it doesn't work.

Can someone throw a new exact algorithm at me? One that doesn't use loops to approximate?

Last edited by Hain; 12-26-2007 at 05:06 AM..
Hain is offline  
Old 12-26-2007, 09:51 AM   #2 (permalink)
has a plan
Hain's Avatar
Location: middle of Whywouldanyonebethere
I know what the problem is.

Both the TI 89 and MATLAB use a shortcut when it comes to roots of numbers.

If you have the cube root of n, that can be written as: n^(1/3)
which can be computed like this: e^((1/3)*ln(n)) ...

That is why this program fails. Any different equations?


With a few handwaves I made the program works.

Last edited by Hain; 12-26-2007 at 01:19 PM..
Hain is offline  
Old 12-27-2007, 05:08 PM   #3 (permalink)
albania's Avatar
So are you saying that itís running into trouble when n is negative, as it would be in this case? Interesting, I wouldnít think something as high-powered as matlab would have that trouble. If that were truly the case then simply flipping the sign on n if itís negative taking the cube root then flipping the sign back would work. I would think your problem is elsewhere perhaps in the code or a negative sign. It may help you if you display the specific values for q, r, s, and t, and see exactly where it goes wrong.

Anywhoo, I was bored so I checked your formula for a couple of cases including the one you said your program didnít handle. For those cases it worked(engineering induction therefore guarantees us that it will always work). It seems as though that the formula you found is the only closed formed solution that is out there, so I donít think youíre going to have any luck on that front. There may be others but I think theyíre probably technical and not closed form, i.e. over my head. Well seeing as I was bored I plugged this into Mathematica which has a root function; this is what I got (my best guess is that it is what you have up there, though the simplify command seemed to do little):
click for computational algebra orgasm   click to show 
albania is offline  
Old 01-03-2008, 01:11 PM   #4 (permalink)
has a plan
Hain's Avatar
Location: middle of Whywouldanyonebethere

I had my friend test my original algorithm out with Mathematica and it worked every time, even when some of the values were negative. I tried the roundabout log function for cube roots and it output the same complex root as taking x^(1/3). I guess the code writers of MATLAB thought they were being oh-so clever using one algorithm that just calls another... I would too in most cases.

I rewrote the code, and added some tolerances (+/- 10^-6 of the imaginary part) so that this does not occur and can be used in most other applications.

If you care, here are my MATLAB codes:
cube root calculator   click to show 

cubic solver, cardano   click to show 

If you need it, for nothing else than educational purposes:
I was backwardly calculating transverse deflections of loaded beam. All of my data was in excel spreadsheets, but Excel is not capable of performing these kinds of calculations. This M-File will take data from XLS and XLSX spreadsheets and compute data, and put it back into the spreadsheets.
strain_analysis   click to show 

Last edited by Hain; 01-03-2008 at 04:24 PM..
Hain is offline  
Old 01-10-2008, 08:36 PM   #5 (permalink)
albania's Avatar
You write very clean code; that was easy to follow. My attempts at the computational side of things seem more like the computer forays of a very angry chimpanzee by comparison. I suppose that has to do with my lack of formal training and the development of skill in the area as necessity instead.
albania is offline  
Old 01-16-2008, 02:08 AM   #6 (permalink)
has a plan
Hain's Avatar
Location: middle of Whywouldanyonebethere
Looks like there might be something wrong with the algorithm... I have no idea where, but I bet it has to do with that cube root algorithm again!

I will try it out again and remove that code that determines which cube root's angle is closest to the original complex input.
Hain is offline  

algorithm, cubic, exact, function, roots

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

All times are GMT -8. The time now is 09:28 AM.

Tilted Forum Project

Powered by vBulletin® Version 3.8.7
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Search Engine Optimization by vBSEO 3.6.0 PL2
© 2002-2012 Tilted Forum Project

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73