Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Find all fixed points of exp[b]
#7
This problem -- how to compute the various branches of -- 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 . 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 :

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);
Reply


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 jaydfox - 09/17/2014, 11:10 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

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)