Post Reply 
(28/48/50) Lambert W Function
04-02-2023, 11:12 PM
Post: #23
RE: (28/48/50) Lambert W Function
(03-26-2023 06:43 PM)Albert Chan Wrote:  Because solving for y = e^W(a) is easier, W code is not recommended.

W code is perfect for other branches! (expW code for W branch 0, -1)

Numbering the branches of the Lambert W function

[Image: unwindw.svg]

Except for 2*k*pi*I term, it is basically same formula for W code!

lua> f = fn'x: x + I.log(x) - I.log(a) - 2*k*pi*I'
lua> df = fn'x: 1 + 1/x'

Let X = |x|, A = |a|

log(x) - log(a) = log(X/A) + I*(arg(x) - arg(a))
log(X/A) = 1 + log(X*(r+err)/A) ≈ 1 + log1p((X*err - (A+r) + r*(1+X)) / A)

This means we can evaluate f accurately around branch point!

f=0 guaranteed we get the correct branch, which is nice.
There is no need to use complicated W guesses.

lua> a, k = 3+4*I, 5
lua> x = 2*k*pi*I -- guess for Wk(a)
lua> repeat h=f(x)/df(x); x=x-h; print(x, I.abs(h)) until x == x+1e-6*h

(-1.815554248884793+30.714634540472343*I)       1.9462907525577031
(-1.817005891456528+30.71333413897562*I)        0.0019489253984594729
(-1.8170058918466274+30.713334137004896*I)      2.0089622485451688e-009
(-1.8170058918466274+30.713334137004896*I)      0

>>> from mpmath import *
>>> lambertw(3+4j, k=5)
mpc(real='-1.8170058918466274', imag='30.713334137004896')
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
(28/48/50) Lambert W Function - John Keith - 03-20-2023, 08:43 PM
RE: (28/48/50) Lambert W Function - Albert Chan - 04-02-2023 11:12 PM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 11:04 AM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 02:47 PM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 06:46 PM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 09:50 PM
RE: (28/48/50) Lambert W Function - Gil - 01-30-2024, 12:33 AM
RE: (28/48/50) Lambert W Function - Gil - 01-30-2024, 12:04 PM
RE: (28/48/50) Lambert W Function - Gil - 01-30-2024, 02:52 PM
RE: (28/48/50) Lambert W Function - Gil - 01-31-2024, 07:10 PM



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