Drawing Mc by Jungreis Algorithm




Jungreis algorithm gives inverse of Boettcher function.



From Douaddy-Hubbard theory it is known that that Mc is connected.
It is proved by constructing a analytic homeomorphism (bijection) F of  Dc onto Mc

c = Psi_M (w)
for c: c belongs to Mc and w: w belongs to Dc
So one can transform complement (exterior) of closed unit disk to complement (exterior) of Mandelbrot set

Douaddy-Hubbard" construction does not lend itself to computation. "
Jungreis "gives alternative construction which does."

It allows to compute points of external rays R(angle) and equipotential lines in w-plane using definitions:
R(angle) = radius * e i * angle
where 1 < radius <= infinity
E(radius) = radius * e i * angle
where 0 < angle <= 1 [ turns]

and after that transform it to c-plane using Psi_M and draw.

Isn't it beautifull ?

(Here I will skip theory and will give only final result)

Because function Psi_M is analytic one can compute it using power series

Psi_M(w)= w + Sum ( bm * w -m)
where 0 <=m <=infinity

It is sufficiant to compute some ( finite and not very large) number of terms.
Jungreis uses 4095 terms of power series to achieve accurate ploting result.
First one should compute coefficients of Psi_M.


list of coefficients of the mapping exterior of unit disk to the exterior of the Mandelbrot set


Definition of Dc

Set D is closed unit disk
Set Dc is a complement of D
Dc= {w: 1 < |w| <= infinity }


Hera are the circles with center=0 and radius: 1, 1.0001, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.4, 1.5
and ninty equaly spaced radial lines ( external rays of angles = p/90 : 1<=p<=90 measured in turns ).
These circles are equipotential lines for ....
It is the initial image for transformation ( see image in Jungreis paper)
and here is a program that draws this image ( in delphi for win32).
Definition of Mc

c,z,w are complex numbers
k is integer number
let Pc is a complex quadratic polynomial function
Pc(z) = z2+c
Pc(k) is the k'th iteration of Pc
( iteration gives set of points {z0, z1, z2, .. ,zn}
named orbit)
Mandelbrot set M = { c : Pc(k) is bounded }
Set Mc is complement of M

On this image Mc is white and M is black.
It is made using boolean escape time ( Mandel ) algorithm.
b(m) = coefficients of the mapping exterior of unit disk to the exterior of the Mandelbrot set. ( Jungreis algorithm )
( translation of maple code  by G. A. Edgar

------------------( function version)-------------------------------
u(n,k):=block([j],
 if 2^n-1=k then 1
            else if 2^n-1>k then sum(u(n-1,j)*u(n-1,k-j),j,0,k)
                            else if 2^(n+1)-1>k then 0
                                                else (u(n+1,k)-sum(u(n,j)*u(n,k-j),j,1,k-1))/2
);

a(j):=u(0,j+1);

w(n,m):=block
([j],
if n=0 then 0
else a(m-1)+w(n-1,m)+sum(a(j)*w(n-1,m-j-1),j,0,m-2)
);

b(m):= if m=0 then -1/2
else -w(m,m+1)/m;



---------------( array version)--------------------------------------

u[n,k]:=block([j],
if 2^n-1=k then 1
else if 2^n-1>k then sum(u[n-1,j]*u[n-1,k-j],j,0,k)
else if 2^(n+1)-1>k then 0
else (u[n+1,k]-sum(u[n,j]*u[n,k-j],j,1,k -1))/2);

a[j]:=u[0,j+1];

w[n,m]:=block
([j],
if n=0 then 0
else a[m-1]+w[n-1,m]+sum(a[j]*w[n-1,m-j-1],j,0,m-2)
);

b[m]:= if m=0 then -1/2
else -w[m,m+1]/m;

thx to G. A. Edgar and Richard J. Fateman for help (:-))

Fast version with numerical backward recursion ( transaltion of maple code by G. A. Edgar

betaF(n,m):=block
(
[nnn:2^(n+1)-1],
if m=0
    then 1.0
    else if ((n>0) and (m < nnn)) then 0.0
                                else (betaF(n+1,m)- sum(betaF(n,k)*betaF(n,m-k),k,nnn,m
 -nnn)-betaF(0,m-nnn))/2.0
);
b(m):=betaF(0,m+1);

The fastest version !!!!!!!!!!!!!!!!!!!:

Robert P. Munafo have noticed that : (above) " maxima code does not gain the efficiency of
remembering previously computed values (as the "option remember" directive does in Maple)"

I have asked Richard J. Fateman for help, and he said:
"replace u(n,k) by u[n,k] everywhere, and it will remember values, because u is now
an array instead of a function." RJF

so here is latest version:

betaF[n,m]:=block
(
[nnn:2^(n+1)-1],
if m=0
then 1.0
else if ((n>0) and (m < nnn)) then 0.0
else (betaF[n+1,m]- sum(betaF[n,k]*betaF[n,m-k],k,nnn,m
-nnn)-betaF[0,m-nnn])/2.0
);

b(m):=betaF[0,m+1];


One can do a list with:
for m:1 thru 20 step 1 do display(b(m));
or
for i:0 thru 64 do ( print(sconcat("b[", i, "]=", b(i))) );
list of coefficients


Main Page
Feel free to e-mail me!
Author: Adam Majewski
adammaj1-at-o2-dot-pl
About

http://republika.pl/fraktal/