Fractal art can be made by testing each point (x,y) in a region of the complex numbers plane to see how soon a given function (in complex numbers) converges on its root, when using the number x+yi as the starting point of a root-finding method. The rate at which it settles on a root determines the color. Here are some fractals created using Halley’s method. Videos are created by having each frame be a fractal of a slightly different formula (change some parameter a tiny bit over a range for each frame, and then assemble the frames into a movie). Some videos look like changing shapes, while others look like a “fly-over” a static image. Code to generate these in Matlab is given below. Typical parameters are:
minx=-5, maxx=1, miny=-7, maxy=-2.5,
formua = z*z
resolution = 1024, frames = 500
You can zoom in infinitely, into ever smaller regions of the image and find novel structure. Note how organic the forms look, and, how this incredible complexity emerges from a tiny seed – the function (plus algorithm) have very small information content; the actual shapes come from where? Where-ever the truths of mathematics come from… Another question is: what would it take to compress such an image to the short formula that produces it (way more compact than what a typical compression algorithm produces)? That is closely related to the inverse problem of development (how to figure out which gene(s) to change to edit a specific anatomical shape in just the way you want).
Here are some videos:
Prior papers on this method of making fractals:
My Matlab code: (it asks you to put in the function’s derivative manually; if you have Matlab’s symbolic computation toolbox, you can easily change the code to have it do that automatically)
function Newtonsmethod
global cancelled;
lambda = 1;
cancelled = 0;
r=0;
e = 2.71828182845904523;
epsilon = 0.00001;
% make the default colormap.
clrmp = [0 0 0; 0 0 1; 0 1 0; 1 0 0; 0 1 1; 1 0 1; 1 1 0; 1 1 1];
fname = input('File to save movie in: ','s');
steps = 0;
r_begin = 0;
r_end = 0;
% (sum(r_end)==0) || sum(r_begin)==0 ||
while (r_end<=r_begin)
r_begin = input('Enter begin value for real axis: ');
r_end = input ('Enter end value for real axis: ');
end
i_begin = 0;
i_end = 0;
% (sum(i_end)==0) || sum(i_begin)==0 ||
while (i_end<=i_begin)
i_begin = input('Enter begin value for imaginary axis: ');
i_end = input ('Enter end value for imaginary axis: ');
end
while (steps<32)
steps = input('Enter resolution in pixels: ');
end
i_step = (i_end-i_begin)/steps;
r_step = (r_end-r_begin)/steps;
formula=input('Enter formula involving the complex variable z and movie parameter t: ','s');
derivative=diff(formula,'z');
q1 = sprintf('First derivative is %s',char(derivative));
disp(q1);
% sanity check for 1st derivative
% t=1.0;
% z = 10.0;
% zy1 = eval(formula);
% z = 10.0+epsilon;
% zy2 = eval(formula);
% slope = (zy2-zy1)/epsilon;
% z = 10.0+epsilon/2.0;
% zz = eval(derivative);
% if (abs(zz-slope)>epsilon)
% disp('That does not seem like the right first derivative!');
% end
derivative2 = diff(derivative,'z');
q2 = sprintf('Second derivative is %s',char(derivative2));
disp (q2);
% sanity check for 2nd derivative
% t=1.0;
% z = 10.0;
% zy1 = eval(derivative);
% z = 10.0+epsilon;
% zy2 = eval(derivative);
% slope = (zy2-zy1)/epsilon;
% z = 10.0+epsilon/2.0;
% zz = eval(derivative2);
% if (abs(zz-slope)>epsilon)
% disp('That does not seem like the right second derivative!');
% end
movie_begin = 0;
movie_end = 0;
movie_step = input('How many frames for movie? ');
if (sum(movie_step)==0 || movie_step<2)
if (length(fname)>0)
filename = sprintf('%s.jpg',fname);
message = sprintf('Will make single frame into file %s.', filename);
disp(message);
else
disp('Will make single frame. ');
end
movie_begin_=0;
movie_end=1;
movie_step=1;
message = sprintf('Estimated completion time = %g hours.', ...
0.00481004 * (steps*steps*movie_step) /(60*60));
disp(message);
else
if (length(fname)>0)
filename = sprintf('%s.avi',fname);
message = sprintf('Will make movie into file %s, estimated size = %d Kbytes.', ...
filename, 0.984*movie_step*steps*steps/1000);
disp(message);
else
disp('Will make movie. ');
end
message = sprintf('Estimated completion time = %g hours.', ...
0.00481004 * (steps*steps*movie_step) /(60*60));
disp(message);
while (movie_end<=movie_begin)
movie_begin = input('Enter begin value for t: ');
movie_end = input ('Enter end value for t: ');
end
end
if (movie_step<0)
movie_step=1;
end
clrs = 0;
while (clrs<1)
clrs = input('Enter colormap to use: 1 = old style (8), N = smooth: ');
end
if (length(fname)>0)
datafile = sprintf('%s.txt',fname);
[fid message] = fopen(datafile,'w');
if (fid==(-1))
disp(message)
else
fprintf(fid,'minx=%g, maxx=%g, miny=%g, maxy=%g, \nt is from %g to %g, \nformua = %s.\n', ...
r_begin,r_end,i_begin,i_end,movie_begin,movie_end,formula);
fprintf(fid,'resolution = %d, frames = %d, method = Halley. \n',steps,movie_step);
end
else
fid = (-1);
end
if (clrs==1)
colormap(clrmp);
clrs = 8;
else
colormap(hsv(clrs));
end
cmap = colormap;
if (movie_step>1 && length(fname)>0)
aviobj=avifile(filename,'colormap',cmap,'compression','None','fps',7);
end
if (clrs==8)
klimit=15; % old-style algorithm
else
klimit=3*clrs; % for bigger colormap; probably takes longer!
end
begint = clock;
w_movie = waitbar(0,'Movie being made...','CreateCancelBtn',@cancel_callback);
% outer movie loop - each one is a frame of the movie
frames_done = 0;
B = zeros(steps+1,steps+1, 'uint8'); % static init for performance gain
axx = r_begin : r_step : r_end;
axy = i_begin: i_step: i_end;
stats = zeros(1,51,'uint16');
for t = movie_begin : ((movie_end - movie_begin)/movie_step) : movie_end
frames_done = frames_done + 1;
time1 = now;
% w_real = waitbar(0,'Real component...');
rr = 1;
for (r=r_begin : r_step : r_end)
if (cancelled==1 || not(ishandle(w_movie)))
break;
end
% plot settings
plothandle = imagesc(axx,axy,B);
ctitle = sprintf('Halley, formula = %s, t = %g ',formula,t);
title(ctitle);
drawnow;
% waitbar((r-r_begin)/(r_end-r_begin),w_real);
% w_imag = waitbar(0,'Imaginary component...');
ii=1;
for (i=i_begin: i_step: i_end)
% waitbar((i-i_begin)/(i_end-i_begin),w_imag);
% generate image
z = complex(r,i);
for (k=1:50)
z = z - eval(formula)/eval(derivative) / (1 - ...
(eval(formula)*eval(derivative2))/(2*power(eval(derivative),2.0)));
g = power(real(z),2)+power(imag(z),2);
if (k>1 && abs(g-gprev)<epsilon)
break;
end
gprev = g;
end
% stats(k) = stats(k) + 1;
if (abs(g-gprev)<0.00001 && mod(k,2)==0)
B(rr,ii) = int8(mod(k,clrs));
else
B(rr,ii) = 0;
end
ii = ii + 1;
end
% close(w_imag);
rr = rr + 1;
end
% delete(w_real);
if (cancelled==1 || not(ishandle(w_movie)))
break;
end
% plot settings
plothandle = imagesc(axx,axy,B);
ctitle = sprintf('Halley, formula = %s, t = %g ',formula,t);
title(ctitle);
drawnow;
if (movie_step>1 && length(fname)>0)
% add the frame to the movie file
% F = getframe(gca);
F = im2frame(B,cmap);
aviobj = addframe(aviobj,F);
elseif (length(fname)>0)
imwrite(B,colormap,filename,'jpg');
end
time2 = now;
timetook = time2 - time1;
time_remaining = timetook * (movie_step - frames_done);
estimate_done = datestr(now + time_remaining);
message = sprintf('Predicted endpoint = %s',estimate_done);
waitbar((t-movie_begin)/(movie_end-movie_begin),w_movie,message);
if (movie_step==1)
break;
end
end
endt = clock;
% close the files
if (movie_step>1 && length(fname)>0)
aviobj = close(aviobj);
end
if (fid==(-1))
disp('Warning: no data file created - remember the parameters!');
else
e = etime(endt,begint);
fprintf(fid,'Computation took %g minutes, or %g seconds/pixel. \n', ...
e/60.0,e/(steps*steps*movie_step));
% this should really convert to hours etc.
fclose(fid);
end
disp('All done. ');
if (ishandle(w_movie)) % otherwise it got cancelled and is already closed
delete(w_movie);
end
end
function n=reverse(z)
rr = imag(z);
ii = real(z);
n = complex(rr,ii);
end
function cancel_callback(h,varargin)
fig_handle=get(h,'parent');
delete(fig_handle);
cancelled = 1;
end
Leave a Reply