This is based on the discussion in #254

This issue seems to affect most phase equilibrium algorithms in feos. We focus on the most important ones

Critical points ✅

The issue is related to the eigenvectors of the 0 components, as well as their composition going to 0 which creates numerical problems. I did not find a way to circumvent this by rewriting the expressions, but we can simply catch the special case (although that might effect the calculation of derivatives)

• Fixed the evaluation of eigenvectors of 2x2 matrices in num-dual Fix the calculation of eigenvectors of 2x2 matrices with zeros as elements itt-ustutt/num-dual#98
• Manually overrode the components which have 0 composition in the calculation of the 2nd criticality condition.

Bubble and dew points ℹ️

There are actually a number of points where x_i=0 leads to failure:

• The definition of the residual in the outer loop |(k_i x_i)/y_i| -> could be changed to absolute error
• The residual in the newton step mu_i^L=mu_i^V due to to the divergence in the ideal gas term -> could be changed to isofugacity
• related: Other divisions by 0/divergences in the ideal gas term in the calculation of the Jacobian -> should be gone when switching to fugacities

Tp flash ✅

An actual bug that set `ln(x/y)` to zero instead of the residual `ln(x phi_l/(y phi_v))`. Fixed by changing the residual (which is only used for determining convergence) from `ln(x phi_l/(y phi_v))` to `x phi_l/phi_v - y`.