Lambert W Function (hp-42s)
|
05-16-2020, 04:07 PM
Post: #1
|
|||
|
|||
Lambert W Function (hp-42s)
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 1E-32 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. |
|||
05-17-2020, 07:56 AM
Post: #2
|
|||
|
|||
RE: Lambert W Function (hp-42s)
Why not use the built-in solver?
Code: 00 { 36-Byte Prgm } Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
05-17-2020, 08:09 AM
Post: #3
|
|||
|
|||
RE: Lambert W Function (hp-42s) | |||
05-17-2020, 08:15 AM
Post: #4
|
|||
|
|||
RE: Lambert W Function (hp-42s)
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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
05-17-2020, 09:29 AM
Post: #5
|
|||
|
|||
RE: Lambert W Function (hp-42s) | |||
05-17-2020, 12:12 PM
Post: #6
|
|||
|
|||
RE: Lambert W Function (hp-42s)
I don't own a hp-42s, line 21 should be 1E-10 or something like that in the hp-42s. My program can take advantage of the fact that the hp-42s 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.3828293673E-1 i 3.605924718714E-1 so i^i^i^i...... ≈ 0.43828293673+ i 0.3605924718714 |
|||
05-18-2020, 08:04 AM
Post: #7
|
|||
|
|||
RE: Lambert W Function (hp-42s)
(05-17-2020 09:29 AM)Gerald H Wrote: What's wrong with this? It has two global labels ;-) Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
05-18-2020, 03:54 PM
(This post was last modified: 05-22-2020 03:08 PM by Albert Chan.)
Post: #8
|
|||
|
|||
RE: Lambert W Function (hp-42s)
A better convergence test might be (x-prev_x)*eps + x - x = 0
With Newton's method, eps can be set to "half" the calculator precision. Example, 12 digits calculator, eps = 1e-6 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.3127784610672455-1.1245675977983851j) >>> g(_,a) →(1.3607126098119278-1.119058042805926j) >>> g(_,a) →(1.3606248500898628-1.1194391815927829j) >>> g(_,a) →(1.3606248702911174-1.1194391662423497j) >>> g(_,a) →(1.3606248702911177-1.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.56641733028546448-0.68845322710770207j) >>> log(y) → (0.56641733028546437-0.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 |
|||
05-18-2020, 10:51 PM
Post: #9
|
|||
|
|||
RE: Lambert W Function (hp-42s)
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.
|
|||
05-19-2020, 12:29 AM
(This post was last modified: 05-25-2020 12:14 PM by Albert Chan.)
Post: #10
|
|||
|
|||
RE: Lambert W Function (hp-42s)
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 |
|||
05-20-2020, 03:25 PM
(This post was last modified: 05-20-2020 05:37 PM by Albert Chan.)
Post: #11
|
|||
|
|||
RE: Lambert W Function (hp-42s)
(05-19-2020 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 |
|||
05-21-2020, 12:09 AM
Post: #12
|
|||
|
|||
RE: Lambert W Function (hp-42s)
Hello Albert, how very interesting your equation, and it looks simpler too.
|
|||
05-22-2020, 11:39 AM
(This post was last modified: 05-22-2020 11:41 AM by Werner.)
Post: #13
|
|||
|
|||
RE: Lambert W Function (hp-42s)
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 HP-35S, 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 1E-15 to 1E-5. 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 { 145-Byte Prgm } define FX: Code: >LBL "FX" "FX" 0 XEQ "CSLV" (0.43829..,0.36059..) XEQ "FX" results in (0,-1e-34) 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 { 115-Byte 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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
05-23-2020, 12:43 AM
(This post was last modified: 05-23-2020 12:48 AM by Albert Chan.)
Post: #14
|
|||
|
|||
RE: Lambert W Function (hp-42s)
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 (34-byte) 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 |
|||
05-23-2020, 04:20 AM
(This post was last modified: 05-23-2020 10:02 AM by Werner.)
Post: #15
|
|||
|
|||
RE: Lambert W Function (hp-42s)
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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
06-11-2020, 03:27 AM
(This post was last modified: 06-11-2020 07:58 AM by Albert Chan.)
Post: #16
|
|||
|
|||
RE: Lambert W Function (hp-42s)
This gives exp(W(a)), which equals a / W(a)
Code: 00 { 33-Byte Prgm } Example, for i^i^i ... -1 SQRT LN +/- XEQ "eW" 1/X → 0.438282936727 + 0.360592471871 j |
|||
06-11-2020, 05:17 AM
Post: #17
|
|||
|
|||
RE: Lambert W Function (hp-42s)
Hi Albert,
something must be missing? Your routine expects two arguments. (the '+' in line 3) Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
06-11-2020, 08:12 AM
(This post was last modified: 06-11-2020 08:14 AM by Albert Chan.)
Post: #18
|
|||
|
|||
RE: Lambert W Function (hp-42s)
Hi, Werner
Thanks. Code corrected. BTW, can the code be done all in the stack (no "a") ? |
|||
06-11-2020, 09:20 AM
Post: #19
|
|||
|
|||
RE: Lambert W Function (hp-42s)
Sure. and if you replace lines 10 and 11 by
X<>Y STO+ ST Z X<> ST Z it is even 41-compatible. Code: 00 { 33-Byte Prgm } 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
06-11-2020, 10:42 AM
Post: #20
|
|||
|
|||
RE: Lambert W Function (hp-42s)
(06-11-2020 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: 6 Guest(s)