A document from MCS 275 Spring 2022, instructor David Dumas. You can also get the notebook file.

Filled Julia sets

MCS 275 Spring 2022 Lecture 27 - David Dumas

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]:
[2, 4, 16, 256, 65536, 4294967296]
In [9]:
orbit(f,1/2,5)
Out[9]:
[0.5, 0.25, 0.0625, 0.00390625, 1.52587890625e-05, 2.3283064365386963e-10]
In [10]:
orbit(g,2,10)
# a=2 is not in K_g
Out[10]:
[2,
 3,
 8,
 63,
 3968,
 15745023,
 247905749270528,
 61457260521381894004129398783,
 3776994870793005510047522464634252677140721938309041881088,
 14265690253996672387291309349232388828298289458234016200317876247121873778287073518355813134107244701354409532063743,
 203509918422975643131625101557493482417516581321177887486174612349689723637228095885608636916042736102424433936617827686729959884103612889121780011527849662844136145473332830301242513425767433130492258922182178369327277600615170048]
In [11]:
orbit(g,1.5,30)
# a=1.5 seems to be in K_g
Out[11]:
[1.5,
 1.25,
 0.5625,
 -0.68359375,
 -0.5326995849609375,
 -0.7162311521824449,
 -0.4870129366434074,
 -0.7628183995419644,
 -0.418108089320236,
 -0.8251856256449815,
 -0.31906868322890036,
 -0.8981951753825757,
 -0.19324542691946411,
 -0.962656204974714,
 -0.07329303102368134,
 -0.9946281316033617,
 -0.010714879823205736,
 -0.9998851913503742,
 -0.00022960411822547933,
 -0.9999999472819489,
 -1.0543609951785271e-07,
 -0.9999999999999889,
 -2.220446049250313e-14,
 -1.0,
 0.0,
 -1.0,
 0.0,
 -1.0,
 0.0,
 -1.0,
 0.0]
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]:
[(0.5+0.5j),
 (-1+0.5j),
 (-0.25-1j),
 (-1.9375+0.5j),
 (2.50390625-1.9375j),
 (1.5156402587890625-9.70263671875j),
 (-92.84399390197359-29.41141365468502j),
 (7753.975950502717+5461.3462200079975j),
 (30297839.50617887+84694294.49462186j),
 (-6255164441199562+5132108283374204j),
 (1.2788546754169306e+31-6.42043624850961e+31j)]
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]:
False
In [15]:
bounded_orbit(f,1/2)
Out[15]:
True
In [16]:
bounded_orbit(g,1.5)
Out[16]:
True
In [17]:
bounded_orbit(g,0.5+0.5j)
Out[17]:
False

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]:
[array([[1, 2, 3],
        [1, 2, 3],
        [1, 2, 3]]),
 array([[ 8,  8,  8],
        [ 9,  9,  9],
        [10, 10, 10]])]
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

(500, 500)
Out[20]:
array([[-2.        , -1.99198397, -1.98396794, ...,  1.98396794,
         1.99198397,  2.        ],
       [-2.        , -1.99198397, -1.98396794, ...,  1.98396794,
         1.99198397,  2.        ],
       [-2.        , -1.99198397, -1.98396794, ...,  1.98396794,
         1.99198397,  2.        ],
       ...,
       [-2.        , -1.99198397, -1.98396794, ...,  1.98396794,
         1.99198397,  2.        ],
       [-2.        , -1.99198397, -1.98396794, ...,  1.98396794,
         1.99198397,  2.        ],
       [-2.        , -1.99198397, -1.98396794, ...,  1.98396794,
         1.99198397,  2.        ]])
In [21]:
yy
Out[21]:
array([[-2.        , -2.        , -2.        , ..., -2.        ,
        -2.        , -2.        ],
       [-1.99198397, -1.99198397, -1.99198397, ..., -1.99198397,
        -1.99198397, -1.99198397],
       [-1.98396794, -1.98396794, -1.98396794, ..., -1.98396794,
        -1.98396794, -1.98396794],
       ...,
       [ 1.98396794,  1.98396794,  1.98396794, ...,  1.98396794,
         1.98396794,  1.98396794],
       [ 1.99198397,  1.99198397,  1.99198397, ...,  1.99198397,
         1.99198397,  1.99198397],
       [ 2.        ,  2.        ,  2.        , ...,  2.        ,
         2.        ,  2.        ]])
In [22]:
zz = xx + 1j*yy
In [23]:
zz # 50 x 50 matrix
Out[23]:
array([[-2.        -2.j        , -1.99198397-2.j        ,
        -1.98396794-2.j        , ...,  1.98396794-2.j        ,
         1.99198397-2.j        ,  2.        -2.j        ],
       [-2.        -1.99198397j, -1.99198397-1.99198397j,
        -1.98396794-1.99198397j, ...,  1.98396794-1.99198397j,
         1.99198397-1.99198397j,  2.        -1.99198397j],
       [-2.        -1.98396794j, -1.99198397-1.98396794j,
        -1.98396794-1.98396794j, ...,  1.98396794-1.98396794j,
         1.99198397-1.98396794j,  2.        -1.98396794j],
       ...,
       [-2.        +1.98396794j, -1.99198397+1.98396794j,
        -1.98396794+1.98396794j, ...,  1.98396794+1.98396794j,
         1.99198397+1.98396794j,  2.        +1.98396794j],
       [-2.        +1.99198397j, -1.99198397+1.99198397j,
        -1.98396794+1.99198397j, ...,  1.98396794+1.99198397j,
         1.99198397+1.99198397j,  2.        +1.99198397j],
       [-2.        +2.j        , -1.99198397+2.j        ,
        -1.98396794+2.j        , ...,  1.98396794+2.j        ,
         1.99198397+2.j        ,  2.        +2.j        ]])

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]:
True
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

(500, 500)
bool
Out[27]:
array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])
In [28]:
orbit_classification[23,23]
Out[28]:
False

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)

Julia set of z^2 + (-1), image radius 1.9, maxiter 500
This Julia set is called 'The Basilica'
Computation took 5.0 seconds
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)

Julia set of z^2 + (0), image radius 1.2, maxiter 500
Computation took 24.1 seconds
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)

Julia set of z^2 + ((-0.52+0.53j)), image radius 1.8, maxiter 1000
Computation took 8.5 seconds
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)

Julia set of z^2 + ((-0-0.635j)), image radius 1.6, maxiter 500
Computation took 8.2 seconds
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)

Julia set of z^2 + (-1.75488), image radius 2, maxiter 500
This Julia set is called 'The Airplane'
Computation took 0.6 seconds
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)

Julia set of z^2 + ((-0.122561+0.744862j)), image radius 2, maxiter 500
This Julia set is called 'The Rabbit'
Computation took 4.0 seconds
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)

Julia set of z^2 + (0.25), image radius 2, maxiter 500
This Julia set is called 'The Cauliflower'
Computation took 8.6 seconds
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)

Julia set of z^2 + ((-0.2+0.2j)), image radius 2, maxiter 500
Computation took 8.2 seconds
Out[84]: