# Forms of life, forms of mind

Dr. Michael Levin

## Halley’s Method fractal art

on

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);
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``````
Next Post

## 2 responses to “Halley’s Method fractal art”

1. Geraldo Medeiros Junior

It’s impossible to read this and not correlate it with the bioelectrical network of a biological structure. I believe that, although fractal theory is often associated with the structure and dynamics of complex systems, it can also be applied to other aspects of non-neural bioelectrical communication in organisms. Some parallels:
Electrical conduction system of the heart: The electrical conduction system of the heart is responsible for generating electrical signals that coordinate the contractions of the cardiac muscle. The way these signals propagate throughout the heart exhibits self-similarity at various scales. For example, electrical signals are generated in the sinoatrial node (SA) and spread through branching, including the atrioventricular node (AV) and conducting fibers, forming a complex fractal-like pattern.
Electrical signals in muscle cells and other cells: In various muscle cells and other cell types, there are bioelectrical processes involving the release of ions and the generation of electrical potentials. The way these signals propagate through cells and tissues can exhibit fractal properties, especially in situations where branching and diffusion processes are relevant.
Propagation of electrical potentials in plants: In plants, bioelectrical communication occurs through the propagation of electrical potentials along conducting cells. Conduction structures in plants, such as xylem and phloem vessels, can exhibit complex branching patterns that resemble fractals.
Cellular communication networks: Communication between cells in an organism can involve the transmission of electrical signals through gap junctions and other communication structures. These cellular communication networks can have a highly branched and self-similar structure, allowing for efficient exchange of information between cells in different parts of the body.
Thus, fractal theory can indeed be applied to non-neural bioelectrical communication in organisms, where patterns of branching and propagation of electrical signals exhibit self-similarity at various scales. These parallels help highlight the complexity and efficiency of bioelectrical systems in the communication and coordination of functions in organisms.

2. Hi Mike!

I’m a fan of your work, and I happen to work on MATLAB at MathWorks. Thanks for sharing your code. I got inspired and started playing around with these calculations myself. I’ll share my code below, in case you’re interested. It only works with polynomials for now.

Ned.

function [rootMap,iterMap,xs,ys] = root_basins(p,xlims,ylims,n,method)
%ROOT_BASINS Determine the root basins for a polynomial using iterative methods.
%
% [ROOTMAP, ITERMAP, XS, YS] = ROOT_BASINS(P) computes the root basins
% for a given polynomial P in the complex plane using Newton’s method.
% P is a vector that contains the polynomial coefficients in decreasing
% order of degree. The function returns the following:
% – ROOTMAP: An NxN matrix where each element indicates to which root of
% the polynomial a particular point in the complex plane converges.
% – ITERMAP: An NxN matrix that indicates the number of iterations taken
% for each point in the complex plane to converge to a root.
% – XS and YS are vectors that define the grid in the complex plane.
%
% ROOT_BASINS(P, XLIMS, YLIMS) specifies the limits for the real
% (XLIMS) and imaginary (YLIMS) parts of the complex plane. By default,
% XLIMS and YLIMS are set to [-1 1].
%
% ROOT_BASINS(P, XLIMS, YLIMS, N) specifies the resolution of the
% grid in the complex plane. N is the number of points in both the x
% and y directions. By default, N is set to 200.
%
% ROOT_BASINS(P, XLIMS, YLIMS, N, METHOD) specifies the iterative
% method to use. METHOD can be either ‘newton’ (default) or ‘halley’.
%
% Example:
% p = [1 0 0 -1];
% [rootMap, iterMap, xs, ys] = root_basins(p);
% imagesc(xs,ys,rootMap)
%

arguments
p (1,:) double
xlims (1,2) double = [-1 1]
ylims (1,2) double = [-1 1]
n (1,1) double = 200
method {mustBeMember(method,[“newton”,”halley”])} = “newton”
end

rt = roots(p);

% Take the necessary derivatives
dp = polyder(p);
ddp = polyder(dp);

% Make some helpful anonymous functions
f = @(q) polyval(p,q);
df = @(q) polyval(dp,q);
ddf = @(q) polyval(ddp,q);

xs = linspace(xlims(1),xlims(2),n);
ys = linspace(ylims(1),ylims(2),n);

[x,y] = meshgrid(xs,ys);
z = x + 1i*y;

tol = 1e-6;
rootMap = zeros(size(z)); % Track roots
iterMap = zeros(size(z)); % Track iteration count
freezeMap = false(size(z)); % Track where you are done calculating

for k = 1:100

% Only calculate on unfrozen regions
zz = z(~freezeMap);

if method == “newton”
zz = zz – f(zz)./df(zz);
elseif method == “halley”
zz = zz – 2*f(zz).*df(zz)./(2*(df(zz).^2 – f(zz).*ddf(zz)));
else
error(“No such method”)
end

% Put the newly calculated values back into z
z(~freezeMap) = zz;

% Loop over the roots to see where we’re close enough to be done
for m = 1:length(rt)

zerosFoundMap = abs(rt(m)-z) < tol;
% Only update when you find zeros in unfrozen regions
zerosFoundMap = zerosFoundMap & ~freezeMap;

if any(zerosFoundMap(:))
rootMap(zerosFoundMap) = m;
iterMap(zerosFoundMap) = k;
% Add this region to the frozen region so it will no longer
% be updated or calculated on
freezeMap = freezeMap | zerosFoundMap;
end

end

if all(freezeMap(:))
% Exit early if all the roots are known
break
end

end

end