Two residuals are
Equating pressure of the two phases
$r_1 = \delta''[1+\delta''\alpha_{\delta}^r(\tau,\delta'')]-\delta'[1+\delta'\alpha_{\delta}^r(\tau,\delta')]$
Equating Gibbs energy of the two phases
$r_2 = \delta''\alpha_{\delta}^r(\tau,\delta'')+\alpha^r(\tau,\delta'')+\ln\delta''-[\delta'\alpha_{\delta}^r(\tau,\delta')+\alpha^r(\tau,\delta')+\ln\delta']$
Independent (unknown) variables are either $\tau,\delta'$ or $\tau,\delta''$. Obtain all the partial derivatives:
$\left.\displaystyle\frac{\partial r_1}{\partial \tau}\right|_{\delta',\delta''} = (\delta'')^2\alpha^r_{\delta\tau}(\tau,\delta'')-(\delta')^2\alpha^r_{\delta\tau}(\tau,\delta')$
$\left.\displaystyle\frac{\partial r_1}{\partial \delta''}\right|_{\tau,\delta'} = [1+\delta''\alpha^r_{\delta}(\tau,\delta'')]+\delta''[\delta''\alpha^r_{\delta\delta}(\tau,\delta'')+\alpha^r_{\delta}(\tau,\delta'')] = 1+2\delta''\alpha^r_{\delta}(\tau,\delta'')+(\delta'')^2\alpha^r_{\delta\delta}(\tau,\delta'')$
$\left.\displaystyle\frac{\partial r_1}{\partial \delta'}\right|_{\tau,\delta''} = -[1+\delta'\alpha^r_{\delta}(\tau,\delta')]-\delta'[\delta'\alpha^r_{\delta\delta}(\tau,\delta')+\alpha^r_{\delta}(\tau,\delta')] = -1-2\delta'\alpha^r_{\delta}(\tau,\delta')-(\delta')^2\alpha^r_{\delta\delta}(\tau,\delta')$
$\left.\displaystyle\frac{\partial r_2}{\partial \tau}\right|_{\delta',\delta''} = \delta''\alpha^r_{\delta\tau}(\tau,\delta'') + \alpha^r_{\tau}(\tau,\delta'')-[\delta'\alpha^r_{\delta\tau}(\tau,\delta') + \alpha^r_{\tau}(\tau,\delta')]$
$\left.\displaystyle\frac{\partial r_2}{\partial \delta''}\right|_{\tau,\delta'} = \delta''\alpha^r_{\delta\delta}(\tau,\delta'')+\alpha^r_{\delta}(\tau,\delta'') +\alpha^r_{\delta}(\tau,\delta'') +\displaystyle\frac{1}{\delta''} = \delta''\alpha^r_{\delta\delta}(\tau,\delta'')+2\alpha^r_{\delta}(\tau,\delta'') +\displaystyle\frac{1}{\delta''} $
$\left.\displaystyle\frac{\partial r_2}{\partial \delta'}\right|_{\tau,\delta''} = -\delta'\alpha^r_{\delta\delta}(\tau,\delta')-\alpha^r_{\delta}(\tau,\delta') -\alpha^r_{\delta}(\tau,\delta') -\displaystyle\frac{1}{\delta'} = -\delta'\alpha^r_{\delta\delta}(\tau,\delta')-2\alpha^r_{\delta}(\tau,\delta') -\displaystyle\frac{1}{\delta'} $ If $\delta'$ is imposed, the Jacobian is given by
$\mathbf{J} = \left[ \begin{array}{cc} \left.\displaystyle\frac{\partial r_1}{\partial \tau}\right|_{\delta',\delta''} & \left.\displaystyle\frac{\partial r_1}{\partial \delta''}\right|_{\tau,\delta'} \\ \left.\displaystyle\frac{\partial r_2}{\partial \tau}\right|_{\delta',\delta''} & \left.\displaystyle\frac{\partial r_2}{\partial \delta''}\right|_{\tau,\delta'}\end{array} \right]$
with the increment $\mathbf{v}$ obtained from
$\mathbf{J}\left[\begin{array}{c}\Delta \tau \\ \Delta\delta''\end{array}\right] = -\left[\begin{array}{c}r_1 \\ r_2\end{array}\right]$
If $\delta''$ is imposed, the Jacobian is given by
$\mathbf{J} = \left[ \begin{array}{cc} \left.\displaystyle\frac{\partial r_1}{\partial \tau}\right|_{\delta',\delta''} & \left.\displaystyle\frac{\partial r_1}{\partial \delta'}\right|_{\tau,\delta''} \\ \left.\displaystyle\frac{\partial r_2}{\partial \tau}\right|_{\delta',\delta''} & \left.\displaystyle\frac{\partial r_2}{\partial \delta'}\right|_{\tau,\delta''}\end{array} \right]$
with the increment $\mathbf{v}$ obtained from
$\mathbf{J}\left[\begin{array}{c}\Delta \tau \\ \Delta\delta'\end{array}\right] = -\left[\begin{array}{c}r_1 \\ r_2\end{array}\right]$ The method needs to start off with decent guesses for densities and temperature, which can be obtained by the ancillaries. First solve $\rho(T)=\rho_{spec}$ to get $T$, then use $T$ to get the other density using the ancillary. And then kick into the full Newton-Raphson solution. Or alternatively, use $\ln\delta'$ and $\ln\delta''$ as independent variables, would bring liquid an vapor density values closer in magnitude. Rather than being 16 orders of magnitude different, they would be much closer to one order of magnitude in difference. For instance:
In [1]:
# Propane at triple point
rhoL = 16625.80554 #mol/m3
rhoV = 0.0000002419454458 #mol/m3
from math import log10
print log10(rhoL), log10(rhoV)
That should help with the truncation error. But make sure to distribute through the \delta in order to avoid truncation in internal calculation
Three residuals are
Equating pressure of the two phases
$r_1 = \delta''[1+\delta''\alpha_{\delta}^r(\tau,\delta'')]-\delta'[1+\delta'\alpha_{\delta}^r(\tau,\delta')]$
Equating Gibbs energy of the two phases
$r_2 = \delta''\alpha_{\delta}^r(\tau,\delta'')+\alpha^r(\tau,\delta'')+\ln\delta''-[\delta'\alpha_{\delta}^r(\tau,\delta')+\alpha^r(\tau,\delta')+\ln\delta']$
Specified value must equal the desired value
$r_3 = \ln(X_S) - \ln(X_{specified})$ for $p$
$r_3 = X_S - X_{specified}$ for $h$, $s$, $u$
Three unknowns are $\tau$, $\ln\delta'$, $\ln\delta''$.
For the first two residuals, the partials can be obtained in the same way as for density imposed
$\left.\displaystyle\frac{\partial r_1}{\partial \tau}\right|_{\delta',\delta''} = (\delta'')^2\alpha^r_{\delta\tau}(\tau,\delta'')-(\delta')^2\alpha^r_{\delta\tau}(\tau,\delta')$
$\left.\displaystyle\frac{\partial r_1}{\partial \delta''}\right|_{\tau,\delta'} = 1+2\delta''\alpha^r_{\delta}(\tau,\delta'')+(\delta'')^2\alpha^r_{\delta\delta}(\tau,\delta'')$
$\left.\displaystyle\frac{\partial r_1}{\partial \delta'}\right|_{\tau,\delta''} = -1-2\delta'\alpha^r_{\delta}(\tau,\delta')-(\delta')^2\alpha^r_{\delta\delta}(\tau,\delta')$
$\left.\displaystyle\frac{\partial r_2}{\partial \tau}\right|_{\delta',\delta''} = \delta''\alpha^r_{\delta\tau}(\tau,\delta'') + \alpha^r_{\tau}(\tau,\delta'')-[\delta'\alpha^r_{\delta\tau}(\tau,\delta') + \alpha^r_{\tau}(\tau,\delta')]$
$\left.\displaystyle\frac{\partial r_2}{\partial \delta''}\right|_{\tau,\delta'} = \delta''\alpha^r_{\delta\delta}(\tau,\delta'')+2\alpha^r_{\delta}(\tau,\delta'') +\displaystyle\frac{1}{\delta''} $
$\left.\displaystyle\frac{\partial r_2}{\partial \delta'}\right|_{\tau,\delta''} = -\delta'\alpha^r_{\delta\delta}(\tau,\delta')-2\alpha^r_{\delta}(\tau,\delta') -\displaystyle\frac{1}{\delta'} $
The third residual's partials when one of h,s,u are specified can be obtained from
$\left.\displaystyle\frac{\partial r_3}{\partial \tau}\right|_{\delta', \delta''} = \left.\displaystyle\frac{\partial X_S}{\partial \tau}\right|_{\delta',\delta''}$
$\left.\displaystyle\frac{\partial r_3}{\partial \delta'}\right|_{\tau,\delta''} = \left.\displaystyle\frac{\partial X_S}{\partial \delta'}\right|_{\tau,\delta''}$
$\left.\displaystyle\frac{\partial r_3}{\partial \delta''}\right|_{\tau,\delta'} = \left.\displaystyle\frac{\partial X_S}{\partial \delta''}\right|_{\tau,\delta'}$
Or when $\ln p$ is specified
$\left.\displaystyle\frac{\partial r_3}{\partial \tau}\right|_{\delta', \delta''} = \dfrac{1}{X_S}\left.\dfrac{\partial X_S}{\partial \tau}\right|_{\delta',\delta''}$
$\left.\displaystyle\frac{\partial r_3}{\partial \delta'}\right|_{\tau,\delta''} = \dfrac{1}{X_S}\left.\dfrac{\partial X_S}{\partial \delta'}\right|_{\tau,\delta''}$
$\left.\displaystyle\frac{\partial r_3}{\partial \delta''}\right|_{\tau,\delta'} = \dfrac{1}{X_S}\left.\dfrac{\partial X_S}{\partial \delta''}\right|_{\tau,\delta'}$