Post Reply 
LambertW function
05-17-2020, 08:48 AM (This post was last modified: 11-07-2020 09:39 AM by Stevetuc.)
Post: #1
LambertW function
This cas program uses fsolve to calc Lambert fn

Code:
#cas
lmb(x):=
fsolve(equal(w*e^w-x,0),w)
#end

Graph from wikipedia
   

If x <0 lmb(x) returns both principal and negative branch solutions :

lmb(-0.5/e) returns [−2.67834699002,−0.231960952987]

lmb(-1/e) returns −0.999999972234 (negative and principal branch converge at -1 when x=-1/e)

lmb(1) returns 0.56714329041

lmb(2) returns 0.852605502014
Find all posts by this user
Quote this message in a reply
11-07-2020, 08:41 AM (This post was last modified: 11-10-2020 08:26 AM by Stevetuc.)
Post: #2
RE: LambertW function
Simplified the code and change to fsolve:

Code:
#cas
fsolve(w*e^w=z,w)▶lmb(z);
#end
Result for lmb(-1.78) in home or cas:
8.92180498562ᴇ−2+1.62562367443*i

Code:
#cas
fsolve(w*e^w=z,w,0)▶lmb(z);
#end

Edit: add initial guess of 0 to avoid terminal screen in cas. This forces going direct to iterative solver rather than first trying and failing with bisection solver.
Edit: it doesn't fail with Bisectional solver. The cas terminal screen is just for information:
Quote:Solving by bisection with change of variable x=tan(t) and t=-1.57..1.57. Try fsolve(equation,x=guess) for iterative solver or fsolve(equation,x=xmin..xmax) for bisection.
Find all posts by this user
Quote this message in a reply
11-07-2020, 11:57 AM (This post was last modified: 11-07-2020 12:37 PM by Albert Chan.)
Post: #3
RE: LambertW function
(11-07-2020 08:41 AM)Stevetuc Wrote:  Edit: add initial guess of 0 to avoid terminal screen in cas.
This forces going direct to iterative solver rather than first trying and failing with bisection solver.

You might want to mention guess of 0 will iterate for W0(x), i.e. principle branch. (*)

Also, guess 0 is same as guess x, but wasted 1 Newton iteration.
We might as well use guess = x

Newton: w - (w*exp(w) - x) / (w*exp(w) + exp(w))

With guess 0, first iteration of w = 0 - (0 - x) / (0 + 1) = x


(*) Assumed W0 is not complex (x ≥ -1/e), see comment below.

Comment: For fsolve, some randomization of guess is going on.
With complex ON, if we fsolve again and again, we got different solutions.

XCas> fsolve(w*e^w = -1.78, w=0)

0.0892180498562+1.62562367443*i
-1.4781113814-7.66344321151*i
-3.68225172433+70.6337502365*i
-2.07259091944+13.9900896316*i
-2.9207293675-32.8981741284*i
...

Without randomization, fsolve should not even converge. (w will not flip to complex)
Maybe this is the reason guess randomization kick in ...
Find all posts by this user
Quote this message in a reply
11-07-2020, 01:50 PM (This post was last modified: 11-07-2020 01:50 PM by Stevetuc.)
Post: #4
RE: LambertW function
(11-07-2020 11:57 AM)Albert Chan Wrote:  
(11-07-2020 08:41 AM)Stevetuc Wrote:  Edit: add initial guess of 0 to avoid terminal screen in cas.
This forces going direct to iterative solver rather than first trying and failing with bisection solver.

You might want to mention guess of 0 will iterate for W0(x), i.e. principle branch. (*)

(*) Assumed W0 is not complex (x ≥ -1/e), see comment below.

Thanks for pointing that out. I want to retain the secondary branch result when 0 >x ≥ -1/e so I've removed the initial guess on my prime.
Find all posts by this user
Quote this message in a reply
Post Reply 




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