numerical solve equation with a list of parameter

Question

I'm trying to numerically solve a equation like $$\exp(-A/x) + \exp(-B/x) = 1$$

there are two known lists of parameters $A=(A_1,A_2,...,A_n)$ and $B=(B_1,B_2,...,B_n)$,

I hope to find a list of $x=(x_1,x_2,...,x_n)$ for each $x_i$ relate to $(A_i,B_i)$, like this:

$$\exp(-A_i/x_i) + \exp(-B_i/x_i) = 1$$

I tried "vpasolve" which can not work with matrix parameters:

syms x;
eqn = exp(-pi$*$A/x) + exp(-pi$*$B/x) == 1;
sol = vpasolve(eqn,x)

syms x;
eqn = exp(-pi$*$A./x) + exp(-pi$*$B./x) == 1;
sol = vpasolve(eqn,x)

syms x;
x = zeros(size(A));
eqn = exp(-pi$*$A./x) + exp(-pi$*$B./x) == 1;
sol = vpasolve(eqn,x)

I've tried those ways, and all of them return an error.

although, the one-by-one for loop is working;

syms x;
for i=1:1:length(A)

eqn = exp(-pi$*$A(i)/x) + exp(-pi$*$B(i)/x) == 1;
sol = vpasolve(eqn,x);

end

and as Claude Leibovici suggested, by replace $y=1/x$ can make it 15~20% faster:

syms x;
S = zeros(length(A),1);

for i=1:1:length(A)

eqn = exp(-pi$*$A(i)$*$x) + exp(-pi$*$B(i)$*$x) == 1;
sol = vpasolve(eqn,x);
S(i) = 1/sol(1);

end

but it will take several hours to solve one group of parameters

Is there any method to solve it faster?

ps, like something solved as a matrix.

---

find this parallel for-loop, and it can run with 8/8 CPU cores at 100% occupation now. it saves a lot of time. still not sure if the "vpasolve" can work with lists A,B and x.

syms x;
S = zeros(length(A),1);

parfor i=1:length(A)

eqn = exp(-pi$*$A(i)$*$x) + exp(-pi$*$B(i)$*$x) == 1;
sol = vpasolve(eqn,x);
S(i) = 1/sol(1);

end

Answer 1

If the problem is to find the zero of function $$f(x)=e^{-\frac{A}{x}}+e^{-\frac{B}{x}}-1$$ more than likely vpasolve faces a problem of initial estimate.

Assuming that there is a solution (read @Lutz Lehmann's comment), first, let $x=\frac 1 y$ and now, try to solve $$g(y)=e^{-A y}+e^{-B y}-1$$ which is bounded by $\left(2e^{-A y}-1\right)$ and $\left(2e^{-B y}-1\right)$ which gives that the solution is somewhere between $a=\frac{\log (2)}{A}$ and $b=\frac{\log (2)}{B}$.

This means that the solution for $x$ is somewhere between $\frac A{\log (2)}$ and $\frac B{\log (2)}$

Just use the initial parameter somewhere in this range and no problem.

Edit

As I wrote in comments, solving for $y$ is probaly better. We have $$g(0)=1 \qquad\qquad g'(0)=-(A+B)\qquad\qquad g''(0)=A^2+B^2$$ So, by Darboux theorem, Newton method should converge without any overshoot of the solution.

Starting with $y_0=0$, the second iterate of Newton method is $$y_2=\frac{1}{A+B}+\frac{e \left(e^{-\frac{A}{A+B}}+e^{-\frac{B}{A+B}}-1\right)}{A e^{\frac{B}{A+B}}+B e^{\frac{A}{A+B}}}$$

Trying for $A=2$ and $B=3$, this gives $y_2=0.27336$ while the solution is $y=0.28120$.

Update (after @Lutz Lehmann's answer)

Suppose that we look instead for the zero of function $$h(y)=\log \left(e^{-A y}+e^{-B y}\right)$$ Using Newton method, we have $$y_0=0 \implies y_1=\frac{2 \log (2)}{A+B}\implies$$ $$ y_2=\frac{2 \log (2)}{A+B}+\frac{\left(2^{-\frac{2 A}{A+B}}+2^{-\frac{2 B}{A+B}}\right) \log \left(2^{-\frac{2 A}{A+B}}+2^{-\frac{2 B}{A+B}}\right)}{ 2^{-\frac{2 A}{A+B}}\, A+ 2^{-\frac{2 B}{A+B}}\,B}$$ Applied to the worked case $y_2=0.281199$ while the solution is $0.281200$.

Simpler, nicer and almost as accurate is the first iterate of Householder method $$y=\frac{ \log(2)}{A+B}\,\, \frac{2 k-\log (2)}{k-\log (2)}\qquad \text{where} \quad k=\left(\frac{A+B}{A-B}\right)^2$$ which, for the worked case, gives $0.281212$.

Another solution is the expand as a Taylor series $h(y)$ and use series reversion to obtain $$y=t+\frac{ (A-B)^2}{4 (A+B)}t^2+\frac{ (A-B)^4}{8 (A+B)^2}t^3+\frac{(A-B)^4 \left(13 A^2-34 A B+13 B^2\right)}{192 (A+B)^3}t^4+O\left(t^5\right)$$ where $t=\frac{2 \log(2)}{A+B}$.

For the worked example, this would give $y=0.2812004$ to be compared to the exact $y=0.2811996$.

Answer 2

Ensure that $A<B$. Then the equation becomes almost linear for $y=1/x$ with the transformation $$ e^{Ay}=1+e^{-(B-A)y}\implies 0<y\le \frac{\ln2}A \\~\\ Ay=\ln(1+e^{-(B-A)y})\approx \ln(2)-\frac{B-A}2y $$ This gives an initial guess value of $y_0=\frac{2\ln2}{A+B}$, for $A=2$, $B=3$ this is $y_0=0.27726$

Apply the Newton step formula to this \begin{align} y_{+1}=N(y)&=y-\frac{ Ay-\ln(1+e^{-(B-A)y}) }{ A-\frac{-(B-A)e^{-(B-A)y}}{1+e^{-(B-A)y}} } \\[.5em] &=y-\frac{ 1+e^{-(B-A)y} }{ A+Be^{-(B-A)y} }(Ay-\ln(1+e^{-(B-A)y})) \end{align}

As the plot shows, the Newton method converges rapidly from around the initial point.

Answer 3

Similar to other answers, let $y=1/x$ to get

$$e^{-Ay}+e^{-By}=1.$$

It was shown by Claude Leibovici through simple inequalities that:

$$\frac{\ln(2)}{\max(A,B)}\le y\le\frac{\ln(2)}{\min(A,B)}.$$

In this answer I improve upon these bounds by taking their geometric and harmonic means and provides another good initial estimate.

---

Substituting $y=\ln(2)/\sqrt{AB}$ results in:

\begin{align}e^{-Ay}+e^{-By}&=e^{-A\ln(2)/\sqrt{AB}}+e^{-B\ln(2)/\sqrt{AB}}\\&=e^{-\ln(2)\sqrt{\frac AB}}+e^{-\ln(2)\sqrt{\frac BA}}\\&=2^\alpha+2^{1/\alpha}\tag{$\alpha=-\sqrt{\frac AB}$}\\&\le1.\end{align}

Substituting $y=\ln(2)/[(A+B)/2]=2\ln(2)/(A+B)=\ln(4)/(A+B)$ results in:

\begin{align}e^{-Ay}+e^{-By}&=e^{-A\ln(4)/(A+B)}+e^{-B\ln(4)/(A+B)}\\&=e^{-\ln(4)/(1+B/A)}+e^{-\ln(4)/(1+A/B)}\\&=4^{-1/(1+\alpha)}+4^{-1/(1+1/\alpha)}\tag{$\alpha=\frac AB$}\\&\ge1.\end{align}

This proves the solution satisfies

$$\frac{2\ln(2)}{A+B}\le y\le\frac{\ln(2)}{\sqrt{AB}}$$

which is a tighter bound than what the other answers show.

---

A few examples further suggest a weighted harmonic mean is a very good initial approximation as well:

$$\boxed{y\approx\frac{7\ln(2)}{A+5\sqrt{AB}+B}}$$

In fact following the above procedure, one derives that this estimate approximately produces the following bounds:

$$0.9974\le e^{-Ay}+e^{-By}\le1.0113$$

which are pretty tight.

---

Let $t=2\ln(2)/(A+B)$ and

$$y_\mathrm{CL}=t+\frac{ (A-B)^2}{4 (A+B)}t^2+\frac{ (A-B)^4}{8 (A+B)^2}t^3+\frac{(A-B)^4 \left(13 A^2-34 A B+13 B^2\right)}{192 (A+B)^3}t^4+\mathcal O\left(t^5\right)$$

be Claude Leibovici's estimate.

With $A=2$ and $B=3$, this gives the following approximations:

\begin{align}\frac{\ln(2)}B&=0.2310\\\frac{2\ln(2)}{A+B}&=0.2773\\y_\mathrm{CL}&=0.2811996\\y&=\boxed{0.2812004}\\\frac{7\ln(2)}{A+5\sqrt{AB}+B}&=0.2813\\\frac{\ln(2)}{\sqrt{AB}}&=0.2830\\\frac{\ln(2)}A&=0.3466\end{align}

With $A=1$ and $B=0.01$, this gives the following approximations:

\begin{align}\frac{\ln(2)}B&=0.6931\\\frac{2\ln(2)}{A+B}&=1.3726\\y_\mathrm{CL}&=2.3522\\\frac{7\ln(2)}{A+5\sqrt{AB}+B}&=3.2133\\y&=\boxed{3.3987}\\\frac{\ln(2)}{\sqrt{AB}}&=6.9315\\\frac{\ln(2)}A&=69.3147\end{align}

Note that although I did not show the Householder's method from Claude Leibovici's answer here, it performed similarly in the above case.

Depending on your values of $A$ and $B$ and how close they are to each other, you may want to choose a different initial estimate to get better convergence. If you're not to sure, then the initial estimate derived here is guaranteed to be rather safe to use.

---

You can feed this initial estimate into vpasolve with its third init_param argument according to the documentation:

syms y;

eqn = log(exp(-A(i) * y) + exp(-B(i) * y)) == 0;
y_init = 7 * log(2) / (A(i) + 5 * sqrt(A(i) * B(i)) + B(i));

y_sol = vpasolve(eqn, y, y_init);
x_sol = 1 / y_sol