HP Forums

Full Version: HP49-50G Lambert function
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Lamber function for HP49-50G
W.Y=>X

Having the value Y,
find the X value such that (...X) × EXP (...X) = Y.

The program gives almost instantaneous answer,
but do no tamper on previously registered variable.

Here is the program code:

« "Lambert W(Y):
Having result Y,
W(Y) finds / gives X
such that
(X * EXP X) = Y
W(X * EXP X) =X

" DROP => Y
« X X TYPE 'X*EXP(X)-Y' 'X' 1 ROOT UNROT 6 ==
IF
THEN DROP 'X' PURGE
ELSE 'X' STO
END "W(" Y + ")" + =>TAG
»
»

Example of use:

What is the Lambert x value for y=(sqrt 2) /2?
1) 2 SQRT 2 / gives. 707
2) 2 SQRT 2 / W.Y=>X gives 0.4506

Application:

Solve x^2*e^x=2

Take square root of the above expression:
x*e^(x/2) = sqrt 2

Divide by 2:
(x/2)*e^(x/2) = (sqrt 2) / 2

That has a form of:
"(X..X)" * e ^"(X..X) " = (sqrt 2)/2

Then W ["(X..X)" * e ^"(X..X)"] = W [(sqrt 2)/2]
Or "X.. X" = W [(sqrt 2)/2]

In our case we have:
W[(x/2)*e^(x/2)] = W[(sqrt 2) / 2]
Or x/2 = W[(sqrt 2) / 2]
Or x = 2 * W[(sqrt 2) / 2] = 2 * 0.4506 = 0.9012

Remarks welcome.

Regards,
Gil Campart
Lambert function for HP49-50G
Version 1.2*
W.Y=>X

Having the value Y,
find the X value such that (...X) × EXP (...X) = Y.

The program gives almost instantaneous answer,
but do no tamper on previously registered variable.

* On that version 2, I slightly modified the code when testing possible existence of variable X: I included here the case of possible existing variable X being a program.

Here is the program code:

« "Lambert W(Y):
Having result Y,
W(Y) finds / gives X
such that
(X * EXP X) = Y
W(X * EXP X) =X

" DROP =>Y
«
IFERR 'X' 'X' RCL
THEN
ELSE
END 'X*EXP(X)-Y' 'X' 1 ROOT UNROT DUP2 SAME
IF
THEN DROP2 'X' PURGE
ELSE SWAP STO
END "W(" Y + ")" + =>TAG
»
»

The same calculation example may apply:

Example of use:

What is the Lambert x value for y=(sqrt 2) /2?
1) 2 SQRT 2 / gives. 707
2) 2 SQRT 2 / W.Y=>X gives 0.4506

Application:

Solve x^2*e^x=2

Take square root, on both sides, of the above expression:
x*e^(x/2) = sqrt 2

Divide both sides by 2:
(x/2)*e^(x/2) = (sqrt 2) / 2

The above result has a[/u] form of:
"(X..X)" * e ^"(X..X) " = (sqrt 2)/2

Then W ["(X..X)" * e ^"(X..X)"] = W [(sqrt 2)/2]
Or "X.. X" = W [(sqrt 2)/2]

In our case we have:
W[(x/2)*e^(x/2)] = W[(sqrt 2) / 2]
Or x/2 = W[(sqrt 2) / 2]
Or x = 2 * W[(sqrt 2) / 2] = 2 * 0.4506 = 0.9012

Remarks welcome.

Regards,
Gil Campart
Lamber function for HP49-50G
W.Z=>X

Version 2.
- That version allows now complex numbers Z as argument (and will give a complex number as answer.

- Having the Z complex number,
find the complex solution small z,
such that (...z) × EXP (...z) = Z.

- The program allows also, as argument, negative numbers < -1/e, in which case the answer will also be a complex number z.

- For the other cases, when entering
Y W(Z) ENTER,
we get the searched real value x=W(Y).

- The program gives almost instantaneous answer,
but do no tamper on previously registered variable.

- Here is the new program code:
« "Lambert W(Z):
Having result Y
Y real in (-879.3 infinity)
or Z complex (a,b)
W(Z) finds / gives X
such that
(X * EXP X) = Z
W(X * EXP X) =X

" DROP =>Y
«
IF Y TYPE 1. == NOT
THEN
IF Y -1. EXP NEG <
THEN Y 0. RC 'Y' STO
END
END
IFERR 'X' 'X' RCL
THEN
ELSE
END
IF Y TYPE 1. ==
THEN DUP2 SAME NOT
IF
THEN 'X' PURGE
END 'a*EXP(a)*COS(b)-b*EXP(a)*SIN(b)' Y OBJ DROP = 'a*EXP(a)*SIN(b)+b*EXP(a)*COS(b)' Y CR SWAP DROP = 2. ARRY 'a' 'b' 2. ARRY [ 1. 1. ] MSLV OBJ DROP RC UNROT DROP2
ELSE 'X*EXP(X)-Y' 'X' 1. ROOT
END UNROT DUP2 SAME
IF
THEN DROP2 'X' PURGE
ELSE SWAP STO
END "W(" Y + ")" + =>TAG
»
»
Example of use:

What is the Lambert x value for y=(sqrt 2) /2?
1) 2 SQRT 2 / gives. 707
2) 2 SQRT 2 / W.Y=>X gives 0.4506

Application:

Solve x^2*e^x=2

Take square root of the above expression:
x*e^(x/2) = sqrt 2

Divide by 2:
(x/2)*e^(x/2) = (sqrt 2) / 2

That has a form of:
"(X..X)" * e ^"(X..X) " = (sqrt 2)/2

Then W ["(X..X)" * e ^"(X..X)"] = W [(sqrt 2)/2]
Or "X.. X" = W [(sqrt 2)/2]

In our case we have:
W[(x/2)*e^(x/2)] = W[(sqrt 2) / 2]
Or x/2 = W[(sqrt 2) / 2]
Or x = 2 * W[(sqrt 2) / 2] = 2 * 0.4506 = 0.9012

Remarks welcome.

Regards,
Gil Campart
Precision

The program gives almost instantaneous answer for non-complex solutions ; as expected, however, longer times are required for the iterations when dealing with complex numbers solutions.

Regards,
Gil[/b]
As some characters were missing in the previous transcription, here is a new format of the exact identical code:

Always version 2

\<< "Lambert W(Z):
Having result Y
Y real in (-879.3 \oo)
or Z complex (a,b)
W(Z) finds / gives X
such that
(X * EXP X) = Z
W(X * EXP X) =X

" DROP \-> Y
\<<
IF Y TYPE 1. == NOT
THEN
IF Y -1. EXP NEG <
THEN Y 0. R\->C 'Y' STO
END
END
IFERR 'X' 'X' RCL
THEN
ELSE
END
IF Y TYPE 1. ==
THEN DUP2 SAME NOT
IF
THEN 'X' PURGE
END 'a*EXP(a)*COS(b)-b*EXP(a)*SIN(b)' Y OBJ\-> DROP = 'a*EXP(a)*SIN(b)+b*EXP(a)*COS(b)' Y C\->R SWAP DROP = 2. \->ARRY 'a' 'b' 2. \->ARRY [ 1. 1. ] MSLV OBJ\-> DROP R\->C UNROT DROP2
ELSE 'X*EXP(X)-Y' 'X' 1. ROOT
END UNROT DUP2 SAME
IF
THEN DROP2 'X' PURGE
ELSE SWAP STO
END "W(" Y + ")" + \->TAG
\>>
\>>

Regards,
Gil Campart
(10-19-2020 09:13 AM)Gil Wrote: [ -> ]The program gives almost instantaneous answer for non-complex solutions ; as expected,
however, longer times are required for the iterations when dealing with complex numbers solutions.

You can speed up convergence by solving y = e^W(x) instead.

W e^W = y * log(y) = x

f = y * log(y) - x, f' = log(y) + 1, then apply Newton's method.

f' is insensitive to guess y, and will converge very fast, even with bad guess.

Code:
from cmath import *
def eW(x, n=5, r=1/e, y=None):
    if not y: y = 0.3*(x+r) + sqrt(2*r*(x+r)) + r
    if y == r: y += 1e-10   # avoid 0 slope
    for i in range(n): print i,y; y = (y+x)/(log(y)+1)
    return y

>>> eW(-1+2j)
0 (0.912470160638+1.60208654199j)
1 (0.985261305638+1.5911210409j)
2 (0.985506513217+1.59184289853j)
3 (0.98550646367+1.59184283455j)
4 (0.98550646367+1.59184283455j)
(0.98550646366951244+1.5918428345475875j)

>>> eW(-1e10+2e10j)
0 (-2999932566.3+6000109109.25j)
1 (-452114322.526+1139389284.48j)
2 (-387942931.588+998560851.502j)
3 (-387765196.033+998129458.041j)
4 (-387765194.562+998129453.635j)
(-387765194.56226903+998129453.63483131j)

https://www.hpmuseum.org/forum/thread-15...#pid136892
http://www.finetune.co.jp/~lyuka/technot...w-42s.html
New version with minor changes not taking into account Chan's post that I just read — for which I thank him — and will study more carefully later.

Regards,
Gil Campart
New version with minor changes not taking into account Chan's post that I just read — for which I thank him — and will study more carefully later.

Regards,
Gil Campart
I tried solving simultaneous equation for w = a+b*1j, using Tk! Solver.
Note: w * exp(w) = xr + 1j*xi

xr = exp(a) * (a*cos(b)-b*sin(b))
xi = exp(a) * (a*sin(b)+b*cos(b))

Rectangular form does not converge well, unless we gives it a good guess, for both a and b.
Converted to polar form, convergence is easier, only required guess of a = b = 0

x = w exp(w) = exp(w.real)*abs(w) * exp(1j*(w.imag + phase(w)))

TK Solver Rule Sheet
Code:
c = xr*xr + xi*xi  "abs(x)**2
d = atan2(xi,xr)   "phase(x)
a*a + b*b = c / exp(2*a)
mod(b+atan2(b,a)-d, 2*pi()) = 0
Thanks for the detailed explanation.

For negative numbers, I get the same complex numbers as when using Wolfram Alpha, and that up to -190.
Beyond that number, I get different — but possible, correct — results.

My program
W(-190.,0.): (3.73813626409,-2.54403205829)
= result given by Wolfram.

But
W(-191.,0.): (3.0807478,-8.212851).
not equal W(-191_Wolfram Œ): (3.742525,2.544493).

Strange?

Could you have a check in your program for w(-191 + 0i)?

Thanks and regards.

By the way, when solving for
x^2 = 2 ^x,
I find only 2 real solutions : -. 766 and 2; no way to find by calculus the missing real solution x=4.

Could you help me?

Regards,
Gil Campart
(10-20-2020 10:44 PM)Gil Wrote: [ -> ]But
W(-191.,0.): (3.0807478,-8.212851).
not equal W(-191_Wolfram Œ): (3.742525,2.544493).

Strange?

Could you have a check in your program for w(-191 + 0i)?

There are many solutions to w*exp(w) = -191

>>> y = eW(-191)
0 (-56.8217567265+11.8431109072j)
1 (-35.6200433622+23.0040034786j)
2 (-34.9016087888+23.7270165869j)
3 (-34.9018254893+23.7292860973j)
4 (-34.9018255003+23.729286094j)
>>> log(y)
(3.742525902068365+2.5444934919969335j)

Besides speed, another benefit of solving for y = e^W:
Recovering W = log(y), |W.imag| ≤ pi

Quote:By the way, when solving for
x^2 = 2 ^x,
I find only 2 real solutions : -. 766 and 2; no way to find by calculus the missing real solution x=4.

Assume x is real, and x > 0:
2 log(x) = x log(2)
(1/x) log(1/x) = -log(2)/2 = a

>>> a = -log(2)/2
>>> 1/eW(a)
0 (0.499474906605+0j)
1 (0.500000902271+0j)
2 (0.500000000003+0j)
3 (0.5+0j)
4 (0.5+0j)
(2+0j)

For W-1, we try guess from below.
[Image: 320px-Mplwp_lambert_W_branches.svg.png]

>>> 1/eW(a, y=0.2)
0 0.2
1 (0.240506189867-0j)
2 (0.24956480343-0j)
3 (0.24999902269-0j)
4 (0.249999999995-0j)
(4.0000000000000018+0j)

Now, for x < 0:
2 log(-x) = x log(2)
(-1/x) log(-1/x) = log(2)/2 = -a

>>> -1/eW(-a)
0 (1.30724304932+0j)
1 (1.30435370354+0j)
2 (1.3043511789+0j)
3 (1.3043511789+0j)
4 (1.3043511789+0j)
(-0.76666469596212317+0j)
Brilliant!
(10-21-2020 02:54 AM)Albert Chan Wrote: [ -> ]>>> a = -log(2)/2
>>> 1/eW(a, y=0.2)
0 0.2
1 (0.240506189867-0j)
2 (0.24956480343-0j)
3 (0.24999902269-0j)
4 (0.249999999995-0j)
(4.0000000000000018+0j)

There is a reason I picked guess for y below -a (≈ 0.346574)
If guess is above -a, Newton's formula, y ← (y+a)/(log(y)+1), might converge to e^W0 (*)
If guess is below -a (but still positive), it always converge to e^W-1

see plot {(y + -ln(2)/2)/(ln(y)+1), y}, y = 0 .. 1/e

(*) it might still converge to e^W-1, with iterations involving complex numbers.
First iteration overshooted to y<0, then all iterations turned complex.
I tried it, the cut-off is guess ≈ 0.365414
Reference URL's