• 0 Vote(s) - 0 Average
• 1
• 2
• 3
• 4
• 5
 Find all fixed points of exp[b] mike3 Long Time Fellow Posts: 368 Threads: 44 Joined: Sep 2009 09/29/2014, 01:49 AM This problem -- how to compute the various branches of $W(z)$ -- is already a solved problem. See here: http://www.aqit.uwo.ca/~rcorless/frames/...index.html http://www.aqit.uwo.ca/~rcorless/frames/...ambertW.ps The following is a simple Pari/GP code snippet I came up with that will compute any branch of $W(z)$. It is based on the asymptotic from the above paper and others found online and guessed, as well as some old Lambert W code of mine. Hopefully it should work for all values of $z$: Code:```/* Computes all branches of the Lambert W function. W(z) is the principal branch W_0(z). Wbr(k, z) is any branch W_k(z). */ stirling1(n, k) = ((2*n-k)!)/((k-1)!) * sum(kk=0,n-k,1/((n+kk)*(n-k-kk)!*(n-k+kk)!) * sum(j=0,kk,((((-1)^j)*(j^(n-k+kk)))/(j!*(kk-j)!)))) c(k, m) = 1/m! * ((-1)^k) * stirling1(k+m, k+1) brlog(k, z) = log(z) + 2*k*Pi*I; Wapprox(k, z) = brlog(k, z) - log(brlog(k, z)) + sum(kk=0,5,sum(m=1,5,c(kk, m)*log(brlog(k, z))^m * brlog(k, z)^(-kk-m))) W0approx(z) = {    local(r);    local(e = exp(1));    local(oneovere = exp(-1));    /* Asymptotic expansion at z = -1/e */    if(abs(z + oneovere) < 0.6,       r = -1 + sqrt(2*e)*sqrt(z + oneovere) - (2*e)/3*(z + oneovere);      );    /* Crude fit for modest values of z */    if((abs(z + oneovere) >= 0.6) && (abs(z + oneovere) < 10),       r = log(z + oneovere);      );    /* Large asymptotic for large values of z */    if(abs(z + oneovere) >= 10,       r = log(z) - log(log(z)) + (log(log(z))/log(z));      );    return(r); } Wbr(k, z) = {    local(r);    local(expr);    if(x == 0, return(0));    if((k == 0) && (z == -1/exp(1)), return(-1.0));    /* Create initial guess. */    if(k == 0,       r = W0approx(z),       r = Wapprox(k, z)      );    /* Use Newton's method to refine to full precision */    rold = r + 1;    while(abs(rold - r) > 10^-(precision(0.01)-1),          expr = exp(r);          rold = r;          r = r - (r*expr - z)/(expr + r*expr);         );    /* Do one last iteration for insurance */    expr = exp(r);    r = r - (r*expr - z)/(expr + r*expr);    /* Done! */    return(r); } W(z) = Wbr(0, z);``` « Next Oldest | Next Newest »

 Messages In This Thread Find all fixed points of exp[b] - by MorgothV8 - 08/30/2014, 01:41 PM RE: Find all fixed points of exp[b] - by MorgothV8 - 09/17/2014, 04:17 PM RE: Find all fixed points of exp[b] - by jaydfox - 09/17/2014, 11:10 PM RE: Find all fixed points of exp[b] - by MorgothV8 - 09/18/2014, 04:31 AM RE: Find all fixed points of exp[b] - by MorgothV8 - 09/24/2014, 03:53 PM RE: Find all fixed points of exp[b] - by jaydfox - 09/26/2014, 04:47 PM RE: Find all fixed points of exp[b] - by mike3 - 09/29/2014, 01:49 AM RE: Find all fixed points of exp[b] - by mike3 - 09/29/2014, 02:41 PM RE: Find all fixed points of exp[b] - by mike3 - 10/07/2014, 01:32 AM RE: Find all fixed points of exp[b] - by mike3 - 10/07/2014, 01:50 AM RE: Find all fixed points of exp[b] - by Gottfried - 10/07/2014, 11:00 AM

 Possibly Related Threads... Thread Author Replies Views Last Post fixed point formula sheldonison 6 9,909 05/23/2015, 04:32 AM Last Post: mike3 Attempt to find a limit point but each step needs doubling the precision... Gottfried 15 16,858 11/09/2014, 10:25 PM Last Post: tommy1729 sexp(strip) is winding around the fixed points Kouznetsov 8 12,182 06/29/2009, 10:05 AM Last Post: bo198214 An error estimate for fixed point computation of b^x bo198214 0 2,444 05/31/2008, 04:11 PM Last Post: bo198214

Users browsing this thread: 1 Guest(s)