HP49-50G Lambert function - Printable Version +- 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. 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 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()) = 0``` RE: 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. >>> 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