I'm dealing with a modelization problem and i've arrived at one point where i need to calculate the roots of a cubic polynomial without specifing their value: $$ Px^{3}-(Pb+RT)x^{2}+ax-ab=0 $$ I know the value of:
$$a=0.244 \\b=0.026 \\R=0.082$$ But the only information i have about T (temperature, K) and P (pressure, atm) is that: $$T>0,P>0$$
With the help of python library sympy i was able to find the roots analytically (result here), this one is the real root (and the only one i'm interested for my purposes), however my professor suggested to use a numerical method rather than an analytical one.
- My question is: is it possible to calculate the roots of this polynomial numerically without specifying their value? And if so, how?
With $x=V$, you are solving Van der Waals cubic equation of state $$P=\frac{R T}{V-b}-\frac{a}{V^2}$$ Considering the values given for $a$ and $b$, I am ready to bet that this is for hydrogen.
If you are sure that there is only one real root, you can have much simpler formulae.
Let $$B=-\frac{R T}{P}-b \qquad \qquad C=\frac a P \qquad \qquad D=- \frac {a b} P$$ Define $$p=C-\frac{1}{3}B^2 \qquad\qquad \text{and} \qquad\qquad q=\frac{2 }{27}B^3-\frac{1}{3}BC+D$$
If $\color{red} {p>0}$, then $$\color{blue} {x=-2\sqrt{\frac{p}{3}}\sinh\left[\frac{1}{3}\operatorname{arsinh}\left(\frac{3q}{2p}\sqrt{\frac{3}{p}}\right)\right]-\frac 13 B}$$
If $\color{red} {p<0}$, then $$\color{blue} {x=-2\frac{|q|}{q}\sqrt{-\frac{p}{3}}\cosh\left[\frac{1}{3}\operatorname{arcosh}\left(-\frac{3|q|}{2p}\sqrt{\frac{-3}{p}}\right)\right]-\frac 13 B}$$
This formulae are simple and, using them, you do not need any iterative method required by a numerical method.
If you want a numerical approach , Newton method would be the simplest but you need a starting guess $x_0$. Taking into account that for the gas phase, the attraction term $\frac{a}{V^2}$ is very small, ignore it and use $$x_0=b+\frac{R T}{P}$$
Let me try an example with $T=345$ and $P=123$. The iterates will be $$\left( \begin{array}{cc} n & x_n \\ 0 & 0.256000000000000 \\ 1 & 0.249242566341334 \\ 2 & 0.248862718778783 \\ 3 & 0.248861544341478 \\ 4 & 0.248861544330268 \end{array} \right)$$ while working with rational values for $(a,b,R)$ the result is $$x=\frac 1{375}\Bigg[32+\sqrt{\frac{152686}{41}} \cosh \left(\frac{1}{3} \cosh ^{-1}\left(\frac{19459933}{305372}\sqrt{\frac{41}{152686}}\right)\right) \Bigg]$$ which is the exact solution.
If you do not need an extravagant accuracy, using the same $x_0$, the first iterate of Halley method will be $$x_1=\frac{b P+R T}{P}-\frac{a R T \left(a P+(b P+R T)^2\right)}{a^2 P^2+2 a b P^2 (b P+R T)+(b P+R T)^4}$$ Applied to the example, this gives $$x_1=\frac{7952127503773}{31951960920500}=0.24887760$$