Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Find all fixed points of exp[b]
#1
Hi, maybe too easy but I have problems.
I waant to find set of all fixed points for complex z^x.
I know that I can do:
L = fixed_point(b^z) for all b in Complex plane.
I'm doing:
fix_point(b) := Lambert_W(a), where a = -ln(b)
Labmert_W(a) := Tania(t), where t = ln(a)-1
Tania(t):
4 different implementations (after citizendium), depending on t values (imag/real).

I can get fixed points quite OK, for:
b=e, b=sqrt(2), b=2.
I'm computing:
L=fix_point(b)
fL = f(L) = pow(b, L)

Then I check |L-fL| and get values < 1e-19, which is OK.

My only problem is:
I want ALL fixed points for ALL complex bases b
f(b) = set(L), set(L) = set of fixed points not just one.

Problem is that for b=sqrt(2) i get L=4, which is OK, but I want to get L=2 --> how to do that?
Anybody will help?

My goal is to have complex map of
f(z) = set { fixed points of pow(z) }
I mean I want multivalued mapping.


EDIT:
It behaved really bad for bases near e^(1/e), specially equal and little below.
This was due to log of values with negative real part and almost zero imag part - which returned really bad results (because log of negative real is cut line), but more itarations helped.
Code updated to "adjust" iterations untill they behave OK. So there is limit 100 iterations, but after each one I check added value (in Tania function), if value drops |s| < EPS (1e-16) then abort iterations and return Tania result.
This helped for e^(1/e) and its vicinity but....
there is NO convergence of Tania function (this is the one from citizendium) when b = PI/3 (I wanted some special PI related value, and less than e^(1/e) and I tested PI/3 and Tania is not converging)
Checked such values:
//z = C(-1.0, 0.);
//z = C(0., 0.);
//z = C(1., 0.);
//z = C(exp(1.), 0.);
//z = C(exp(2.), 0.);
//z = C(exp(-exp(1.)), 0.0);
//z = C(exp(1./exp(1.)), 0.0);
//z = C(exp(1./exp(1.)-0.00000001), 0.0);
//z = C(exp(1./exp(1.)+0.00000001), 0.0);
//z = C(sqrt(2.), 0.0);
//z = C(M_PI/3., 0.0);
Updated code


Attached Files
.zip   math4.zip (Size: 1.54 KB / Downloads: 121)
Fuji GSW690III
Nikon D3, Nikkors 14-24/2.8, 24/1.4, 35/2, 50/1.4, 85/1.4, 135/2, 80-200/2.8
Reply
#2
No one can help how to find all fixed points of log(b), for all complex bases b.

I know exp() have infinite numbe rof fixed points, but two of them are "main" and are the same as log(), I want to find them using some algorithm and elementary functions. I know there is no closed for for them, but I want to find algorithm that approximates them for all bases.

Any ideas?
My program ffrom 1st post inds fixed point - but only one, how to find another one?
Fuji GSW690III
Nikon D3, Nikkors 14-24/2.8, 24/1.4, 35/2, 50/1.4, 85/1.4, 135/2, 80-200/2.8
Reply
#3
(09/17/2014, 04:17 PM)MorgothV8 Wrote: No one can help how to find all fixed points of log(b), for all complex bases b.

I know exp() have infinite numbe rof fixed points, but two of them are "main" and are the same as log(), I want to find them using some algorithm and elementary functions. I know there is no closed for for them, but I want to find algorithm that approximates them for all bases.

Any ideas?
My program ffrom 1st post inds fixed point - but only one, how to find another one?

Well, unfortunately, I don't have c++ code for anything related to tetration. I don't know how easily you can work with PARI/gp or SAGE/python, but those are the languages I've been working with. Don't get me wrong; I like the speed of c++ (I've been working in c/c++ for over 20 years).

But when it comes to tetration and similar mathematical projects, I've been working primarily with "arbitrary precision" math lately, and that would require linking GMP or MPFR or similar libraries. It can be done, and the speed benefits would be nice, but the extra time to write the code nullifies any benefits. I'd only resort to c++ with MPFR if I needed to solve a particularly slow/complex problem. For example, some day I might port my accelerated slog solution.

Anyway, if you can use PARI/gp, the following code might be useful to you. (Or to anyone else, like Sheldon?)


.gp   jsloglib_v0.5.gp (Size: 5.47 KB / Downloads: 122)

To test the fixed point stuff:

.gp   jslog_fptest.gp (Size: 305 bytes / Downloads: 130)

To test the Schroeder functions (this is more for Sheldon's benefit):

.gp   jslog_test.gp (Size: 933 bytes / Downloads: 128)

Here's an example of the output for the jslog_fptest.gp attachment:
   
~ Jay Daniel Fox
Reply
#4
I have pari/GP.
Thanks for reply.
I will try it in some spare time and try to convert to c++.
If I manage to do that, I will post C++ versions.

Thanks!
Fuji GSW690III
Nikon D3, Nikkors 14-24/2.8, 24/1.4, 35/2, 50/1.4, 85/1.4, 135/2, 80-200/2.8
Reply
#5
(09/18/2014, 04:31 AM)MorgothV8 Wrote: I have pari/GP.
Thanks for reply.
I will try it in some spare time and try to convert to c++.
If I manage to do that, I will post C++ versions.

Thanks!

Seems to woark really good, I've edited a little and played with different bases, even close to 0 or 1 or even negative ones.

THANKS a lot.
If only I find time, I will convert it to C++ and if You agree repost.
C++ have really limited built-in functions, Mathematica makes it so easy.

Anyway - I want to understand how it is working, not just copy it.




Attached Files
.gp   jslog_fptest.gp (Size: 1.34 KB / Downloads: 105)
Fuji GSW690III
Nikon D3, Nikkors 14-24/2.8, 24/1.4, 35/2, 50/1.4, 85/1.4, 135/2, 80-200/2.8
Reply
#6
(09/17/2014, 11:10 PM)jaydfox Wrote:
(09/17/2014, 04:17 PM)MorgothV8 Wrote: No one can help how to find all fixed points of log(b), for all complex bases b.

I know exp() have infinite numbe rof fixed points, but two of them are "main" and are the same as log(), I want to find them using some algorithm and elementary functions. I know there is no closed for for them, but I want to find algorithm that approximates them for all bases.

Any ideas?
My program ffrom 1st post inds fixed point - but only one, how to find another one?

(...) if you can use PARI/gp, the following code might be useful to you. (Or to anyone else, like Sheldon?)

So the code I posted previously will generate a Taylor series. I don't know for sure if the Taylor series is entire, or if it has a finite radius of convergence. So far I've managed to calculate the first 1024 terms (took several hours), and the root test sure looks like there's a radius of about 2.5.

The root test climbs up to about 2.77 by the 44th term, then begins descending (in a very cool pattern which I'll explore later).

The 77th term is the last term with a root test over 2.70
The 125th term is the last term with a root test over 2.65
The 233rd term is the last term with a root test over 2.60

I don't have the larger data set (1024th degree), but eyeballing a graph I made, it looks like the root test is a little under 2.49 by that point in the Taylor series (terms 1 mod 4, e.g., 1017th term, 1021st term, etc.).

I don't know how long the root test continues to get smaller, before (or if?) it gets bigger. It could be tens of thousands of terms, which is beyond my ability to calculate. Using my algorithm, we'd need a super-computer. So either we need a better algorithm, or we need some other way of approaching the problem.

By the way, it's a shame that the series is effectively limited to a radius of about 2.5. I found that I can find non-primary fixed points using the same formula. For example, the primary fixed points of base e are 0.3181315 +/- 1.3372357*I. Well, when I plugged sqrt((log(log(e)) + 2*Pi*I)+1) into the formula with 320 terms, I got a value of 2.064031 - 7.59095*I, which isn't too far from a fixed point, 2.06227773 - 7.58863118*I. I assume with more terms it will be even closer to that fixed point. (My 1024-term solution is temporarily unavailable, so I can't test with it at the moment.)

So my working hypothesis is that if we had the full Taylor series, we could calculate primary and at least some (if not all) of the non-primary fixed points for any base (other than 1 or 0). Since the full Taylor series is out of reach, we might look at the idea of analytically extending the function and then developing the Taylor series somewhere that has a larger effective radius (as indicated by the root test).
~ Jay Daniel Fox
Reply
#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
#8
ADD: code seems to have a convergence problem on some bases for the principal branch. Will need to look into this.
Reply
#9
That's strange ... can't reproduce the apparent problem. Anyways, this code:

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 = (0.639 + 0.152*I)*log(z + oneovere) + (0.390 - 0.343*I);
     );
   /* Large asymptotic for large values of z */
   if(abs(z + oneovere) >= 10,
      r = log(z) - log(log(z)) + (log(log(z))/log(z));
     );

   /* ensure purely real result for positive and not-too-negative real inputs */
   if((imag(z) == 0) && (real(z) >= -oneovere), r = real(r));

   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 + 1) - ((r + 2)*(r*expr - z))/(2*r + 2));
        );

   /* 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);

fix(n, b) = Wbr(n, -log(b))/(-log(b));

has an improved initial guess for the principal branch. Use "fix(n, b)" to find the nth fixed point of base b.

This method is highly efficient, being able to compute the fixed points to thousands of digits in seconds on my 2.4GHz Core 2 machine. E.g. fix(4, 2) gives

Code:
5.2738486587969421374330315218571433423534146415852330401066294382423556546501676001880571055973904853663679469274813764620000693663488878468878992575133848048931917567327477719613549216850000740670739245408605739020265724724505221922149165357780497139655091703495699587502138683716726663081682755519618348414500732021219035238961163860776131773637452101798940688378851967487213699777902485555983113779509422469571944880078361570690436015178078493425261482043433724684114279337464744511646424440703191308122783279 -
38.327787228780575754183706624252395882735744965637005425264867291957073919505982312610787328357650749246081503608230774425033102729695628681018727886229348965817703601243873934949849312466874132625859591878503716732102631968981649986915683789960123279085687250569540219407639546682348231751332415437790604182399983123646641624765667027493935467428011331004110282142632651966995879343798360074353334590171964725501188490595138754927696831217983525656552968160304620509175861943833510599624314954439880911317940197*I

512 digits of precision almost instantly.
Reply
#10
Addendum: You can replace Wapprox with this:

Code:
Wapprox(k, z) = {
   local(brl = brlog(k, z));
   local(lbrl = log(brl));

   return(brl - lbrl + sum(kk=0,5,sum(m=1,5,c(kk, m)*lbrl^m * brl^(-kk-m))));
}

to get even more speed (saves some repeated log calculations).
Reply


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



Users browsing this thread: 1 Guest(s)