Función de Lambert aproximación$W_0$ rama

Question

Estoy buscando una aproximación simple, económica y muy precisa de la función de Lambert ($W_0$ rama) ($-1/e < x < 0$).

Answer 1

Aquí hay un post que hice en el sci.matemáticas hace un tiempo, con respecto a un método que yo he utilizado.

Análisis de $we^w$

Para $w\gt0$, $we^w$ aumenta monótonamente de $0$ a $\infty$. Cuando $w\lt0$, $we^w$ es negativo.

Por lo tanto, para $x\gt0$, $\mathrm{W}(x)$ es positivo y bien definido y aumenta monótonamente.

Para $w\lt0$, $we^w$ llega a un mínimo de $-1/e$ a $w=-1$. En $(-1,0)$, $w e^w$ aumenta monótonamente de $-1/e$ a $0$. En $(-\infty,-1)$, $w e^w$ disminuye monótona de $0$ a $-1/e$. Por lo tanto, en $(-1/e,0)$, $\mathrm{W}(x)$ puede tener uno de dos valores, uno en $(-1,0)$ y otro en $(-\infty,-1)$. El valor en $(-1,0)$ es llamado el principal valor.

La iteración

Utilizando el método de Newton para resolver $we^w$ produce el siguiente proceso iterativo paso para encontrar $\mathrm{W}(x)$: $$w_{\text{nuevo}}=\frac{xe^{w}+w^2}{w} {+1}$$

Los valores iniciales de $w$

Para el principal valor, al $-1/e\le x\lt0$, y al $0\le x\le10$, el uso de $w=0$. Al $x\gt10$, el uso de $w=\log(x)-\log(\log(x))$.

Para el no-valor principal, al $x$ es de $[-1/e,-.1]$, el uso de $w=-2$; y si $x$ es de $(-.1,0)$, el uso de $w=\log(-x)-\log(-\log(-x))$.

Esto nos dice que para la rama que desea, utilice la iteración con un valor inicial de $w=0$.

Answer 2

Corless, et al. De 1996, es una buena referencia para la aplicación de los Lambert $W$. Una de las claves para ser eficiente es una buena estimación inicial para el refinamiento iterativo. Aquí es donde la serie de aproximaciones son útiles. Para $(-1/\text{e}, 0)$ es muy útil el uso de dos diferentes expansiones: acerca de $0$ y sobre el punto de ramificación de la $-1/\text{e}$. De lo contrario, se podrían perder un montón de tiempo de ejecución de las iteraciones.

He aquí algunos unvectorized (para mayor claridad) código de Matlab para $W_0$. He indicado la ecuación números de Corless, et al. en el que varias líneas de base.

function w = fastW0(x)
% Lambert W function, upper branch, W_0
% Only valid for real-valued scalar x in [-1/e, 0]

if x > -exp(-sqrt(2))
    % Series about 0, Eq. 3.1
    w = ((((125*x-64)*x+36)*x/24-1)*x+1)*x;
else
    % Series about branch point, -1/e, Eq. 4.22
    p = sqrt(2*(exp(1)*x+1));
    w = p*(1-(1/3)*p*(1-(11/24)*p*(1-(86/165)*p*(1-(769/1376)*p))))-1;
end

tol = eps(class(x));
maxiter = 4;
for i = 1:maxiter
    ew = exp(w);
    we = w*ew;
    r = we-x; % Residual, Eq. 5.2
    if abs(r) < tol
        return;
    end

    % Halley's method, for Lambert W from Corless, et al. 1996, Eq. 5.9
    w = w-r/(we+ew-0.5*r*(w+2)/(w+1));
end

Aquí está una parcela de doble precisión, tolerancia a la máquina epsilon) de número de iteraciones hasta que la convergencia si no se utiliza ningún inicial de la serie, supongo, la serie acerca de la $0$, la serie sobre el punto de rama, o de ambas series:

Por supuesto, la evaluación de las series de aproximación tiene un costo en relación a la realización de una iteración. Uno puede ajustar el número de términos que se utilizan con el fin de sintonizar este.

Answer 3

El arce produce $$LambertW (x) = x- {x} ^ {2} + {\ frac {3} {2}} {x} ^ {3} - {\ frac {8} {3}} {x} ^ {4} + {\ frac {125} {24}} {x} ^ {5} + O \ left ({x} ^ {6} \ right), x \ a 0.$$ Comparemos$LambertW(-.2)=-.2591711018$ con $-0.2-(-0.2)^{2}+3/2(-0.2)^{3}-8/3(-0.2)^{4}+{\frac {125}{24}}(-0.2)^{5}= -.2579333334$.

Answer 4

Un procedimiento iterativo, dado en wikipedia, q & d-traducido a Pari / GP, que se adapta bien a mis necesidades:

  LW(x, prec=1E-80, maxiters=200) = local(w, we, w1e); 
   w=0;    
   for(i=1,maxiters,  
           we=w*exp(w);
           w1e=(w+1)*exp(w); 
           if(prec>abs((x-we)/w1e),return(w));
           w = w-(we-x)/(w1e-(w+2)*(we-x)/(2*w+2)) ;
        );
    print("W doesn't converge fast enough for ",x,"( we=",we);
    return(0);
 

Ejemplo:

  default(precision,200)
   x = -0.99999/exp(1)
 %382 = -0.367875762377

   y=LW(x)
 %383 = -0.995534517079

  [chk=y*exp(y);x;chk-x]
 %384 = 
[-0.367875762377]
[-0.367875762377]
[2.49762470622 E-163]