In [1]:
import numpy # will use Numpy package to represent grid of points
In [2]:
a = numpy.zeros((3,3)); a # makes a 3x3 array of all zeros
Out[2]:
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
In [3]:
a[0,0] # can access elements similarly to lists
Out[3]:
0.0
In [4]:
a[1,0] = 5 # can assign values to particular elements
In [5]:
a
Out[5]:
array([[0., 0., 0.],
       [5., 0., 0.],
       [0., 0., 0.]])
In [7]:
height = 3
width = 3
In [8]:
test_array = numpy.zeros((height,width))
In [9]:
test_array
Out[9]:
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
In [10]:
test_array[1,1]=1; test_array[2,1]=0.5; test_array
Out[10]:
array([[0. , 0. , 0. ],
       [0. , 1. , 0. ],
       [0. , 0.5, 0. ]])
In [11]:
matrix_plot(test_array) # Sage comand to plot array of numbers between 0 and 1
Out[11]:

Julia sets

Below is the code to draw Julia set of quadratic polynomial $f(z)=z^2+c$.

The Julia set of a polynomial $f$ is defined as the boundary of the escaping set $$\{z: \lim_{n\to\infty} f^{\circ n}(z) = \infty\}$$

In [3]:
# Code to draw Julia set of quadratic polynomial f(z)=z^2+c
# (Should only use with |c|<2)

import numpy
height = 600  # height in pixels of image
width = 600

c=complex(0.3,0.4) # python complex number 0.3 + 0.4i
          # polynomial is z^2 + c
max_iter=30 # max number of iterations 
max_abs=2 # if |z|>2 then the iterates f(z),f(f(z)),... escape to infinity

xmin=-2; xmax=2  #smallest, largest x-coordinates to check
ymin=-2; ymax=2
xinc=(xmax-xmin)/width   # width of each pixel
yinc=(ymax-ymin)/height

julia = numpy.zeros((height,width)) # initialize a numpy array of all zeros

for x in range(width):
    for y in range(height):
        z = complex((xmin + xinc*x),(ymin + yinc*y))
        iters = 0 # count of number of iterations
        while iters<max_iter:
            z = z*z + c
            if abs(z)>max_abs: 
                # in this case, known escape by Claim in class
                julia[y,x] = 1  # color black in this case
                    # need this order of x,y to plot correctly
                break
            else: 
                iters += 1
matrix_plot(julia, origin='lower')
Out[3]: