Filled Julia sets¶
MCS 275 Spring 2022 Lecture 27 - Emily Dumas
Orbit related functions¶
In [6]:
# this is the one whose filled Julia set I want to draw in this notebook
def g(z):
return z*z-1
# we discussed this one in the slides
def f(z):
return z*z
In [7]:
def orbit(func,a,n):
"""
Takes a function `func` of one paramater, and a starting
value `a`. Returns a list of length `n+1` consisting of
[a, f(a), f(f(a)), ... , f^n(a)]
"""
L = [a]
for _ in range(n):
# apply f to the last element computed
# and append the result
L.append( func(L[-1]) )
return L
In [8]:
# try orbit of 1.5 under z^2-1
# and 2 under z^2-1
orbit(f,2,5)
Out[8]:
In [9]:
orbit(f,1/2,5)
Out[9]:
In [10]:
orbit(g,2,10)
# a=2 is not in K_g
Out[10]:
In [11]:
orbit(g,1.5,30)
# a=1.5 seems to be in K_g
Out[11]:
In [12]:
# 0.5 + 0.5i is written as 0.5+0.5j in Python
orbit(g,0.5+0.5j,10)
# 0.5+0.5j is not in K_g
Out[12]:
In [13]:
def bounded_orbit(func,a,maxiter=10):
"""
Return True if the first `maxiter+1` elements of the orbit of `a` under `f`
are less than 2 in absolute value. Otherwise return false.
"""
z = a
for _ in range(maxiter):
if abs(z) >= 2:
return False
z = func(z)
# if we make it here, then f was applied maxiter times
return abs(z) < 2 # big z -> False (unbounded), small z -> True (bounded?)
In [14]:
bounded_orbit(f,2)
Out[14]:
In [15]:
bounded_orbit(f,1/2)
Out[15]:
In [16]:
bounded_orbit(g,1.5)
Out[16]:
In [17]:
bounded_orbit(g,0.5+0.5j)
Out[17]:
Making a grid¶
np.meshgrid takes a vector of x values and a vector of y values. It returns two matrices (one of x values, which change only in the rows, and one of y values, which change only in columns).
It's primarily used to make matrices that give coordinates of a rectangular grid in the plane, for example to make calculations using functions of two variables.
In [18]:
import numpy as np
np.meshgrid([1,2,3],[8,9,10])
Out[18]:
In [19]:
x = np.linspace(-2,2,500)
y = np.linspace(-2,2,500)
xx, yy = np.meshgrid(x,y)
In [20]:
print(xx.shape)
xx
Out[20]:
In [21]:
yy
Out[21]:
In [22]:
zz = xx + 1j*yy
In [23]:
zz # 50 x 50 matrix
Out[23]:
Classifying orbits¶
In [24]:
# Make bounded_order into a unary function by partial eval.
def has_bounded_orbit_under_g(a):
return bounded_orbit(g,a,maxiter=100)
In [25]:
has_bounded_orbit_under_g(1.5)
Out[25]:
In [26]:
# np.vectorize turns a func into a ufunc
orbit_classification = np.vectorize(has_bounded_orbit_under_g)(zz)
In [27]:
print(orbit_classification.shape)
print(orbit_classification.dtype)
orbit_classification
Out[27]:
In [28]:
orbit_classification[23,23]
Out[28]:
Making an image¶
PIL.Image.fromarray turns an array of uint8 values into a grayscale image.
In [31]:
from PIL import Image
In [40]:
imdata = 255*(1-orbit_classification.astype("uint8"))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[40]:
The filled Julia set of $g(z) = z^2-1$, often called "the basilica"
Several more examples¶
Assuming the cells above have been evaluated!
In [52]:
# We had three steps
# * Define function of z (e.g. g)
# * Turn it into a function that classifies orbits (e.g. has_bounded_orbit_under_g)
# * Vectorize that function so it can act on numpy arrays
# Let's combine them into one function, where you give a value `c` and it
# returns the vectorized orbit classifier for z^2+c
import time
def vectorized_classifier_factory(c,maxiter=100):
"""Returns a numpy ufunc that classifies orbits
of z^2+c as bounded or not"""
def classifier(a):
return bounded_orbit(lambda z:z*z+c, a, maxiter)
return np.vectorize(classifier)
In [70]:
size = 800
radius = 1.9
c = -1
maxiter = 500
jsname = "The Basilica"
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[70]:
In [54]:
size = 800
radius = 1.2
c = 0
maxiter = 500
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[54]:
In [81]:
size = 800
radius = 1.8
c = -0.52+0.53j
maxiter = 1000
jsname = None
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[81]:
In [76]:
size = 800
radius = 1.6
c = -0.635j
maxiter = 500
jsname = None
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[76]:
In [73]:
size = 800
radius = 2
c = -1.75488
maxiter = 500
jsname = "The Airplane"
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[73]:
In [75]:
size = 800
radius = 2
c = - 0.122561 + 0.744862j
maxiter = 500
jsname = "The Rabbit"
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[75]:
In [82]:
size = 800
radius = 2
c = 0.25
maxiter = 500
jsname = "The Cauliflower"
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[82]:
In [84]:
size = 800
radius = 2
c = -0.2+0.2j
maxiter = 500
jsname = None
# -- the rest of this cell is the same as all the other cells in this section --
t0 = time.time()
x = np.linspace(-radius,radius,size)
y = np.linspace(-radius,radius,size)
xx, yy = np.meshgrid(x,y)
zz = xx + 1j*yy
F = vectorized_classifier_factory(c,maxiter)
imdata = 255*(1-F(zz)).astype("uint8")
print("Julia set of z^2 + ({}), image radius {}, maxiter {}".format(
c,
radius,
maxiter
))
if jsname:
print("This Julia set is called '{}'".format(jsname))
print("Computation took {:.1f} seconds".format(time.time()-t0))
Image.fromarray(imdata).transpose(method=Image.FLIP_TOP_BOTTOM)
Out[84]:
