I have to solve an equation $$\sum_{i=1}^N x_i = \sum_{i=1}^N y_i,$$ where
$$x_i = \frac{z_i}{1 + (K_i - 1) w}$$ and $$y_i = \frac{K_i z_i}{1 + (K_i - 1) w}.$$

The $z_i$ are all positive and add to $1$; the $K_i$ are positive but range from very small to very large (at least, one is strictly greater than $1$ and one is strictly smaller than $1$).

The solution ($w$) looked for is between two vertical asymptotes corresponding to $-1 / (K_{\text{max}} - 1)$ (which is negative) and $1 / (1 - K_{\text{min}})$ (which is greater than $1$). Considering the fact that, at solution, all $x_i$ and $y_i$ must be positive and smaller than $1$ (since they both must add to $1$), I have been able to find a left bound corresponding to the maximum value of $\dfrac{K_i z_i - 1}{K_i-1}$ (considering only the $K_i > 1$) and a right bound corresponding to the minimum value of $\dfrac{1 - z_i}{1 - K_i}$ (considering only the $K_i < 1$).

I have not been able to go further and would appreciate some help.

  • $\begingroup$ I think it is useful to post the following reformulation: we are looking for an approximation of $w$ such that $$\left.\frac{d}{dx}\prod_i\left(1+(K_i-1)\,w\, e^{z_i x}\right)\right|_{x=0}=0,$$ i.e. an approximation of $w$ for which $x=0$ is a stationary point of $f_w(x)=\prod_i\left(1+w\,(K_i-1)\,e^{z_i x}\right).$ $\endgroup$ – Jack D'Aurizio Jan 10 '14 at 17:09
  • $\begingroup$ It is essential that some $K_i$ is below $1$, otherwise, by translation super-additivity of the geometric mean, $$f_w(x)\geq 1 + w\,e^x\prod_i(K_i-1).$$ $\endgroup$ – Jack D'Aurizio Jan 10 '14 at 17:15
  • $\begingroup$ Another equivalent form is: $$\left.\frac{d}{dx}\prod_i\left(1+w(K_i-1)(1+x z_i)\right)\right|_{x=0}=0,$$ or $$\left.\frac{d}{dx}\sum_i\log\left(1+w(K_i-1)(1+x z_i)\right)\right|_{x=0}=0.$$ Now we can expand the logarithm as a Taylor series, or use a suitable polynomial approximation of $\log(1+x)$ in order to approximate $w$. $\endgroup$ – Jack D'Aurizio Jan 10 '14 at 17:25
  • $\begingroup$ Could you post the details of the original problem for the sake of completeness? I strongly believe that this question comes from an eigenvalues approximation problem with preconditioning for a certain (maybe stochastic) matrix. $\endgroup$ – Jack D'Aurizio Jan 10 '14 at 17:38
  • $\begingroup$ Also consider that the algebraic equation $$\sum_i\frac{(K_i-1)z_i}{1+(K_i-1)w}=0$$ can be solved through Newton's method, so (computationally speaking) we only need a reasonable starting point. $\endgroup$ – Jack D'Aurizio Jan 10 '14 at 18:48

We are looking for a root of $$f(w) = \sum_{i=1}^{n}\frac{(K_i-1)z_i}{1+(K_i-1)w}.$$ For the sake of brevity, we put $A_i=\frac{1}{K_i-1}$ and assume $A_1<\ldots<A_n$. $$f(w) = \sum_{i=1}^{n}\frac{z_i}{A_i+w}$$ changes its sign in every $w_i=-A_i$ and $\lim_{w\to\infty}f(w)=0$, so it has a root (and only one) in every interval $(-A_{j+1},-A_j)$. Moreover, the roots of $f(w)$ are the same of the roots of the $(n-1)$-th degree polynomial $$p(w) = \sum_{i=1}^{n}z_i\prod_{k\neq i}(w+A_k), $$ so we can find our $w$ by looking for the only root of $p(w)$ in $(-A_{j+1},-A_j)$. Many computational methods can accomplish this task: bisection, Newton's, secants or a suitable combination of these ones.

Update. Since the $A_i$s can be only positive or less than $-1$, let $I^{+}$ be the set of $i$ such that $A_i>0$ and $I^{-}$ be the set of $i$ such that $A_i<0$. Moreover, let $a=-\frac{1}{K_{max}-1}<0$ and $b=\frac{1}{1-K_{min}}>1$. We are looking for the only solution in $(a,b)$ of: $$L(x)=-\sum_{i\in I^{-}}\frac{z_i}{x+A_i}=\sum_{i\in I^+}\frac{z_i}{x+A_i}=R(x),$$ where $L(x)$ is a convex increasing function and $R(x)$ is a convex decreasing function on $(a,b)$. We can perform a few bisection steps in order to ensure that $R(x)$ and $L(x)$ are intersecting convex, monotonic and bounded functions over $[c,d]\subset(a,b)$. The problem (*) is now to find the point of intersection between the graphics of two bounded convex functions, with opposite monotonicity, over $[c,d]$. For this kind of problem, we can tweak the secant-tangent method in order to ensure fast convergence, due to the following simple geometric fact: Intersecting the graphics of two convex bounded function with opposite monotonicity

The point of intersection $P=(p_x,p_y)$ belong to the interior part of the triangle having two tangents and one secant of $R(x)$ as sides. The same holds for $L(x)$, so by intersecting the $L(x)$-secant with the $R(x)$-tangent in $d$, and the $R(x)$-secant with the $L(x)$-tangent in $d$, we can find a sub-interval $[c',d']\subset[c,d]$ for which $p_x\in[c',d']$. Moreover, one between $c'$ and $d'$ can be replaced with the abscissa of the intersection between the $L(x)$-secant and the $R(x)$-secant, as clearly depicted above. (Other choices are the abscissa of the intersection of the tangents in $c$, or in $d\ldots$)

To avoid the first bisection steps, we can multiply both $R(x)$ and $L(x)$ by $(x-a)(x-b)$ in order to remove the singularities in the endpoints. Convexity is preserved, monotonicity may be not, but the intersection point still is unique, and a good approximation for it seems to be the intersection of the $L(x)(x-a)(x-b)$-tangent in $b$ with the $R(x)(x-a)(x-b)$-tangent in $a$. Since:

$$L(x)(x-a)(x-b)= -\sum_{i\in I^-,\,A_i\neq -b}\frac{z_i(x-b)(x-a)}{x+A_i}-z_b(x-a),$$ $$R(x)(x-a)(x-b)= \sum_{i\in I^+,\,A_i\neq -a}\frac{z_i(x-b)(x-a)}{x+A_i}+z_a(x-b),$$ The approximation we get this way is the solution of $$z_b+(x-b)\sum_{i\in I^-,\,A_i\neq -b}\frac{z_i}{b+A_i}=z_a+(x-a)\sum_{i\in I^+,\,A_i\neq -a}\frac{z_i}{a+A_i},$$ i.e. $$\tilde{w}=\frac{(z_a-z_b)+\sum_{i\in I^-,\,A_i\neq -b}\frac{b\,z_i}{b+A_i}-\sum_{i\in I^+,\,A_i\neq -a}\frac{a\,z_i}{a+A_i}}{\sum_{i\in I^-,\,A_i\neq -b}\frac{z_i}{b+A_i}-\sum_{i\in I^+,\,A_i\neq -a}\frac{z_i}{a+A_i}}.$$ If now $R(\tilde{w})-L(\tilde{w})$ is positive we replace $[a,b]$ with $[\tilde{w},b]$, otherwise we replace $[a,b]$ with $[a,\tilde{w}]$ and iterate the argument. After two steps, both the endpoints of the original interval are replaced, so we do not need to "normalize" $L(x)$ and $R(x)$ anymore, we are back to (*) and formulas, if we want, return to be good-looking and simple.

Another idea is to approximate $L(x)$ and $R(x)$ (or their regularized version) with homographic functions (Padé approximants) with a 3-points $\left(c,d,\frac{c+d}{2}\right)$ fit (with respect to four parameters), then solve $\tilde{L}(x)=\tilde{R}(x)$. This may or may not have stability issues, I need to investigate further.

Another important fact that I believe exploitable is that $L(x)$ and $R(x)$ are much more than convex on $(a,b)$. $L(x)$ is an analytic function with non-negative derivatives of any order, $R(x)$ is an analytic function with non-negative derivatives of even order and non-positive derivatives of odd order over $(a,b)$. The same holds for $1/R(x)$ and $1/L(x)$, that are also non-negative and bounded over $[a,b]$. The intersection of tangents in the endpoints gives the approximation: $$\tilde{w} = \frac{R(a)L(b)(R(a)-L(b))+(bL'(b)R^2(a)-aR'(a)L^2(b))}{R^2(a)L'(b)-L^2(b)R'(a)}.$$

This is an old idea, took from simultaneous roots approximation algorithms, but may be really effective. Compose $f$ with a linear map such that $[a,b]=[-1,1]$. Consider the map $h$ that sends $z$ to $\frac{1}{2}\left(z+\frac{1}{z}\right)$. For any $r>0$, $h^{(n)}(r)$ converges pretty fast to $1$, and for any $r<0$, $h^{(n)}(r)$ converges pretty fast to $-1$. Moreover, $h$ is easy to invert. So, by composing with $h$ or its inverse function, we have the choice to "condensate" the singularities of $f,L,R$ near the endpoints of $[-1,1]$ or to "push them away". In both cases, it is very likely to have the possibility to provide extremely tight bounds for $L(x)$ and $R(x)$, so for $\tilde{w}$.

  • $\begingroup$ You discovered by yourself (your first equation) that this is the so-called Rachford-Rice equation which was published in 1952. As I said in the post, we just look for the root located between the first left and the first right asymptotes. In 1993, I found that there are tighter bounds which make the problem better conditioned. However, the equation is not very well conditioned because the K(i) can range between epsilon and infinity. My goal is to reduce the number of iterations since this equation has to be solved billions of times in a single run. $\endgroup$ – Claude Leibovici Jan 11 '14 at 6:08
  • 1
    $\begingroup$ Your answer is not naive at all and I shall work from here. As you saw, I spent a large part of my life working this problem (defining bounds inside the asymptotes, fonction transforms to remove singularities, convex transforms, ...). I must say that I do not know where your idea will take me but, for first time in years, I see some possibly interesting perspectives. Thanks. If,by chance, you find anything new, please let me know. $\endgroup$ – Claude Leibovici Jan 12 '14 at 7:11
  • 1
    $\begingroup$ As you understood, I dedicated a lot of time to this problem. You came with new ideas which, I am sure, are worth to be explored and investigated. Would you like we continue together ? Cheers. $\endgroup$ – Claude Leibovici Jan 13 '14 at 9:47
  • 1
    $\begingroup$ Glad to hear that ! This "simple" problem is of major importance in chemical engineering. Any new idea for tighter bounds "easy" to obtain could have a significant impact. I think that, over my life time, I probably spent more that five years full time working it ! Cheers. $\endgroup$ – Claude Leibovici Jan 13 '14 at 11:50
  • 1
    $\begingroup$ The number of asymptotes is equal to the number of components in the mixture and can be on the order of thousands. Z(i) are irrational numbers (mole fractions). Please forget any approach based on Fourier series. $\endgroup$ – Claude Leibovici Jan 14 '14 at 5:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.