iconEuler Examples

Julia and Mandelbrot Sets

by R. Grothmann

In this notebook, I want to study simple Julia sets and the Mandelbrot
set.

First we start by defining an iteration procedure iter(z,c), which
iterates

Mandelbrot

for some numbers in the complex plane until the maximal value of the
numbers z gets very large. The function returns the values of z_n and
the number of iterations.
>function iter (z,c,n=100) ...
h=z;
loop 1 to n;
  h=h^2-c;
  if totalmax(abs(h))>1e20 then 
    m=#; break; 
  endif;
  end;
return {h,m};
endfunction
Now we apply the iteration to numbers in [-2,2]x[-2,2] in the
complex plane.
>c=0.8; ..
>x=-2:0.05:2; y=x'; z=x+I*y; ...
>{w,n}=iter(z,c);
We scale the numbers, determining the logarithmic growth of the
numbers.
>wr=max(0,log(abs(w)))/2^n;
We get an approximation of the logarithmic potential of the Julia set
of z^2-c. This set is the set of all points, where the iteration
remains bounded.
>plot3d(re(z),im(z),wr,hue=1,zoom=4); insimg();

Mandelbrot

Here is another plot of the same set (as accurately as our grid
was).
>density(wr>0.0001); setplot(-2,2,-2,2); xplot(); insimg;

Mandelbrot

The Mandelbrot Set

Now we show the Mandelbrot set. This is the set of all c, such that
the iterate z=z^2-c is bounded, when starting with z=0. I.e., the
orbit of 0 is bounded.
>x=-2:0.05:3; y=(-2.5:0.01:2.5)'; z=x+I*y;
We can use the same iteration, since c ist just replaced by z.
>{w,n}=iter(zeros(size(z)),z); ...
>wr=2*max(0,log(abs(w)))/2^n; ...
>plot3d(re(z),im(z),wr,hue=1,light=[1,1,0],zoom=4,height=45°,>insimg);

Mandelbrot

The density is plot of this coarse approximation.
>density(wr>0.0001); setplot(-2,2,-2,2); xplot(); insimg;

Mandelbrot

The following function iterates a 100 times. The z larger than 100 in
absolute value are set to 100 in each iteration. This is an effective
way to use the Euler matrix language.
>function mandelbrot (a,b,c,d,n=500,m=100) ...
z=linspace(a,b,n)+1i*linspace(c,d,n)'; h=zeros(size(z));
loop 1 to m
  h=h^2-z;
  h=(abs(h)>10)*10+(abs(h)<=10)*h;
end
density(log(max(abs(h),1))>0.001); 
setplot(a,b,c,d); 
xplot();
return h;
endfunction
The following takes about 10 seconds on my computer. We iterate on the
default 500x500 grid.
>mandelbrot(-0.1,1,-0.1,1); insimg;

Mandelbrot

Using Compiled Code

In the following library function mandel(z), we don't iterate all
values of the matrix 100 times, but each value only until it exceeds
10 in absolute value. The return value is the number of iterations it
took to exceed the maximal value.

You cannot compile this on the installation directory. You need to
copy the notebook and mandelbrot.c to a directory with write access,
if you want to try the compilation.
>// closedll("mandelbrot"); tccompile mandelbrot
Load the library function.
>dll("mandelbrot","mandel",2);
>function mandelbrot (a,b,c,d,n=500,m=100) ...
z=linspace(a,b,n)+1i*linspace(c,d,n)'; 
res=mandel(z,m);
density(m-res/m); 
setplot(a,b,c,d); 
xplot();
endfunction
>mandelbrot(-0.1,1,-0.1,1); insimg;

Mandelbrot

A closer look.
>mandelbrot(0.4,0.8,0.4,0.8); insimg;

Mandelbrot

Here is a vary intimate detail with more iterations for more sharpness
in details.
>mandelbrot(0.746,0.75,0.072,0.076,m=1000); insimg;

Mandelbrot

Zooming in once more.
>mandelbrot(0.747,0.749,0.074,0.076,m=1000); insimg;

Mandelbrot

Here is a more colorful version of the plot, using the plotrgb
function, and some color manipulation.
>function cmandelbrot (a,b,c,d,n=500,m=100) ...
z=linspace(a,b,n)+1i*linspace(c,d,n)'; 
res=mandel(z,m);
res=1-(res-totalmin(res))/(totalmax(res)-totalmin(res));
res=sqrt(res);
window(0,0,1024,1024);
plotrgb(rgb(res*0.5,res,(res+res*res)/2)); 
shrinkwindow;
endfunction
>cmandelbrot(0.747,0.749,0.074,0.076,m=1000); insimg;

Mandelbrot

Mandelbrot Curves

The following curve builds the Mandelbrot set.
>t=linspace(0,2*pi,500); z=(exp(2*I*t)-2*exp(I*t))/4;
>plot2d(re(z),im(z),r=2); insimg;

Mandelbrot

It forms the main chunk of bounded iterations. The other parts are
parts, which happen to fall into this set after some iterations.
>mandelbrot(-1,1,-1,1); ...
>  plot2d(re(z),im(z),r=1,add=1,thickness=2,color=yellow);  ...
>  insimg;

Mandelbrot

There is a second part, which is completely coverted with bounded
iterations.
>mandelbrot(0,2,-1,1); ...
>  t=linspace(0,2*pi,500); z=1+exp(t*I)/4; ...
>  plot2d(re(z),im(z),a=0,b=2,c=-1,d=1,>add,color=yellow,thickness=2);  ...
>  insimg;

Mandelbrot

Examples Homepage