+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: HP49-50G Lambert function (/thread-15747.html)
HP49-50G Lambert function - Gil - 10-17-2020 08:42 PM
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
RE: HP49-50G Lambert function - Gil - 10-18-2020 10:22 AM
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
RE: HP49-50G Lambert function W(z) - Gil - 10-19-2020 09:06 AM
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. RC '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 CR SWAP DROP = 2. ARRY 'a' 'b' 2. ARRY [ 1. 1. ] MSLV OBJ DROP RC 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
RE: HP49-50G Lambert function - Gil - 10-19-2020 09:13 AM
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]
RE: HP49-50G Lambert function - Gil - 10-19-2020 09:40 AM
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
RE: HP49-50G Lambert function - Albert Chan - 10-20-2020 01:48 AM
(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-15004-post-136892.html#pid136892
http://www.finetune.co.jp/~lyuka/technote/lambertw/lambertw-42s.html
RE: HP49-50G Lambert function - Gil - 10-20-2020 08:28 AM
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
RE: HP49-50G Lambert function - Gil - 10-20-2020 08:42 AM
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
RE: HP49-50G Lambert function - Albert Chan - 10-20-2020 04:40 PM
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()) = 0RE: HP49-50G Lambert function - Gil - 10-20-2020 10:44 PM
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
RE: HP49-50G Lambert function - Albert Chan - 10-21-2020 02:54 AM
(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]](https://upload.wikimedia.org/wikipedia/commons/thumb/b/b9/Mplwp_lambert_W_branches.svg/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)
RE: HP49-50G Lambert function - Gil - 10-21-2020 11:44 AM
Brilliant!
RE: HP49-50G Lambert function - Albert Chan - 10-21-2020 12:13 PM
(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