La función W de Lambert

Consideremos la función, f(x)=x·ex

f=@(x) x.*exp(x);
fplot(f,[-4,1])
xlabel('x')
ylabel('y')
axis([-4, 1, -.5, 1])
grid on
line([-1,-1],[-1/exp(1),-1/exp(1)],'marker','.','markersize',15,'color','k')
set(gca,'yTick',[-1/exp(1),0,0.5,1])
set(gca,'yTickLabel',{'-1/e','0','0.5','1'})

La función f(x) tiene un mínimo en x=-1, y el valor de este mínimo es -1/e juega un papel importante en la definición de la función W de Lambert.

La función W de Lambert

La función W de Lambert, W(x) es la inversa de la función f(x).

W(x)· e W(x) =x

Se representa gráficamente, intercambiando los ejes horizontal y vertical mediante el comando view.

f=@(x) x.*exp(x);
hold on
fplot(f,[-4 -1])
fplot(f,[-1,1])
hold off
xlabel('y')
ylabel('x')
axis([-4 1 -.5 1])
view([90,-90])
grid on
line([-1,-1],[-1/exp(1),-1/exp(1)],'marker','.','markersize',15,'color','k')
set(gca,'ytick',[-1/exp(1),0,0.5,1])
set(gca,'yticklabel',{'-1/e','0','0.5','1'})

La función W(x) está definida en el intervalo [-1/e, ∞), con dos ramas:

La función W(x) tiene un valor en el intervalo [0,∞) y dos valores en el intervalo [-1/e,0)

Valores significativos de W(x)

Como vemos en la gráfica, W(0)=0 y W(-1/e)=-1. Existen otros valores significativos derivados de que W(x) es la función inversa de f(x)=x·ex.

En la primera gráfica, al valor de abscisa 2 le corresponde 2·e2 en la ordenada. Luego, en la función inversa, a la abscisa x=2e2 la corresponde la ordenada W(2e2)=2. En general, W(x·ex)=x. Valores importantes son

Otros valores de de W(x)

Otros valores de W(x) hay que calcularlos de forma numérica, mediante el procedimiento de Newton u otros procedimientos que convergen más rápidamente como el de Halley.

El cálculo deW(x) para un valor dado de x requiere resolver la ecuación transcendente, f(w)=wew-x=0.

El procedimiento de Newton calcula la raíz de la ecuación mediante un proceso iterativo de la siguiente forma

w n+1 = w n f( w n ) f'( w n )

f'(w) es la derivada de f(w) respecto de w. f'(w)=ew+w·ew=(1+w)ew.

El procedimiento de Newton precisa un valor de partida w0 en el proceso de iteracción:

x=1;
f=@(w) w*exp(w)-x; %función
f_prima=@(w) (w+1)*exp(w); %derivada
w0=1; %valor inicial
while(1)
        w=w0-f(w0)/f_prima(w0);
        if abs((w-w0)/w)<1e-5 
            break
        end
        w0=w;
    end
disp([x,w]);

Corremos el script y en la ventana de comandos, aparece el valor W(1)=0.5671

>>   0.5671

Creamos una función denominada lambert_newton, tal que dado el valor de x y la rama 0 ó, -1, correspondientes a las funciones W0(x) o W-1(x), devuelve el valor de w.

function w =lambert_newton(rama,x)
    if (rama==0)
        w0=1;
    else
        w0=-2;
    end
    while(1)
        w=w0-(w0*exp(w0)-x)/((w0+1)*exp(w0));
        if abs((w-w0)/w)<1e-5 
            break
        end
        w0=w;
    end
end

Llamamos a la función lambert_newton para obtener una tabla de valores, x, W(x). Tomamos valores de x comprendidos entre -1/e e ∞ y calculamos el valor correspondiente a la primera rama W0(x)

x1=[2*exp(2),exp(1),2,1,1/exp(1),-1/4,-2/exp(1)^2,-1/exp(1)];
for x=x1
    w=lambert_newton(0,x);
    disp([x,w])
end
   14.7781    2.0000
    2.7183    1.0000
    2.0000    0.8526
    1.0000    0.5671
    0.3679    0.2785
   -0.2500   -0.3574
   -0.2707   -0.4064
   -0.3679   -1.0000

Llamamos a la función lambert_newton para obtener una tabla de valores, x, W(x). Tomamos valores de x comprendidos entre -1/e y cero y calculamos el valor correspondiente a la segunda rama W-1(x)

x2=[-1/4,-2/exp(1)^2,-1/exp(1)];
for x=x2
    w=lambert_newton(-1,x);
    disp([x,w])
end
   -0.2500   -2.1533
   -0.2707   -2.0000
   -0.3679   -1.0000

En vez de utilizar la función propia, lambert_newton prodríamos haber utilizado la función lambertw de MATLAB, para obtener los mismos resultados. Ahora la utilizaremos para representar gráficamente la función W(x)

hold on
x=-1/exp(1):0.01:1;
plot(x,lambertw(0,x))

x=-1/exp(1):0.01:0;
plot(x,lambertw(-1,x))
hold off

axis([-.5, 1, -4, 1])
line([-1/exp(1), -1/exp(1)],[-1, -1],'marker','.','markersize',18,'color','black')
set(gca,'xtick',[-1/exp(1), 0, 0.5, 1])
set(gca,'xticklabel',{'-1/e','0','0.5','1'})
set(gca,'ytick',-4:1)
grid on
xlabel('x')
ylabel('w')
title('Lambert W')
legend({'W_0','W_{-1}'},'location','southeast')

Ecuaciones transcendentes

Muchas ecuaciones que contienen exponenciales se pueden resolver en términos de la función W Lambert, agrupando los términos de la forma f(x)·exp(f(x)).

Ecuación x+ex=0

Multiplicamos los dos miembros de la ecuación -x=ex por e-x

-x·e-x=1

La solución exacta es -x=W0(1).

>> -lambertw(0,1)
ans =   -0.5671

Ecuación, xx=8

Resolvemos esta ecuación siguiendo los pasos

x x =8 xlnx=ln8 lnx· e lnx =ln8 z· e z =ln8

Como ln8>0 solamente hay una solución

z= W 0 ( ln8 ) x=exp( W 0 ( ln8 ) )

>> exp(lambertw(0,log(8)))
ans =    2.3884

La ley del desplazamiento de Wien

En la página titulada El cuerpo negro, calculamos el máximo de la función de la distribución de Planck expresada en términos de la longitud de onda, la denominada ley de desplazamiento de Wien

5( e x 1)x e x =0x= hc λkT

Una solución de esta ecuación es x=0. Buscamos otras soluciones de esta ecuación transcendente

(x5) e x =5 (x5) e (x5) =5 e 5 z e z =5 e 5 z= W 0 ( 5 e 5 ) x=5+ W 0 ( 5 e 5 )

Como -5e-5=-0.034 está en el intervalo que va desde -1/e=-0.3679 a cero, hay dos posibles soluciones

>> 5+lambertw(-1,-5*exp(-5))
ans =     0
>> 5+lambertw(0,-5*exp(-5))
ans =    4.9651

Se obtuvo este valor utilizando procedimientos numéricos

Ecuación x2=2x

Representamos las dos funciones y=x2, e y=2x

f=@(x) x.^2;
g=@(x) 2.^x;
hold on
fplot(f,[-2,5])
fplot(g,[-2,5]),
hold off
grid on
axis([-2,5,-0.5,20])
xlabel('x')
ylabel('y')
title('Ecuación 2^x=x^2')

Vemos que hay tres puntos de intersección, x=2, x=4 y otro próximo a x=-1.

Hallando la raíz cuadrada

x 2 = 2 x x=± 2 x/2

Tomamos la raíz positiva para calcular las abscisas de los dos primeros puntos de intersección, siguiendo los pasos

lnx= xln2 2 x=exp( xln2 2 ) x·exp( xln2 2 )=1 xln2 2 exp( xln2 2 )= ln2 2 z·exp(z)= ln2 2

Como -ln2/2=-0.3466, está en el intervalo que va desde -1/e=-0.3679 a cero, hay dos posibles soluciones.

z 1 = W 0 ( ln2 2 ) z 2 = W 1 ( ln2 2 ) x 1 = 2 ln2 W 0 ( ln2 2 )=2 x 2 = 2 ln2 W 1 ( ln2 2 )=4

Calculamos ahora la abscisa del tercer punto de intersección, para el caso negativo

x= 2 x/2 ln(x)= xln2 2 x=exp( xln2 2 ) x·exp( xln2 2 )=1 xln2 2 exp( xln2 2 )= ln2 2 z·exp(z)= ln2 2

Como ln2/2=0.3466, es mayor que cero, solamente es posible una solución.

z 3 = W 0 ( ln2 2 ) x 3 = 2 ln2 W 0 ( ln2 2 )=0.7667

Así pues, las raíces exactas de la ecuación x2-2x=0, son, 2, 4 y

2 ln2 W 0 ( ln2 2 )

>> -2*lambertw(0, -log(2)/2)/log(2)
ans =    2.0000
>> -2*lambertw(-1, -log(2)/2)/log(2)
ans =    4.0000
>> -2*lambertw(0,log(2)/2)/log(2)
ans =   -0.7667

Ecuación ab·x+c=d·x+f

Consideremos en general, la ecuación ab·x+c=d·x+f con a>0, b≠0, d≠0

Efectuamos el cambio de variable, -t=b·x+b·f/d y seguimos los pasos

bx=t bf d dx+f=d( t bf d b )+f= td b a bx+c =dx+f a t a cbf/d = td b b d a cbf/d =t a t b d a cbf/d =t e tlna

Multiplicando ambos miembros por lna

b d ( lna ) a cbf/d =t( lna ) e tlna

Para conseguir la forma u=w·exp(w) ó w=W(u)

tlna=W( b d ( lna ) a cbf/d )

Deshaciendo el cambio de variable

x= 1 blna W( b d ( lna ) a cbf/d ) f d

Sea la ecuación 2x=10x

a=2, b=1, c=0, d=10, f=0.

La solución exacta es

x= 1 ln2 W( ln2 10 )

Como ln2/10=-0.069 está comprendido entre -1/e y 0, hay dos posibles soluciones, tal como puede verse en la figura

f=@(x) 10*x;
g=@(x) 2.^x;
hold on
fplot(f,[0,6])
fplot(g,[0,6]),
hold off
grid on
xlabel('x')
ylabel('y')
title('Ecuación 2^x=10x')

>> -lambertw(0,-log(2)/10)/log(2
ans =    0.1078
>> -lambertw(-1,-log(2)/10)/log(2)
ans =    5.8770

Alcance de un proyectil bajo la acción de una fuerza de rozamiento proporcional a la velocidad

En el estudio del movimiento de un proyectil sometido a la acción de una fuerza de rozamiento proporcional a la velocidad, obtuvimos las ecuaciones de la posición (x, y) en función del tiempo t.

x= v 0x b ( 1exp(bt) ) y= 1 b ( g b + v 0y )( 1exp(bt) ) g b t

donde g es la aceleración de la gravedad, b es la constante de proporcionalidad de la fuerza de rozamiento, v0x·cosθ, y v0y·sinθ, son las componentes de la velocidad de disparo v0, con ángulo de tiro θ

El tiempo de vuelo es el tiempo que tarda el proyectil desde que es disparado hasta que llega al suelo y=0.

g b + v 0 sinθgt=( g b + v 0 sinθ )exp(bt) 1 g g b + v 0 sinθ t= e bt

Se trata de un caso particular de la ecuación cuya solución obtuvimos en el apartado anterior. Haciendo la sustitución

xt ae bb c0 d g g b + v 0 sinθ f1

El tiempo t de vuelo es exactamente

t= 1 b ( u+W( u e u ) ) u=1 b v 0 g sinθ

Dado el tiempo de vuelo t, calculamos el alcance x del proyectil que denominamos R

R= v 0 cosθ b ( 1exp( uW( u e u ) ) )

Teniendo en cuenta que z=W(z)exp(W(z)), con z=ueu, el alcance R es

R= v 0 cosθ b ( 1 W( u e u ) u )

En la página titulada Fuerza de rozamiento proporcional a la velocidad, calculamos utilizando procedimientos numéricos el valor aproximado del alcance R, obteniendo una gráfica similar del alcance del proyectil R en función del ángulo de tiro θ, fijada la velocidad de disparo v0 y la constante de proporcionalidad b de la fuerza de rozamiento

b=0.05; %parámetro
v0=60;  %velocidad incial
R=zeros(1,70);
i=0;
for angulo=(10:80)*pi/180;
    i=i+1;
    u=-1-b*v0*sin(angulo)/9.8;
    R(i)=(v0*cos(angulo)/b)*(1-lambertw(0,u*exp(u))/u); %alcance 
end
plot(10:80,R)
grid on
ylabel('Alcance (m)')
xlabel('\theta')
title('Tiro parabólico con rozamiento')

Ejemplos en el Curso Interactivo de Física

La descarga de un depósito de líquido a través de un tubo capilar

Referencias

Seán Stewart. A new elementary function for our curricula?. The Petroleum Institute, United Arab Emirates. Australian Senior mathematics 19(2)

Thomas P. Dence. A Brief Look into the Lambert W Function Applied Mathematics, 2013, 4, 887-892.

Edward W. Packel, David S. Yuen. Projectile Motion with Resistance and the Lambert W Function. VOL. 35, NO. 5, NOVEMBER 2004 THE COLLEGE MATHEMATICS JOURNAL