In [13]:
# Math libraries
import numpy as np
from math import sqrt
# Matplotlib plotting libraries
import matplotlib.pyplot as plt
%matplotlib inline
# Bokeh plotting libraries
import bokeh.plotting as blt
blt.output_notebook()
In [2]:
class JuliaSet(object):
def __init__(self,c,n=100):
self.c = c
self.n = n
self._d = 0.001
self._complexplane = np.array([])
self.set = np.array([])
def juliamap(self,z):
return z**2 + self.c
def iterate(self,z):
m=0
while True:
z=self.juliamap(z)
m=m+1
if abs(z)>2:
return m
if m>=self.n:
return 0
def set_spacing(self,d,l_inf=-2,l_sup=2):
self._d=d
self.l_inf = l_inf
self.l_sup = l_sup
t=np.arange(self.l_inf,self.l_sup,self._d)
x,y=np.meshgrid(t,t)
self._complexplane=x+y*1j
def generate(self):
iterate_vector=np.vectorize(self.iterate)
self.set=iterate_vector(self._complexplane)
return self.set
In this last version of JuliaSet, the implementation of 'set_spacing' has changed. Instead of creating 'complexplane' element by element, now 'set_spacing' prepares the matrix x with repetitions of the vector composed by the elements of the "arange" and its transposed matrix y. After this, the method sums x and 1j times y. It produces a more efficient algorithm.
Another change in that method is the creation of the variables l_inf and l_sup, that allows the alteration of the inferior and superior limits for the range. Because of this, some parts of the class JuliaSetPlot were modified, replacing the fixed limits -2 and 2 for the new variables.
The last modification was made in the method 'generate' with the vectorization of the function 'iterate'. This change is not related to the efficience of the algorithm, but it helps to get a cleaner code.
In [3]:
class JuliaSetPlot(JuliaSet):
"""Extend JuliaSet to add plotting functionality"""
def __init__(self, *args, **kwargs):
# Invoke constructor for JuliaSet first, unaltered
super(JuliaSetPlot, self).__init__(*args, **kwargs)
# Add one more attribute: a rendered image array
self.img = np.array([])
def get_dim(self):
# get what should be an attribute
return int(sqrt(self.img.size)) #int(4.0 / self._d)
def render(self):
if not self.set: self.generate()
# Convert inefficient list to efficient numpy array
self.img = np.array(self.set)
dim = self.get_dim()
# Reshape array into a 2d complex plane
self.img = np.reshape(self.img, (dim,dim)).T
def show(self):
if not self.img.size: self.render()
# Specify complex plane axes efficiently
xy = np.linspace(self.l_inf,self.l_sup,self.get_dim())
# Use matplotlib to plot image as an efficient mesh
plt.figure(1, figsize=(12,9))
plt.pcolormesh(xy,xy,self.img, cmap=plt.cm.hot)
plt.colorbar()
plt.show()
def interact(self):
from matplotlib.colors import rgb2hex
if not self.img.size: self.render()
# Mimic matplotlib "hot" color palette
colormap = plt.cm.get_cmap("hot")
bokehpalette = [rgb2hex(m) for m in colormap(np.arange(colormap.N))]
# Use bokeh to plot an interactive image
f = blt.figure(x_range=[self.l_inf,self.l_sup], y_range=[self.l_inf,self.l_sup], plot_width=600, plot_height=600)
f.image(image=[j.img], x=[self.l_inf], y=[self.l_inf], dw=[self.l_sup-self.l_inf], dh=[self.l_sup-self.l_inf], palette=bokehpalette)
blt.show(f)
In [20]:
j = JuliaSetPlot(1)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [9]:
j = JuliaSetPlot(1j)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [17]:
j = JuliaSetPlot(-1)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [16]:
j = JuliaSetPlot(-1j)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
(Why are the plots for c=1 and c=-1 so differents from each other?)
It could be interesting to check if there is some kind of symmetry in the axes rotated by $\pi/4$:
In [13]:
j = JuliaSetPlot(-1/sqrt(2) - (1/sqrt(2))*1j)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [12]:
j = JuliaSetPlot(1/sqrt(2) + (1/sqrt(2))*1j)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [18]:
j = JuliaSetPlot(1/sqrt(2) - (1/sqrt(2))*1j)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [19]:
j = JuliaSetPlot(-1/sqrt(2) + (1/sqrt(2))*1j)
%time j.set_spacing(0.006)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
The following two plots were used to analyse the stability of the solution for imaginary entrances:
In [4]:
j = JuliaSetPlot(9j)
%time
j.set_spacing(0.006,-5,5)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [5]:
j = JuliaSetPlot(0.1+9j)
%time
j.set_spacing(0.006,-5,5)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
It is very interesting to observe the behavior of the following plots too:
In [15]:
j = JuliaSetPlot(1+1j)
%time
j.set_spacing(0.006,-10,10)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [5]:
j = JuliaSetPlot(2+2j)
%time
j.set_spacing(0.006,-10,10)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [8]:
j = JuliaSetPlot(4+4j)
%time
j.set_spacing(0.006,-10,10)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [7]:
j = JuliaSetPlot(8+8j)
%time
j.set_spacing(0.006,-10,10)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [8]:
j = JuliaSetPlot(32+32j)
%time
j.set_spacing(0.006,-10,10)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [14]:
j = JuliaSetPlot(64+64j)
%time
j.set_spacing(0.006,-10,10)
print len(j._complexplane)
print sqrt(len(j._complexplane))
%time j.show()
In [ ]: