# Reports on Julia Set Activities

## Contents

# Reports on Julia Set Activities#

## Julia Set Activity 1#

Write a Python program to draw the Sierpinski gasket (perhaps by using binomial coefficients mod 2).

This turns out to be a very popular question, and there are several web tutorials and examples on the subject. For instance, here is a YouTube video and here is another (maybe more detailed) with a program that uses the Turtle package. Some links show *recursive* programs.

Sierpinski triangles are also part of this free book by Miller and Radum. Again we don’t know those people; these links are the result of a Google search. We just wanted to show you that there was a lot of stuff online about this.

The Wikipedia page is particularly good and has many references. It even connects to the Chaos Game, which we cover in our final unit. Compared to the animations of that page, and of some of the external links, our little code below is pretty basic. We used some ideas from Robin Truax’s video and we chose not to do a *full* Sierpinski triangle but rather a random sampling of it. For basic code, we think that the result is pretty cool.

```
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(2022)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_axes([0,0,1,1])
level = 7
scales = np.zeros(level)
for i in range(level):
scales[i] = 2**(-i-1)
numberpossible = 3**(level)
print("Using {} levels".format(level))
Nsample = 1729 # Ramanujan's number, why not
print("We will sample {} out of {} possible top vertices.".format(Nsample,numberpossible))
waterfall = np.transpose(np.array([[0,0], [1,0], [0,1]], dtype=float))
v0 = np.array([-0.5, -np.sqrt(3)/2]) # These directions give equilateral triangles
v1 = np.array([ 0.5, -np.sqrt(3)/2])
p0 = np.array([0,np.sqrt(3)/2])
for k in range(Nsample):
# Choose a random number base 3 with "level" trinary digits
pick = rng.integers(low=0, high=3, size=level) # exclusive top range. Sigh.
# Matrix whose columns are the allowed waterfall coordinate choices
A = np.transpose(np.array([waterfall[:,pick[i]] for i in range(level)], dtype=float))
xy = A.dot(scales) # Waterfall coordinates in binary
vertx = p0 + xy[0]*v0 + xy[1]*v1
lft = vertx + v0*2**(-level-1)
rgt = vertx + v1*2**(-level-1)
# Draw a triangle below the vertex
plt.plot([vertx[0],lft[0]], [vertx[1],lft[1]], 'k')
plt.plot([vertx[0],rgt[0]], [vertx[1],rgt[1]], 'k')
plt.plot([lft[0],rgt[0]], [lft[1],rgt[1]], 'k')
ax.set_xlim([-0.5,0.5 ])
ax.set_ylim([0, np.sqrt(3)/2])
plt.gca().set_aspect('equal')
plt.show()
```

```
Using 7 levels
We will sample 1729 out of 2187 possible top vertices.
```

## Julia Set Activity 2#

Explore pictures of the binomial coefficients mod 4 (and then consult Andrew Granville’s paper previously referenced). Here we didn’t have anything to add to that paper.

```
import math as mth
ans = mth.comb(15,8) # How to compute 15 choose 8
ans4 = ans % 4 # How to get its remainder mod 4
print(ans, ans4)
```

```
6435 3
```

## Julia Set Activity 3#

Investigate pictures (mod 2 or otherwise) of other combinatorial families of numbers, such as Stirling Numbers (both kinds). Try also “Eulerian numbers of the first kind” mod 3.

Talk about an open-ended activity. There is an *infinite* amount of information about Stirling numbers and other numbers; for instance the Digital Library of Mathematical Functions Chapter 26 lists a remarkable number of facts and links many references. As one fractal fact, computing the Stirling cycle numbers mod 2 gives a Sierpinski gasket, but tilted to one side. As just one thing, we chose to implement Stirling cycle numbers in Python in a way similar to the way Fibonacci numbers were naively implemented. For a reference, we used this handy paper but we really recommend that you look at Concrete Mathematics, by Graham, Knuth, and Patashnik. The recurrence relation we implemented is

with initial conditions \([n,0] = 1\) if \(n=0\) and \(0\) otherwise, and \([n,n] = 1\). This recurrence relation gives a triangular array. These numbers are the *Stirling cycle numbers* which count the number of permutations of \(n\) objects having \(k\) cycles.

```
nmax = 8
StirlingCycles = [[1],[0,1]]
for n in range(2,nmax):
prevrow = StirlingCycles[n-1]
nextrow = [0]
for k in range(1,n):
tmp = (n-1)*prevrow[k] + prevrow[k-1]
nextrow.append(tmp)
nextrow.append(1)
StirlingCycles.append(nextrow)
for n in range(nmax):
print(StirlingCycles[n])
```

```
[1]
[0, 1]
[0, 1, 1]
[0, 2, 3, 1]
[0, 6, 11, 6, 1]
[0, 24, 50, 35, 10, 1]
[0, 120, 274, 225, 85, 15, 1]
[0, 720, 1764, 1624, 735, 175, 21, 1]
```

```
# Now do the same mod 2
nmax = 21
StirlingCycles = [[1],[0,1]]
for n in range(2,nmax):
prevrow = StirlingCycles[n-1]
nextrow = [0]
for k in range(1,n):
tmp = ((n-1)*prevrow[k] + prevrow[k-1]) % 2 # The only change to the process
nextrow.append(tmp)
nextrow.append(1)
StirlingCycles.append(nextrow)
for n in range(nmax):
print(StirlingCycles[n])
```

```
[1]
[0, 1]
[0, 1, 1]
[0, 0, 1, 1]
[0, 0, 1, 0, 1]
[0, 0, 0, 1, 0, 1]
[0, 0, 0, 1, 1, 1, 1]
[0, 0, 0, 0, 1, 1, 1, 1]
[0, 0, 0, 0, 1, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 1, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1]
[0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1]
[0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1]
```

We’ve really just scratched the surface, here. Note that brackets are needed for modular computations: `(blah) % 2`

is needed, or if the modulus is complicated, `(blah) % (bug-ugly-modulus)`

.

## Julia Set Activity 4#

Write a Python program to animate Newton’s method for real functions and real initial estimates in general (the animated GIF above was produced by a Maple program, Student:-Calculus1:-NewtonsMethod, which is quite a useful model). This activity asks you to “roll your own” animation. At least, write Python code to draw the basic lines in the figure above.

Animation turns out to be complicated, when you think about the three things “Python”, “Jupyter notebook”, and “Jupyter Book”, together with all the different combinations of packages that are possible, and all the different kinds of machines to run them on. Here is a link to the documentation, and a blog post and some more information.

The documentation example alone has several new aspects to it. To help you get started, we do some preliminary programming. We consider (as with the figure) the function \(y = x^2 + 1\), and we write down explicitly the \(3\)-periodic orbit that we want to animate. If we start at \(x_0 = \cot \pi/7 \approx 2.07652139657233\), then it turns out that \(x_1 = \cot 2\pi/7\) and \(x_2 = \cot 4\pi/7\) and \(x_3 = \cot 8\pi/7 = \cot \pi/7 = x_0\) and we have a periodic orbit. So we want to draw a vertical dashed line from \((x_0,0)\) to \((x_0,x_0^2+1)\). Then draw a line from \((x_0, x_0^2+1)\) to \((x_1,0)\); then a vertical dashed line from \((x_1,0)\) to \((x_1, x_1^2+1)\); then a line from \((x_1,x_1^2+1)\) to \((x_2, 0)\); then a vertical dashed line to \((x_2,x_2^2+1)\); then a line from there to \((x_0,0)\), and then repeat. We can begin by simply drawing all those lines. Now, we need to figure out how to get “animate” to draw these lines in sequence.

We continued by following the ideas at the link above, but it turned out to be fussy to get the right packages loaded to make the animation work. We hope that you managed to find a combination that worked.

```
import numpy as np
import matplotlib.pyplot as plt
fig,ax = plt.subplots()
f = lambda x: x*x + 1
# x0 = np.cot(np.pi/7) # Pathetically, there is no cotangent in NumPy
cot = lambda x: 1.0/np.tan(x)
x0 = cot(np.pi/7)
y0 = f(x0)
x1 = cot(2*np.pi/7)
y1 = f(x1)
x2 = cot(4*np.pi/7)
y2 = f(x2)
x = np.arange(x2-0.5,x0+0.5,0.01)
y = f(x)
basecurve = ax.plot(x, y, color='red')
dash1 = ax.plot([x0,x0], [0,y0], linestyle='dashed', color='blue')
line1 = ax.plot([x0,x1], [y0,0], color='black' )
dash2 = ax.plot([x1,x1], [0,y1], linestyle='dashed', color='blue')
line2 = ax.plot([x1,x2], [y1,0], color='black')
dash2 = ax.plot([x2,x2], [0,y2], linestyle='dashed', color='blue')
line3 = ax.plot([x2,x0], [y2,0], color='black')
plt.show()
```

The following might be the first instance in this OER of a magic command.

```
%matplotlib inline
```

```
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML, Image
```

```
rc('animation', html='html5')
fig,ax = plt.subplots()
```

```
f = lambda x: x*x + 1
cot = lambda x: 1.0/np.tan(x)
x0 = cot(np.pi/7)
y0 = f(x0)
x1 = cot(2*np.pi/7)
y1 = f(x1)
x2 = cot(4*np.pi/7)
y2 = f(x2)
xdata = [[x0,x0],[x0,x1],[x1,x1],[x1,x2],[x2,x2],[x2,x0]]
ydata = [[0,y0], [y0,0], [0,y1], [y1,0], [0,y2], [y2, 0]]
c = ['blue', 'black']
sty = ['dashed', 'solid']
```

```
fig, ax = plt.subplots()
ax.set_xlim((x2-0.5, x0+0.5))
ax.set_ylim((0, y0+1))
```

```
x = np.arange(x2-0.5,x0+0.5,0.01)
y = f(x)
ax.plot(x, y, 'r')
line, = ax.plot([], [], lw=2)
```

```
def init():
line.set_data([],[])
return (line,)
def animate(i):
line.set_data( xdata[i%6], ydata[i%6]) # update the data.
line.set_color(c[i % 2])
line.set_linestyle(sty[i%2])
return (line,)
```

```
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=100, interval=50, blit=False)
```

```
anim.save('../Supplementary Material/Fractals and Julia Sets/animation.gif', writer='imagemagick', fps=1)
```

## Julia Set Activity 5#

Take a few moments and write out (maybe using pen and ink, in a nice notebook) a bunch of questions about this material. As always, don’t worry about *answers* just now. The goal is *questions*, and more questions.

What is the *dimension* of the boundary of the Newton fractal for \(z^3-1\)? How do you compute a dimension of a boundary, anyway? Which colouring scheme gives the most pleasing results? (We think the scheme in the notebook is pretty ugly; we hope you can do better). Which fractal is hardest to compute? Which is the simplest? Which is the most interesting to look at? What is the history of the fractal idea? Are there fractal *sounds*? Fractals seem to be good models for trees, rivers, clouds, blood vessel networks, algae mats, lichen, ice crystals growing in water as it freezes, frost on the windowpane, lots of things; do they really tell us anything interesting about those real-life things? There is lots of fractal art (again, see the Bridges Math Art site. How is it made? What makes it *artistic*? Is it art, if it’s computer-generated? (Some of the Deep Dream art makes us shudder) What are the open mathematical questions about fractals?

## Julia Set Activity 6#

Write your own code for computing Newton fractals, perhaps based on the code above (but at least improve the colour scheme).

```
import numpy as np
import sympy
import cmath
from sympy.abc import x
from sympy.utilities.lambdify import lambdify, implemented_function
import colorsys
import matplotlib
import matplotlib.pyplot as plt
from PIL import Image
np.seterr(all='ignore')
def gen_colours(n):
degree = 0
start = 60
space = 120
flag = 0
colours = []
for i in range(n):
if degree >= 360:
if flag == 0:
flag = 1
degree = start
start /= 2
else:
degree = start
start /= 2
space /= 2
colours.append(degree/360.0)
degree += space
return colours
def gen_shade(colours, n_shade):
col_len = len(colours)
shade = np.linspace(170, 0, n_shade)
list_shade = []
for c in colours:
temp = []
for s in shade:
if s <= 90.0:
t = colorsys.hsv_to_rgb(c, 0.1+s/100.0, 1.0)
temp.append(tuple(int(round(255*x)) for x in t))
else:
t = colorsys.hsv_to_rgb(c, 1.0, 1.0-(s - 90.0)/100.0)
temp.append(tuple(int(round(255*x)) for x in t))
list_shade.append(temp)
return list_shade
def colours(n, n_shade):
col = gen_colours(n)
return gen_shade(col, n_shade)
def newton(func, dfunc, z, tol, max_iter):
count = 0
while (abs(func(z)) > tol) and (abs(func(z)/dfunc(z))/abs(z)) and (count < max_iter):
if abs(func(z)) > 1e16:
return (None, max_iter)
z = z - func(z)/dfunc(z)
count += 1
if (count == max_iter):
return (None, max_iter)
else:
return(z, count)
def fractal(func, dfunc, bl, ur, size, tol, max_i, converge_tol):
roots = []
max_iter = 0
A = np.zeros((size[1], size[0], 3))
for i in range(size[1]):
y = i*(ur[1] - bl[1])/(size[1] - 1) + bl[1]
for j in range(size[0]):
x = j*(ur[0] - bl[0])/(size[0] - 1) + bl[0]
(root, count) = newton(func, dfunc, complex(x, y), tol, max_i)
if count > max_iter:
max_iter = count
if root:
alreadyFound = False
for test_root in roots:
if abs(test_root - root) < converge_tol:
root = test_root
alreadyFound = True
break
if not alreadyFound:
roots.append(root)
if root:
A[i][j][0] = root.real
A[i][j][1] = root.imag
A[i][j][2] = count
return (A, roots, max_iter)
def f_embed(M, **kwargs):
filename = kwargs.get('filename', None)
A = M[0]
(size_y, size_x, temp) = A.shape
roots = M[1]
n_shade = M[2]
C = colours(len(roots), n_shade)
img = Image.new("RGB", (size_x, size_y), (255, 255, 255))
for y in range(size_y):
for x in range(size_x):
if A[y][x][0]:
i = int(A[y][x][2])
if i > n_shade - 1:
i = n_shade - 1
img.putpixel((x, size_y - y - 1), C[roots.index(complex(A[y][x][0], A[y][x][1]))][i])
display(img)
if filename != None:
img.save(filename + ".png", "PNG")
def Newt_Fractal(func, bl, ur, size, **kwargs):
i = kwargs.get('iterationlimit', 25)
filename = kwargs.get('filename', None)
tolerance = kwargs.get('tolerance', 1e-6)
converge_tol = kwargs.get('converge_tolerance', 1e-3)
dfunc = kwargs.get('derivative', None)
if dfunc == None:
dfunc = lambdify(x, sympy.diff(func(x), x), np)
else:
dfunc = lambdify(x, dfunc(x), np)
func = lambdify(x, func(x), np)
if type(size) == int:
size = (size, size)
M = fractal(func, dfunc, bl, ur, size, tolerance, i, converge_tol)
print('This program found ' + str(len(M[1])) + ' roots.')
roots = kwargs.get('printroots', False);
if (roots == True):
print('The roots are')
print(M[1])
if filename == None:
f_embed(M)
else:
f_embed(M, filename = filename)
```

## Julia Set Activity 7#

Compute Newton fractals for several functions of your own choosing. Test your code on the function \(f(z) = z^{8}+4 z^{7}+6 z^{6}+6 z^{5}+5 z^{4}+2 z^{3}+z^{2}+z\) used above.

```
f = lambda z: z**8+3*z**7+5*z**6+5*z**5+4*z**4+2*z**3+z*z+z
Newt_Fractal(f, (-1.8, -1.4), (1.0, 1.4), 500)
```

```
This program found 8 roots.
```

## Julia Set Activity 8#

Compute Halley fractals for the same functions.

```
df = lambda z: 8*z**7+21*z**6+30*z**5+25*z**4+16*z**36*z**2+2*z+1
halley = lambda z: f(z)/sympy.sqrt(df(z))
Newt_Fractal(halley, (-1.8, -1.4), (1.0, 1.4), 500)
```

## Julia Set Activity 9#

Compute secant fractals for the same functions, using the \(x_1 = x_0 - f(x_0)/f'(x_0)\) rule to generate the needed second initial estimate. Try a different rule for generating \(x_1\) and see if it affects your fractals.

```
# A very short hacky program to draw the edges of a secant fractal
# RMC 2022.2.6
import numpy as np
from matplotlib import pyplot as plt
# We will take an N by N grid of initial estimates
N = 200 # 800 by 800 is a lot and it takes a few seconds to draw
x = np.linspace(-2,2,N)
y = np.linspace(-2,2,N)
F = np.zeros((N,N))
# Here is the function and its derivative whose zeros we are looking for
f = lambda x: x ** 3 - 2;
df = lambda x: 3*x**2 ;
# SirIsaac performs one Newton iteration
SirIsaac = lambda x, fx: x - fx/df(x);
secantIter = lambda x,y,fx,fy: x - fx*(x-y)/(fx-fy)
for k in range(N):
for i in range(N):
# We range over all initial estimates in the grid
z0 = x[i]+1j*y[k];
f0 = f(z0)
z1 = SirIsaac(z0, f0)
f1 = f(z1)
# Hard-wire in 20 iterations (maybe not enough)
for m in range(20):
z = secantIter(z1, z0, f1, f0)
z0 = z1
f0 = f1
z1 = z
if z1==z0:
break
f1 = f(z1)
# After twenty iterations we hope the iteration has settled down, except on
# the boundary between basins of attraction.
# The phase (angle) is a likely candidate for a unique identifier for the root
F[i,k] = np.angle(z)
# A magic incantation
X,Y = np.meshgrid(x, y)
# I admit to being quite puzzled as to why if we say X,Y in the call
# this contour plot has the real axis vertical, pointing up
# and the imaginary axis horizontal (not sure which way it's pointing)
# so I have just empirically input Y,X to get the graph looking right.
plt.figure(figsize=(10,10))
plt.contourf(Y,X, F, levels=[-3,-2,0,2,3], colors=['brown','red','black','yellow','black','blue','black'])
plt.gca().set_aspect('equal', adjustable='box')
```

## Julia Set Activity 10#

Try a few different values of “c” in the Mandelbrot example above, and generate your own “Julia sets”.

```
import numpy as np
from numpy.polynomial import Polynomial as Poly
N = 10001
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print(History[0])
here = 0
there = 1
c = [-0.73, 0, 1] # Try c=0.8 (Solution to Problem 1)
d = len(c)-1
while there <= N-d:
cc = c.copy()
cc[0] = c[0] - History[here]
p = Poly( cc );
rts = p.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [e.real for e in History]
y = [e.imag for e in History]
plt.scatter(x, y, s=0.5, marker=".")
plt.show()
```

```
N = 10001
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print(History[0])
here = 0
there = 1
c = [-0.2-0.8j, 0, 1] # Try c=0.2+0.8j (Solution to Problem 1)
d = len(c)-1
while there <= N-d:
cc = c.copy()
cc[0] = c[0] - History[here]
p = Poly(cc);
rts = p.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [e.real for e in History]
y = [e.imag for e in History]
plt.scatter(x, y, s=0.5, marker=".")
plt.show()
```

## Julia Set Activity 11#

These are not really Julia sets; they include too much of the history! Alter the program so that it plots only (say) the last half of the points computed; increase the number of points by a lot, as well. Compare your figure to (say) the Maple Julia set for c=1.2.

```
N = 64001
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print(History[0])
here = 0
there = 1
c = [-0.2-0.8j, 0, 1] # Try c=0.2+0.8j (Same as olution to Problem 1)
d = len(c)-1
while there <= N-d:
cc = c.copy()
cc[0] = c[0] - History[here]
p = Poly(cc);
rts = p.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
import math
x = [History[j].real for j in range(math.floor(63*N/64), N)] # Not just the last half, the last sixty-fourth
y = [History[j].imag for j in range(math.floor(63*N/64), N)]
plt.scatter(x, y, s=0.5, marker=".")
plt.show()
```

## Julia Set Activity 12#

Change the function F to be a different polynomial; find places where both F and F’ are zero (if any). If necessary, change your polynomial so that there is such a “critical point”. Start your iteration there, and go backwards. Plot your “Julia set”.

```
N = 3**10+1
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print( History[0])
here = 0
there = 1
# Newton's method on f(z) = z^3-1 gives z_{n+1} = (2z_n^3 + 1)/(3z_n^2)
p = [1, 0, 0, 2]
q = [0, 0, 3]
d = max(len(p), len(q)) - 1
while there <= N-d:
cc = Poly(p) - History[here]*Poly(q) # Python has polynomial arithmetic so this is simple
rts = cc.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [e.real for e in History]
y = [e.imag for e in History]
fractal = plt.figure(figsize=(8,8))
fractalax = plt.scatter(x, y, s=0.1, marker=".");
plt.xlim(-2,2);
plt.ylim(-2,2);
```

```
N = 3**10+1
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print(History[0])
here = 0
there = 1
# Halley's method on f(z) = z^3-1 gives z_{n+1} = z_n(z_n**3+2)/(2z_n^3+1)
p = [0, 2, 0, 0, 1]
q = [1, 0, 0, 2]
d = max(len(p), len(q)) - 1
while there <= N-d:
cc = Poly(p) - History[here]*Poly(q) # Python has polynomial arithmetic so this is simple
rts = cc.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [e.real for e in History]
y = [e.imag for e in History]
fractal = plt.figure(figsize=(8,8))
fractalax = plt.scatter(x, y, s=0.01, marker=".");
plt.xlim(-2,2);
plt.ylim(-2,2);
```

## Julia Set Activity 13#

Extend the program so that it works for *rational* functions F, say \(F(z) = p(z)/q(z)\). This means solving the polynomial equation \(p(z_n) - z_{n+1}q(z_n)=0\) for \(z_n\). Try it out on the rational functions you get from Newton iteration on polynomial (or rational!) functions; or on Halley iteration on polynomial functions. Try any of the that arise from the methods that you can find listed in Revisiting Gilbert Strang’s “A Chaotic Search for *i*”.

```
N = 90001
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print(History[0])
here = 0
there = 1
# Halley's method on a Mandelbrot polynomial gives
p = [0,0,0,-1,-9,-9,32,198,525,927,1236,1278,1014,570,192,28]
q = [1,3,9,26,69,186,474,948,1497,1874,1842,1404,766,252,36]
d = max(len(p), len(q)) - 1
while there <= N-d:
cc = Poly(p) - History[here]*Poly(q) # Python has polynomial arithmetic so this is simple
rts = cc.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [e.real for e in History]
y = [e.imag for e in History]
fractal = plt.figure(figsize=(8,8))
fractalax = plt.scatter(x, y, s=0.05, marker=".");
plt.xlim(-2,2);
plt.ylim(-2,2);
```

```
N = 15001
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print(History[0])
here = 0
there = 1
# Halley's method on the Wilkinson N=8 polynomial gives
p = [2209213440,-7144139520,8138672640,-1235703640,-6860019600,9045979866,-6258158928,2854960371,-921502440,216337751,-37195200,4642407,-409752,24255,-864,14]
q = [3622946688,-15347414304,29274362808,-33337014480,25349245110,-13634595702,5357600661,-1564677432,341763501,-55631772,6651225,-567000,32613,-1134,18]
d = max(len(p), len(q)) - 1
while there <= N-d:
cc = Poly(p) - History[here]*Poly(q) # Python has polynomial arithmetic so this is simple
rts = cc.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [e.real for e in History]
y = [e.imag for e in History]
fractal = plt.figure(figsize=(8,8))
fractalax = plt.scatter(x, y, s=0.5, marker=".");
plt.xlim(0,9);
plt.ylim(-1.75,1.75);
```

```
N = 180001
History = np.zeros(N,dtype=complex)
# History[0] is deliberately 0.0 for this example
#print( History[0 ])
here = 0
there = 1
# Halley's method on the Schroeder iteration for x**2 + 1 gives
p = [-1, 0, -6, 0, 3]
q = [0, 0, 0, 8]
d = max(len(p), len(q)) - 1
while there <= N-d:
cc = Poly(p) - History[here]*Poly(q) # Python has polynomial arithmetic so this is simple
rts = cc.roots();
#print(here, rts)
for j in range(d):
#print(j, rts[j], there)
History[there] = rts[j];
there += 1;
here += 1;
import matplotlib.pyplot as plt
x = [History[j].real for j in range(math.floor(44*N/45), N)] # Not just the last half, the last forty-fifth
y = [History[j].imag for j in range(math.floor(44*N/45), N)]
fractal = plt.figure(figsize=(8,8))
fractalax = plt.scatter(x, y, s=0.5, marker=".");
plt.xlim(-2,2);
plt.ylim(-1.,1.);
```

## Julia Set Activity 14#

Read the Wikipedia entry on Julia sets; it ought to be a little more intelligible now (but you will see that there are still lots of complications left to explain). One of the main items of interest is the theorem that states that the Fatou sets all have a *common boundary*. This means that if the number of components is \(3\) or more, then the Julia set (which is that boundary!) *must be a fractal*.

We still find that amazing.