Lambert W Function (hp42s)

05162020, 04:07 PM
Post: #1




Lambert W Function (hp42s)
The Lambert W Function can be implemented in your calculator using Newton's Method. To find W(A), the equation to solve is:
f(x) = x*exp(x)A = 0 (1) f”(x) = (x+1)*exp(x) (2) A recursive formula for Newton's Method is: X = x+(A/exp(x)x)/(x+1)) (3) From (1) we have x*exp(x) = A or ln(x) + x = ln(A); x = ln(A)ln(x); x ~ ln(A); ln(A) is a good initial value for the Newton Method when A>1, ln(A+1.4) is a good initial value for all reals ≥ 1/e, for complex numbers we use ln(A). The next program implements the Lambert W Function For the hp 42s: 01 LBL “LWF” 02 1.4 03 X<>Y 04 ENTER 05 REAL? 06 RCL+ ST Z 07 LN 08 STO “X” 09 LBL 00 10 X<>Y 11 RCL ST X 12 RCL “X” 13 E↑X 14 ÷ 15 RCL “X” 16 1 17 RCL+ “X” 18 ÷ 19 STO+ “X” 20 ABS 21 1E32 22  23 X>0? 24 GTO 00 25 RCL “X” 26 END The program doesn't check if you entered a value where the function its not defined. Values of x < 1/e should be entered as complex. The strategy to solve equations with exponential using the Lambert W Function is to use algebraic manipulations to get an equation of the form x*exp(x) = a; then x = W(a) or x*ln(x) = a; then ln(x) = W(a) and x = exp(W(a)). Example: to solve x^x = π, we have ln(x^x) = ln(π), x*ln(x) = ln(π); x = exp(W(ln(π))); in the calculator: [■][π][LN][XEQ][LWF][■][e^x] give us a value of x ~ 1.85410596792. 

05172020, 07:56 AM
Post: #2




RE: Lambert W Function (hp42s)
Why not use the builtin solver?
Code: 00 { 36Byte Prgm } Cheers, Werner 

05172020, 08:09 AM
Post: #3




RE: Lambert W Function (hp42s)  
05172020, 08:15 AM
Post: #4




RE: Lambert W Function (hp42s)
I'm afraid Free42/DM42 has taken over for me ;)
(actually, the version that I have stored is with STO; but I wanted to see if it worked with LSTO, and of course it does. I will not doubt you again, Thomas.) Cheers, Werner 

05172020, 09:29 AM
Post: #5




RE: Lambert W Function (hp42s)  
05172020, 12:12 PM
Post: #6




RE: Lambert W Function (hp42s)
I don't own a hp42s, line 21 should be 1E10 or something like that in the hp42s. My program can take advantage of the fact that the hp42s can handle complex numbers. For example let's consider a infinite tower of the imaginary number i, i^i^i^i......
We have i^i^i^i...... = x; i^(i^i^i^i)...... = x; i^x = x or i = x^(1/x) taking ln on both sides of the equation ln(i) = 1/x ln(x); ln(i) = 1/x ln(x); ln(i) = 1/x ln(1/x), calling W to the Lambert W function. W(ln(i)) = W(1/x ln(1/x)); W(ln(i)) = ln(1/x); 1/x = exp(W(ln(i))) and finally x=exp(W(ln(i))). In the calculator: 0 [ENTER] 1 [COMPLEX] [LN] [+/] [XEQ] LWF [+/][■][E↑X] You get 4.3828293673E1 i 3.605924718714E1 so i^i^i^i...... ≈ 0.43828293673+ i 0.3605924718714 

05182020, 08:04 AM
Post: #7




RE: Lambert W Function (hp42s)
(05172020 09:29 AM)Gerald H Wrote: What's wrong with this? It has two global labels ;) Werner 

05182020, 03:54 PM
(This post was last modified: 05222020 03:08 PM by Albert Chan.)
Post: #8




RE: Lambert W Function (hp42s)
A better convergence test might be (xprev_x)*eps + x  x = 0
With Newton's method, eps can be set to "half" the calculator precision. Example, 12 digits calculator, eps = 1e6 is sufficient. Quote:For example let's consider a infinite tower of the imaginary number i, i^i^i^i...... Instead of using LambertW, we may also do the direct route. x = k^k^k ... = k^x log(x) = x log(k) Let f = log(x)  x log(k), f' = 1/x  log(k) Newton's iteration x = x  f(x)/f'(x): x = (1  log(x)) / (1/x  log(k) y = 1/x = (y  log(k)) / (log(y) + 1) >>> from cmath import * >>> a = log(1j) # k = 1j >>> g = lambda y,a: (y+a)/(log(y)+1) # g(1,a) = 1+a >>> g(1+a,a) →(1.31277846106724551.1245675977983851j) >>> g(_,a) →(1.36071260981192781.119058042805926j) >>> g(_,a) →(1.36062485008986281.1194391815927829j) >>> g(_,a) →(1.36062487029111741.1194391662423497j) >>> g(_,a) →(1.36062487029111771.11943916624235j), pass convergence test, eps=2**26 >>> y = _ >>> x = 1/y →(0.43828293672703211+0.36059247187138549j) >>> 1j ** x →(0.43828293672703211+0.36059247187138554j) x = 1/y = W(a)/a → (a/y) exp(a/y) = a → exp(a/y) = y → W(a) = a/y = log(y) >>> a/y → (0.566417330285464480.68845322710770207j) >>> log(y) → (0.566417330285464370.68845322710770218j) Comment: Solving for x = k^x is not quite the same as x = k^k^k ... Example: 3 = (³√3)^3, but (³√3) ^ (³√3) ^ (³√3) ... ≈ 2.47805268029 ≠ 3 

05182020, 10:51 PM
Post: #9




RE: Lambert W Function (hp42s)
Hello Albert, I know that any equation that can be solved by Lambert W Function, can be solved using a numeric method, but I found satisfying using the algebraic manipulation so I can use Lambert W Function.


05192020, 12:29 AM
(This post was last modified: 05252020 12:14 PM by Albert Chan.)
Post: #10




RE: Lambert W Function (hp42s)
Hi, Juan14
I was experimenting another way to get W(a), via infinite tetration route. Solving for y = exp(W(a)) maybe better than directly going for W(a): Code: from cmath import * >>> test(1) y>w = ((0.56714329040978384+0j), 4) w = ((0.56714329040978384+0j), 5) >>> test(1e10) y>w = ((20.028685413304952+0j), 5) w = ((20.028685413304952+0j), 8) >>> test(1e100) y>w = ((224.84310644511851+0j), 4) w = ((224.84310644511851+0j), 10) Note: both versions use the same guess, y = e^x = 1+a Note: y>w uses a/y instead of log(y) to recover W(a), to reduce errors when y ≈ 1 

05202020, 03:25 PM
(This post was last modified: 05202020 05:37 PM by Albert Chan.)
Post: #11




RE: Lambert W Function (hp42s)
(05192020 12:29 AM)Albert Chan Wrote: g = lambda y,a: (y+a)/(log(y)+1) Both Newton's iteration formula are equivalent, but going for y converge faster. (try plotting x*exp(x) vs x*ln(x), 2nd curve slope = 1+ln(x), does not vary as much) W(a) = x x * exp(x) = a Let y = exp(x), we have y * ln(y) = a Newton's iteration: y ← y  (y*ln(y)a) / (ln(y)+1) = (y+a) / (ln(y)+1) When y converged, y = exp(x) = a/x, we have W(a) = ln(y) = a/y 

05212020, 12:09 AM
Post: #12




RE: Lambert W Function (hp42s)
Hello Albert, how very interesting your equation, and it looks simpler too.


05222020, 11:39 AM
(This post was last modified: 05222020 11:41 AM by Werner.)
Post: #13




RE: Lambert W Function (hp42s)
Ah yes, my simple program of course only works for reals.
Okay then, let's use a general purpose complex solver. I converted the one written by Valentin Albillo for the HP35S, to Free42/DM42. (you can find it here on his site) Of course, it is equally usable on the 42S if you change all LSTO statements to STO, and change the 1E15 to 1E5. I changed three things in the program:  a more accurate version of the root of the quadratic  a simpler test to see whether the number is sufficiently close to a real (that is not possible on a 35S or Valentin would've used it)  and just taking the real part of the current estimate if the test is positive Usage: ALPHA: name of function to solve X: guess XEQ "CSLV" The guess may be real or complex, and the resulting root as well. Code: 00 { 145Byte Prgm } define FX: Code: >LBL "FX" "FX" 0 XEQ "CSLV" (0.43829..,0.36059..) XEQ "FX" results in (0,1e34) Implement Lambert's W: Code: >LBL "W" x = i^x x = e^(x*ln(i)) x*ln(i) = ln(i)*e^(x*ln(i)) ln(i) = x*ln(i)*e^(x*ln(i)) Let y = x*ln(i), then solve y^*e^y = ln(i) or calculate W(ln(i)) 1 SQRT LN +/ XEQ "W" 1 SQRT LN +/ / results in the almost the the same number as above. However, Houston, we have a problem here. The definition of WX must not use the local variables used by CSLV, namely X, S, U and V. If we had used variable X instead of W, the routine would not have worked. To remove this dependency, here's a version of CSLV that does not use alphanumeric variables, except REGS: (making it less usable on a real 42S because it will overwrite REGS) 00 { 115Byte Prgm } 01▸LBL 01 02 RCL 00 03 + 04 ASTO ST L 05 GTO IND ST L 06▸LBL "CSLV" 07 2 08 ENTER 09 NEWMAT 10 ENTER 11 COMPLEX 12 + 13 LSTO "REGS" 14 1ᴇ15 15 STO 01 16▸LBL 02 17 CLX 18 XEQ 01 19 STO+ ST X 20 STO 02 21 RCL 01 22 +/ 23 XEQ 01 24 STO 03 25 RCL 01 26 XEQ 01 27 ENTER 28 RCL+ 03 29 RCL 02 30  31 RCL 01 32 X↑2 33 ÷ 34 X<>Y 35 RCL 03 36  37 RCL 01 38 STO+ ST X 39 ÷ 40 ÷ 41 RCL 02 42 LASTX 43 ÷ 44 STO ST Z 45 × 46 1 47  48 +/ 49 SQRT 50 1 51 + 52 ÷ 53 STO 00 54 RCL 00 55 ÷ 56 ABS 57 RCL 01 58 ABS 59 X↑2 60 X<Y? 61 GTO 02 62 RCL 00 63 SIGN 64 COMPLEX 65 X<>Y 66 R↓ 67 ABS 68 X>Y? 69 GTO 00 70 RCL 00 71 COMPLEX 72 R↓ 73 STO 00 74▸LBL 00 75 RCL 00 76 END Hope you like it, Werner 

05232020, 12:43 AM
(This post was last modified: 05232020 12:48 AM by Albert Chan.)
Post: #14




RE: Lambert W Function (hp42s)
Hi, Werner
Very nice complex solver. Regarding LambertW code, it missed LN before entering solver. Instead of guess of LN(1.4+X) or LN(X), I found that guess LN(1+X) is simpler. Here is the revised LambertW (34byte) Code: 01▸LBL "W" Example, i^i^i ... = W(ln(i))/(ln(i)) = exp(W(ln(i))) 1 SQRT LN +/ XEQ "W" +/ EXP → 0.4382829367270321116269751635512648 + 0.3605924718713854859529405269060007i 

05232020, 04:20 AM
(This post was last modified: 05232020 10:02 AM by Werner.)
Post: #15




RE: Lambert W Function (hp42s)
Thanks Albert. You can then simply use LN1+X ;)
Update: no you can't! LN accepts complex arguments, but LN1+X does not. Another example of the incredible attention to detail in Free42/DM42. Cheers, Werner 

06112020, 03:27 AM
(This post was last modified: 06112020 07:58 AM by Albert Chan.)
Post: #16




RE: Lambert W Function (hp42s)
This gives exp(W(a)), which equals a / W(a)
Code: 00 { 33Byte Prgm } Example, for i^i^i ... 1 SQRT LN +/ XEQ "eW" 1/X → 0.438282936727 + 0.360592471871 j 

06112020, 05:17 AM
Post: #17




RE: Lambert W Function (hp42s)
Hi Albert,
something must be missing? Your routine expects two arguments. (the '+' in line 3) Werner 

06112020, 08:12 AM
(This post was last modified: 06112020 08:14 AM by Albert Chan.)
Post: #18




RE: Lambert W Function (hp42s)
Hi, Werner
Thanks. Code corrected. BTW, can the code be done all in the stack (no "a") ? 

06112020, 09:20 AM
Post: #19




RE: Lambert W Function (hp42s)
Sure. and if you replace lines 10 and 11 by
X<>Y STO+ ST Z X<> ST Z it is even 41compatible. Code: 00 { 33Byte Prgm } 

06112020, 10:42 AM
Post: #20




RE: Lambert W Function (hp42s)
(06112020 09:20 AM)Werner Wrote: Thanks ! This is more than I expected Now, we can get W both ways, directly from the stack 1e100 XEQ "eW" ; eW = 4.44754573894e97 LN ; W = LN(eW) = 224.843106445 R↓ ; a = 1e100 ÷ ; W = a / eW = 224.843106445 Or, to check its accuracy, from definition W eW  a = 0 1e100 XEQ "eW" LN × − → 0 

« Next Oldest  Next Newest »

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