Tilted Forum Project Discussion Community

Tilted Forum Project Discussion Community (https://thetfp.com/tfp/)
-   Tilted Knowledge and How-To (https://thetfp.com/tfp/tilted-knowledge-how/)
-   -   Algorithm for exact roots of Cubic Function (https://thetfp.com/tfp/tilted-knowledge-how/129371-algorithm-exact-roots-cubic-function.html)

Hain 12-25-2007 09:26 AM

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)*(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:
Code:

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?

Hain 12-26-2007 09:51 AM

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.

albania 12-27-2007 05:08 PM

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 

Hain 01-03-2008 01:11 PM

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

albania 01-10-2008 08:36 PM

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.

Hain 01-16-2008 02:08 AM

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.


All times are GMT -8. The time now is 02:57 PM.

Powered by vBulletin® Version 3.8.7
Copyright ©2000 - 2024, 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 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360