LambertW function
05-17-2020, 08:48 AM (This post was last modified: 11-07-2020 09:39 AM by Stevetuc.)
Post: #1
 Stevetuc Senior Member Posts: 300 Joined: Jan 2014
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
11-07-2020, 08:41 AM (This post was last modified: 11-10-2020 08:26 AM by Stevetuc.)
Post: #2
 Stevetuc Senior Member Posts: 300 Joined: Jan 2014
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.
11-07-2020, 11:57 AM (This post was last modified: 11-07-2020 12:37 PM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
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 ...
11-07-2020, 01:50 PM (This post was last modified: 11-07-2020 01:50 PM by Stevetuc.)
Post: #4
 Stevetuc Senior Member Posts: 300 Joined: Jan 2014
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.
 « Next Oldest | Next Newest »

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