S.O.S. Mathematics CyberBoard

Hi all,

I need to write an algorithm which finds the real roots to a cubic equation.

I am applying a cubic equation of state to a chemical engineering problem, and a specific cubic equation needs to be solved thousands of times over the course of one routine.

Therefore it needs to be COMPLETELY robust, and here's where my problem lies.

It seems analytical methods are prone to errors due to the floating point precision of computers, and I am tending towards several numerical approximation methods (Newton-Raphson, Halley's method, Brent's Method). The problems are:

Newton-Raphson / Halley's - Not sure whether these will be guaranteed to converge, they may diverge if the initial guess is not good enough (which it may well not be)

Brent's Method - although convergence is guaranteed, I am unsure of how to bracket the roots and whether it will be affected by the presence of complex roots.

ANy help or pointers much appreciated. As I say I am not simply trying to solve a cubic equation (fairly easy) but write a code which will NEVER fail no matter what values are used. Thanks.

If you want it to be completely robust, you need to split it into cases. For real coefficients, at the very least you will have three cases to consider, according to the sign of the discriminant, and make sure you avoid subtracting nearly equal quantities in the intermediate calculations.

_________________