In [1]:
# Draw some orbits
In [2]:
def draw_orbit(p, c, iters):
    """Compute list of points in orbit of point p under f(z)=z^2+c.
    Up to iters iterations."""
    list_orbit = [p] # list of points in orbit
    z = p # z will keep track of point in the orbit
    for n in range(iters):
        z = z*z + c  # apply one additional iterate
        list_orbit.append(z)
    return list_orbit
In [3]:
draw_orbit(complex(0.5,0),complex(0,0),5)
Out[3]:
[(0.5+0j),
 (0.25+0j),
 (0.0625+0j),
 (0.00390625+0j),
 (1.52587890625e-05+0j),
 (2.3283064365386963e-10+0j)]
In [4]:
# f(z)=z^2, orbit of i
list_plot(draw_orbit(complex(0,1),complex(0,0),5))
Out[4]:
In [5]:
draw_orbit(complex(1,0),complex(0,0),5)
Out[5]:
[(1+0j), (1+0j), (1+0j), (1+0j), (1+0j), (1+0j)]
In [6]:
# 1 is not an attracting fixed point; points nearby get pushed away
list_plot(draw_orbit(complex(1.1,0),complex(0,0),5))
Out[6]:
In [7]:
# 0 is a fixed point, and is attracting
list_plot(draw_orbit(complex(-0.11,-0.4),complex(0,0),6))
Out[7]:
In [8]:
# find fixed points for f(z)=z^2 + c
# solve z^2+c = z for z
# z^2-z+c = 0
def find_fp(c):
    p1 = (1+sqrt(1-4*c))/2
    p2 = (1-sqrt(1-4*c))/2
    return (p1,p2)
In [23]:
(p1,p2) = find_fp(1.0); (p1,p2)
Out[23]:
(0.500000000000000 + 0.866025403784439*I,
 0.500000000000000 - 0.866025403784439*I)
In [10]:
# confirm to make sure p1 is actually fixed
(0.500000000000000 + 0.866025403784439*I)^2 + 1.0
Out[10]:
0.499999999999999 + 0.866025403784439*I
In [11]:
# test whether f.p. is attracting
# for z^2 + 1
# compute multiplier abs(f'(p)) = abs(2p)
abs(2*p1)
# since this is >= 1, the fp is not attracting
Out[11]:
2.00000000000000
In [12]:
abs(2*p2)
# also not attracting
Out[12]:
2.00000000000000
In [13]:
# neither f.p. of z^2+1 is attracting
In [14]:
(p1,p2) = find_fp(0); (p1,p2)
Out[14]:
(1, 0)
In [15]:
(abs(2*p1), abs(2*p2))
Out[15]:
(2, 0)
In [16]:
# we already know 0 is an attracting fp for z^2
In [17]:
# try f(z) = z^2 + c for c small
(p1,p2)=find_fp(complex(0.1,0.06))
(abs(2*p1), abs(2*p2))
# p2 is attracting
Out[17]:
(1.7958190258468067, 0.2597567710742235)
In [18]:
# try f(z) = z^2 + c for c somewhat bigger
(p1,p2)=find_fp(complex(0.2499999,0))
(abs(2*p1), abs(2*p2))
# p2 is attracting
Out[18]:
(1.0006324555320427, 0.9993675444679573)
In [19]:
# try f(z) = z^2 + c for c somewhat bigger
(p1,p2)=find_fp(complex(0.25,0))
(abs(2*p1), abs(2*p2))
# p2 is no longer attracting
# correspond to cusp of main cardiod of mandlebrot
Out[19]:
(1, 1)
In [20]:
# try f(z) = z^2 + c for c somewhat bigger
(p1,p2)=find_fp(complex(0.26,0))
(abs(2*p1), abs(2*p2))
# p2 is no longer attracting
Out[20]:
(1.019803902718557, 1.019803902718557)
In [21]:
# outside of main cardiod, still have two f.p.
# but neither is attracting
(p1,p2)=find_fp(complex(0,1.5))
(abs(2*p1), abs(2*p2))
# p2 is no longer attracting
Out[21]:
(3.2933985694551016, 1.821826260461609)
In [22]:
# outside of main cardiod, still have two f.p.
# but neither is attracting
(p1,p2)=find_fp(complex(0,-1.3))
(abs(2*p1), abs(2*p2))
# p2 is no longer attracting
Out[22]:
(3.1374510494983605, 1.6573963762180182)
In [ ]:
# Upshot: main cardiod of Mandlebrot set 
# (the biggest heart-shaped region)
# parametrizes those f_c(z) that have an attracting fixed point.