3 replies to this topic

[adw07]

Hi all,

I am having a lot of trouble trying to solve the cubic equation of state (for compressibility factor):

f(Z) = Z^3 + (B-1)Z^2 + (A - 2B - 3B^2)Z + (B^3 + B^2 - AB) = 0

Cardano may have worked out the first way to do this in 1545 but I cannot seem to find a COMPLETELY robust way of doing it even with the wealth of knowledge I can find on the internet, so please any suggestions appreciated.

Analytical methods such as trigonometric, etc. seem to suffer from the floating point precision accuracy of computers and so I have ruled these out. (They work for a lot of cases, but are not 100% foolproof)

Newton-Raphson & similiar methods work well for a lot of situations, however they need a good estimate (not necessarily known) and can still fail to converge in many situations and there are problems when the derivative is near zero, etc.

My current line of thinking is to use Brent's method (combining bisection, secant & inverse quadratic interpolation methods) to guarantee a solution. However, I am unsure about which initial interval limits to set.

The program I am ultimately trying to create will have to solve the cubic equation of state thousands of times over the course of one routine, so I need a completely robust way of solving it. Any help or pointers very much appreciated, thanks.

[kkala]

I think cubic equation can be solved analytically, that is straightforward using Algebra. A Fortran program was elaborated by me in 1973 to solve it, and it was most probably right, since no comments were received. Though I cannot recall details, it was application of the published formulas but under a lot of conditional cases. http://en.wikipedia..../Cubic_function may help to elaborate this solution. Legth of formulas can be disappointing at first look, but big parts of them are repeatedly used making solving easier.

Edited by kkala, 01 September 2011 - 06:45 AM.

[suhas]

Dear Adw
If you want to solve numerically use the Newton raphson method to solve. Just ensure two things, use analytical derivative, i.e differentiate above eqn wrt z and use this eqn to calculate f'(z) in newton raphson method.
For initial estimate use the Z= 1 as this means gas is ideal. Actual value of z differs depends on the non-ideality.

I hope this will solve your problem

[adw07]

Thanks guys, I have developed a custom method which uses all three householder's methods (Newton-Raphson, Halley's & Third Order) simultaneously, with multiple initial guesses(0.5, 1,0,5,10). In practice this technique usually converges in less than 10 iterations. I have put in a bisection backup routine just in case all householder's methods fail (extremely unlikely), and so overall should be a pretty robust routine.