The Gamma function is an extension of the factorial (n!) function, with its argument shifted down by 1, to real and complex numbers. And there is a fractal associated to it. I just wanted to learn how to make it because I want to learn how to plot the Mandelbrot set, and the Gamma function fractal seems easier to catch the idea of recursion. So here is how to make it!
First of all, I have been following the great explanation at Mathistopheles website, but I needed to study by myself some details to make it work properly.
Basically, when the gamma function is calculated for a complex number z repeatedly (nested, and shifting down 1 the complex number according to the definition), this is: Gamma(Gamma(...(Gamma(z-1)...-1)-1), the final value result of the nested calculations either tends to 1 or grows without limit. So it is possible to map each z (=a+bi) value of the complex plane (in Cartesian coordinates) with a color code showing if the value tends to 1 (this the same as saying that it tends to the complex number 1+0i) or not in a specific numbers of calculations of the nested Gamma function.
The condition to tend to 1 is as follows: let us say we want to know if the nested calculation of the Gamma function applied to z = (a+bi) tends to 1. In a first step, calculate Gamma(z-1) = c+di. If the distance of the real part a to 1 is greater or equal to the distance of c to 1 and the absolute value of the imaginary part b is greater than the absolute value of the imaginary part d (abs(b)>=abs(d)) then it is possible to say that z tends to 1+0i when the Gamma function is applied over z (nested as it was explained above).
It is important to understand that it can happen that the first calculation for a given complex number z1 of Gamma(z1-1) provides a complex number z2 which is not closer to 1 than the original z1. And that does not mean that z1 does not tend to 1 when the nested Gamma function is applied again to z2. It could happen that when Gamma(z2-1) = z3 is calculated, z3 is closer to 1 than z2, and in the moment that this happens then the result is that from that point z1 tends to 1 when the nested calculations are applied. So sometimes we just need one calculation of Gamma for the initial complex number z1 to know if it tends to 1 and sometimes we need two nested Gamma calculations, sometimes 3, 4, 5, ... etc sometimes hundreds or thousands. But the calculation of the computer is very slow when the Gamma function is applied hundreds of times for each complex number, so usually a limit of for instance 7 nested Gamma calculations is applied to decide if a complex number tends to 1 or not. It is not totally accurate, but good enough to start to see the fractal patterns of the Gamma function. The more nested calculations the more accurate the patterns will be, and of course, the more complex numbers are calculated, the more we can zoom into the fractal pattern.
And here is the result:
It is the typical image of the Gamma fractal function in the complex plane interval [-5,0]+[-1,1]i. I tried to find a good combination of colors. The complex numbers that tend to 1 in just one Gamma calculation are in black color. Then Dark blue two nested calculations, then Pink... finally up to Yellow (7 nested loops). The complex numbers that do not tend to 1 in less than 8 nested calculations are in the Blue color (if the reader can rememeber the good old Fractint, it was the popular color used for that purpose!).
Just to finish here is another version in polar coordinates:
And here is the Python code, modifying the value of the "method" variable it is possible to render Cartesian, Polar or Spherical coordinates. Use it freely:
def gamma_fractal():
import matplotlib.pyplot as plt
from scipy.special import gamma
from numpy import real, imag, around, arange, float, cos, sin, pi
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
#0 Cartesian, 1 Polar, 3 Spherical
method = 0
if method==3:
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
arraydim = 1000000
xrange = [arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float)]
yrange = [arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float)]
zrange = [arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float),arange(arraydim,dtype=float)]
hexcolor=["#000000","#131667","#F028B3","#42677C","#3EA670","#6DD02D","#F9B90F","#0000FF"]
pos=[0,0,0,0,0,0,0,0]
ival = 0+1j
if method!=3:
firstx = -5
lastx = 0
firsty = -1
lasty = 1
incx = 0.0005
incy = 0.0005
else:
firstx = -pi
lastx = pi
firsty = -pi
lasty = pi
incx = 0.002
incy = 0.002
posx = firstx
while posx < lastx:
posy = firsty
while posy < lasty:
counter = 0
mycomplex = gamma((posx+1)+(posy*ival))
nextcomplex = gamma(mycomplex+1)
while counter < 7:
if (abs(mycomplex.imag) > abs(nextcomplex.imag)) and (abs(1-mycomplex.real) > abs(1-nextcomplex.real)):
break
mycomplex = nextcomplex
nextcomplex = gamma(mycomplex+1)
counter = counter + 1
if method == 0:
xrange[counter][pos[counter]] = posx
yrange[counter][pos[counter]] = posy
elif method == 1:
xrange[counter][pos[counter]] = posx*cos(posy)
yrange[counter][pos[counter]] = posx*sin(posy)
else:
xrange[counter][pos[counter]] = cos(posy)*cos(posx)
yrange[counter][pos[counter]] = sin(posy)*cos(posx)
zrange[counter][pos[counter]] = sin(posx)
pos[counter] = pos[counter] + 1
if pos[counter]==arraydim:
newxrange = arange(pos[counter],dtype=float)
newyrange = arange(pos[counter],dtype=float)
newzrange = arange(pos[counter],dtype=float)
for i in range(0,pos[counter]):
newxrange[i]=xrange[counter][i]
newyrange[i]=yrange[counter][i]
if method == 3:
newzrange[i]=zrange[counter][i]
if method != 3:
plt.plot(newxrange,newyrange,",",color=hexcolor[counter])
else:
ax.plot(newxrange, newyrange, newzrange,",",color=hexcolor[counter])
pos[counter]=0
posy = posy + incy
posx = posx + incx
for c in range (0,8):
newxrange = arange(pos[c],dtype=float)
newyrange = arange(pos[c],dtype=float)
newzrange = arange(pos[c],dtype=float)
for i in range(0,pos[c]):
newxrange[i]=xrange[c][i]
newyrange[i]=yrange[c][i]
if method == 3:
newzrange[i]=zrange[c][i]
if method != 3:
plt.plot(newxrange,newyrange,",",color=hexcolor[c])
else:
ax.plot(newxrange, newyrange, newzrange,",",color=hexcolor[c])
pos[c]=0
if method!=3:
ax = plt.gca()
ax.set_axis_bgcolor((0, 0, 0))
figure = plt.gcf()
figure.set_size_inches(30, 20)
if method == 0:
plt.axis([firstx,lastx,firsty,lasty])
plt.savefig("gamma_"+str(method)+".png")
plt.show()
else:
#ax.legend()
#ax.dist=8
for ii in range(0,361):
ax.view_init(elev=-178., azim=ii)
mpl.pyplot.savefig("movie"+str(ii)+".png")
#ax.view_init(elev=28., azim=ii)
#plt.show()
gamma_fractal()
Enjoy the fractals!
No comments:
Post a Comment