Forms of life, forms of mind

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. Below 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.

There are a few aspects of this that go beyond their aesthetic appeal. First, this is actually a map of the Halting Problem, in a way – one of the colors indicates places where the root-finding method has not converged in many iterations (and is unknown when or if it will ever do so). What you are seeing is the spatial representation of the fact that between any two points that lead the algorithm to halt might exist ones for which it does not. Second, these patterns are an amazing example of complex form that is determined neither by any aspects of physics nor by any kind of history (e.g., an evolutionary explanation as used in biology). These images illustrate mathematical objects whose specific structure is not explainable by the usual methods of physics or the life sciences, giving us a glimpse into a non-physical source of order: the world of mathematics. See this talk for the implications of this for understanding biology.

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