Solution Techniques for Cubic Expressions and Root Finding

Hopefully, we have convinced you that the use of cubic equations of state can represent a very meaningful and advantageous way of modeling the PVT behavior of petroleum fluids. What we need now are the tools that will allow us to get the information that we want out of them. Even though cubic equations of state are explicit in pressure, pressure is not the common unknown to be calculated in the typical problem. In the most common problem, pressure and temperature are known and we want either molar volume (or its reciprocal, molar density) or compressibility factor (the most likely case). Therefore, we are faced very often with the need to solve for the roots of a cubic expression. Here we present a number of approaches that may be followed.

Analytical scheme

Given the cubic polynomial with real coefficients: x³ + ax² + bx + c = 0, the first step is to calculate the parameters:

$Q = \frac{a^{2} - 3 b}{9} and R= \frac{2 a^{3} - 9 a b + 27 c}{54}$ (11.5)

Now let M = R² – Q³ be the discriminant. We then consider the following cases:

If M < 0 (R² < Q³), the polynomial has three real roots. For this case, compute $θ = \arccos (\frac{R}{\sqrt{Q^{3}}})$ and calculate the three distinct real roots as:

$x_{1} = - (2 \sqrt{Q} \cos \frac{θ}{3}) - \frac{a}{3}$ (11.6a)

$x_{2} = - (2 \sqrt{Q} \cos \frac{θ + 2 π}{3}) - \frac{a}{3}$ (11.6b)

$x_{3} = - (2 \sqrt{Q} \cos \frac{θ - 2 π}{3}) - \frac{a}{3}$ (11.6c)

Note that x₁, x₂, x₃ are not given in any special order, and that $θ$ has to be calculated in radians.
If M > 0 (R² > Q³), the polynomial has only one real root. Compute:

$S = \sqrt[3]{- R + \sqrt{M}}$ (11.7a)

$T = \sqrt[3]{- R - \sqrt{M}}$ (11.7b)
and calculate the real root as follows:
$x_{1} = S + T - \frac{a}{3}$ (11.7c)

Two complex roots (complex conjugates) may be found as well. However, they are of no interest for our purposes, and thus no formulas are provided. Such formulas can be found in the following suggested readings:

W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, 2^nd Edition, Cambridge Univ. Press, (1992), p. 179.
Spiegel, M., Liu, J., Mathematical handbook of formulas and tables, 2^nd Edition, Schaum’s Outline Series, McGraw Hill, p.10.

Sometimes, the equations for S and T listed above cause problems while programming. This usually happens whenever the computer/calculator performs the cubic root of a negative quantity. If you want to avoid such a situation, you may compute S’ and T’ instead:

$S' = - s i g n (R) \sqrt[3]{a b s (R) + \sqrt{M}}$ $T' = Q / S'$ (making T’=0 when S’=0)

where:

abs(R) = Absolute value of R

and

sign(R) = (+1) or (– 1) if R is positive or negative respectively.

It may be defined as:

sign(R) = R/Abs(R)

and then the real root is:

$x_{1} = S' + T' - \frac{a}{3}$

Keep in mind the following useful relationships among the roots of any cubic expression:

$x_{1} + x_{2} + x_{3} = - a$ (11.8a)

$x_{1} x_{2} + x_{2} x_{3} + x_{3} x_{1} = + b$ (11.8b)

$x_{1} x_{2} x_{3} = - c$ (11.8c)

Self-Check

Test your understanding of cubic root calculations by analytical means by solving the following examples.

With general cubic expressions,

\begin{array}{l} x^{3} - 6 x^{2} + 11 x - 6 = 0 \Rightarrow \\ x_{1} = 1 \\ x_{2} = 2 \\ x_{3} = 3 \end{array}

\begin{array}{l} x^{3} + 7 x^{2} + 49 x + 343 = 0 \\ x_{1} = - 7 \\ x_{2} = 0 + 7 i \\ x_{3} = 0 - 7 i \end{array}

\begin{array}{l} x^{3} + 2 x^{2} + 3 x + 4 = 0 \Rightarrow \\ x_{1} = - 1.65063 \\ x_{2} = - 0.174685 + 1.54687 i \\ x_{3} = - 0.174685 - 1.54687 i \end{array}

with EOS in terms of volume (v), for a pure component,

$v^{3} - 7.8693 v^{2} + 13.3771 v - 6.5354 = 0$

=> One real root (one phase),

$v_{1} = 5.7357$

Another example,

$v^{3} - 15.6368 v^{2} + 30.315 v - 14.8104 = 0$

=> Two possible phases (three real roots),

$\begin{array}{l} v_{1} = 0.807582 (liquid phase) \\ v_{2} = 1.36174 (rejected) \\ v_{3} = 13.4675 (gas phase) \end{array}$

with EOS in terms of the compressibility factor (z), for a pure component,

$z^{3} - 1.0595 z^{2} + 0.2215 z - 0.01317 = 0$

=> One phase,

$z = 0.8045$

Another example,

$z^{3} - z^{2} + 0.089 z - 0.0013 = 0$

=> Two possible phases,

$\begin{array}{l} z_{l i q u i d} = 0.0183012 \\ z_{x} = 0.0786609 (useless) \\ z_{vapor} = 0.903038 \end{array}$

Numerical Scheme

The Newton-Raphson method provides a useful scheme for solving for a non-explicit variable from any form of equation (not only cubic ones). Newton-Raphson is an iterative procedure with a fast convergence, although it is not always capable of providing an answer — because a first guess close enough to the actual answer must be provided.

In solving for “x” in any equation of the type f(x)=0, the method provides a new estimate (new guess) closer to the actual answer, based on the previous estimate (old guess). This is written as follows:

$x_{n e w} = x_{o l d} - \frac{f (x_{o l d})}{f' (x_{o l d})}$ (11.9)

Considering a cubic equation $f (x) = x^{3} + a x^{2} + b x + c = 0$ , the previous equation takes the form:

$x_{n e w} = x_{o l d} - \frac{x_{o l d}^{3} + a x_{o l d}^{2} + b x_{o l d} + c}{3 x_{o l d}^{2} + 2 a x_{o l d} + b}$ (11.10)

The iterations continue until no significant improvement for “x_new” is achieved, i.e., | x_new – x_old | < tolerance. An educated guess must be provided as the starting value for the iterations. If you are solving a cubic equation in Z (compressibility factor), it is usually recommend to take Z = bP/RT as the starting guess for the compressibility of the liquid phase and Z = 1 for the vapor root.

Semi-analytical Scheme

If you used the previous numerical approach to calculate one of the roots of the cubic expression, the semi-analytical scheme can give you the other two real roots (when they exist). By using the relationships given before, with the value ‘x₁’ as the root already known, the other two roots are calculated by solving the system of equations:

$x_{2} + x_{3} = - a - x_{1}$ (11.11a)

$x_{2} x_{3} = - c / x_{1}$ (11.11b)

which leads to a quadratic expression.

This procedure can be reduced to the following steps:

Let $x^{3} + a x^{2} + b x + c = 0$ be the original cubic polynomial and “E” the root which is already known (x₁ = E). Then, we may factorize such a cubic expression as:

$(x - E) (x^{2} + F x + G) = 0$ (11.12a)
where:

$F = a + E$ (11.12b)

$G = - c / E$ (11.12c)
Solve for x₂, x₃ by using the quadratic expression formulae,

$x_{1} = E$ (11.13a)

$x_{2} = \frac{- F + \sqrt{F^{2} - 4 G}}{2}$ (11.13b)

$x_{3} = \frac{- F - \sqrt{F^{2} - 4 G}}{2}$ (11.13c)

Analytical scheme

Self-Check

Numerical Scheme

Navigation

EMS

Programs

Related Links