The positive solution of $a x e^{2x} + b x e^x - c = 0$ where $a,c > 0$

Question

I've stumbled upon the equation $$f(x) := a x e^{2x} + b x e^x - c = 0$$ where $a,c > 0$ while doing my research. I'm only interested in the positive solution, i.e. $x>0$. It is easy to show a unique positive solution exists since $f(0)=-c < 0$ and if $f(x_0) > 0$ then $f(x)>0, \forall x\geq x_0$ (obviously such $x_0$ exists). This shows the existence of the unique positive solution because of continuity.

I'm interested in the analytical expression of the solution (if such an expression exists). I tried to put this equation in one of the forms in https://arxiv.org/pdf/1902.08910.pdf without any luck.

Is there any known solution to this kind of equation or am I stuck with numerical solutions?

Answer 1

Since $a>0$, WLOG we can fix it as $a=1$ and
consider zeros of the function with two parameters instead of three, \begin{align} f(x) &= x \exp(2x) + \alpha x \exp(x) - \beta = 0 \tag{1} \end{align}
where $\alpha=\tfrac{b}{a}$, $\beta=\tfrac{c}{a}$.

Unfortunately, I don't think there is a closed form solution for (1), but the numeric can be just fine. As a heuristic for initial guess of the positive root of (1) the root $x_0$ of a closely related equation

\begin{align} f_2(x) &= x^2 \exp(2x) + \alpha x \exp(x) - \beta = 0 , \tag{2} \end{align}
which has a positive root, \begin{align} x_0&= \operatorname{W}\Big(-\tfrac12\,\alpha+\tfrac12\,\sqrt{\alpha^2+4\beta}\Big) , \end{align} can be used.

Here is an example python code, which uses Halley's method to iteratively find the approximation of the root of (1) for $\alpha=1.2$, $\beta=1.6$, which is approximately $x\approx 0.39903258116278510968$:

import decimal
decimal.getcontext().prec = 20
exp=decimal.getcontext().exp

from scipy.special import lambertw as W
from numpy.lib.scimath import sqrt as sqrt

a=1.2; b=1.6
alpha = decimal.Decimal(a)
beta = decimal.Decimal(b)
 
def f(x):
  return x*exp(2*x) + alpha*x*exp(x) - beta

def df(x):
  return exp(x)*((1+2*x)*exp(x)+alpha*(1+x))

def ddf(x):
  return exp(x)*(4*exp(x)*(1+x)+alpha*(2+x))

def F(x):
  fx=f(x)
  dfx=df(x)
  ddfx=ddf(x)
  return x-2*fx*dfx/(2*dfx**2-fx*ddfx)

def fx0():
  return W(1/2*(-a+sqrt(a*a+4*b))).real

x=decimal.Decimal(fx0()); print(x)
# 0.49006785880157982537

x=F(x); print(x)
x=F(x); print(x)
x=F(x); print(x)
x=F(x); print(x)

# 0.39948996502974891810
# 0.39903258122349558881
# 0.39903258116278510968
# 0.39903258116278510968

And this is a contour plot of solutions of (1) for $\alpha\in[-3,3]$, $\beta\in[0,6]$, $x=0.1,\dots,1.2$: