12252007, 08:26 AM  #1 (permalink) 
has a plan
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: Code:
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)*(st)*j; x3 = (1/2)*(s+t)  (b / (3*a))  (sqrt(3)/2)*(st)*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: Code:
x1 = 0.5000  0.2887i x2 = 0.5000  0.2887i x3 = 1.0000 + 0.5774i Can someone throw a new exact algorithm at me? One that doesn't use loops to approximate?
__________________
Last edited by Hain; 12262007 at 04:06 AM.. 
12262007, 08:51 AM  #2 (permalink) 
has a plan
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; 12262007 at 12:19 PM.. 
12272007, 04:08 PM  #3 (permalink) 
Psycho

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 highpowered 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 
01032008, 12:11 PM  #4 (permalink) 
has a plan
Location: middle of Whywouldanyonebethere

Solutions
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 ohso 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 MFile 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; 01032008 at 03:24 PM.. 
01102008, 07:36 PM  #5 (permalink) 
Psycho

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.

01162008, 01:08 AM  #6 (permalink) 
has a plan
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.
__________________

Tags 
algorithm, cubic, exact, function, roots 
Thread Tools  

