HP4950G Lambert function

10172020, 08:42 PM
Post: #1




HP4950G Lambert function
Lamber function for HP4950G
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 

10182020, 10:22 AM
Post: #2




RE: HP4950G Lambert function
Lambert function for HP4950G
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 

10192020, 09:06 AM
Post: #3




RE: HP4950G Lambert function W(z)
Lamber function for HP4950G
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 

10192020, 09:13 AM
Post: #4




RE: HP4950G Lambert function
Precision
The program gives almost instantaneous answer for noncomplex solutions ; as expected, however, longer times are required for the iterations when dealing with complex numbers solutions. Regards, Gil[/b] 

10192020, 09:40 AM
Post: #5




RE: HP4950G Lambert function
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 

10202020, 01:48 AM
(This post was last modified: 10212020 01:12 AM by Albert Chan.)
Post: #6




RE: HP4950G Lambert function
(10192020 09:13 AM)Gil Wrote: The program gives almost instantaneous answer for noncomplex solutions ; as expected, 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 * >>> 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/thread15...#pid136892 http://www.finetune.co.jp/~lyuka/technot...w42s.html 

10202020, 08:28 AM
Post: #7




RE: HP4950G Lambert function
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 

10202020, 08:42 AM
Post: #8




RE: HP4950G Lambert function
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 

10202020, 04:40 PM
(This post was last modified: 10212020 12:52 AM by Albert Chan.)
Post: #9




RE: HP4950G Lambert function
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 

10202020, 10:44 PM
Post: #10




RE: HP4950G Lambert function
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 

10212020, 02:54 AM
Post: #11




RE: HP4950G Lambert function
(10202020 10:44 PM)Gil Wrote: But 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 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.2405061898670j) 2 (0.249564803430j) 3 (0.249999022690j) 4 (0.2499999999950j) (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) 

10212020, 11:44 AM
Post: #12




RE: HP4950G Lambert function
Brilliant!


10212020, 12:13 PM
Post: #13




RE: HP4950G Lambert function
(10212020 02:54 AM)Albert Chan Wrote: >>> a = log(2)/2 There is a reason I picked guess for y below a (≈ 0.34657359) If guess is above a, Newton's formula, y ← (y+a)/(log(y)+1), might converge to e^W_{0} (*) If guess is below a (but still positive), it always converge to e^W_{1} see plot (y + log(2)/2)/(log(y)+1), y = 0 .. 0.5 (*) 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 cutoff is guess = a + 0.0188407 ... = 0.3654143 ... It is safer to pick a guess from the side without the "overshoot". 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)