numpy - Simple Mandelbrot set in Python -

i've looked @ other questions regarding this, can't seem figure out going wrong. goal "write program make image of mandelbrot set performing iteration values of c = x + iy on n × n grid spanning region −2 ≤ x ≤ 2 , −2 ≤ y ≤ 2. make density plot in grid points inside mandelbrot set colored black , outside colored white." being in mandelbrot set considered magnitude of z never greater 2 when iterated z' = z^2 + c.

my code produces grid entirely black 1 tiny notch of white.

from pylab import imshow,show,gray numpy import zeros,linspace  z = 0 + 0j n=100  m = zeros([n,n],int) xvalues = linspace(-2,2,n) yvalues = linspace(-2,2,n)  x in xvalues:     y in yvalues:         c = complex(x,y)         in range(100):             z = z*z + c             if abs(z) > 2.0:                 m[y,x] = 1                 break  imshow(m,origin="lower") gray() show() 

for future readers, here how new code ended looking:

from pylab import imshow,show,gray numpy import zeros,linspace  n=1000  m = zeros([n,n],int) xvalues = linspace(-2,2,n) yvalues = linspace(-2,2,n)  u,x in enumerate(xvalues):     v,y in enumerate(yvalues):         z = 0 + 0j         c = complex(x,y)         in range(100):             z = z*z + c             if abs(z) > 2.0:                 m[v,u] = 1                 break  imshow(m,origin="lower") gray() show() 

there couple of problems code.

firstly, using xvalues , yvalues index m, indices should pixel index integers in range 0..(n-1). numpy convert floats in xvalues , yvalues integers, resulting numbers in -2..2, there won't many pixels plotted, image tiny, , you'll wrapping due way negative indices work in python.

a simple way required pixel indices use built-in python function enumerate, there may way re-organize code numpy functions.

the other problem need reset z 0 each pixel. currently, code reuses whatever last z previous pixel, , if pixel in mandelbrot set z big.

here's repaired version of code. don't have pylab, wrote simple bitmap viewer using pil. can save them image file calling img.save(filename) inside show function; pil figure out correct file format filename extension.

import numpy np pil import image  def show(data):     img = image.frombytes('1', data.shape[::-1], np.packbits(data, 1))     img.show()  n = 100 maxiter = 100  m = np.zeros([n, n], np.uint8) xvalues = np.linspace(-2, 2, n) yvalues = np.linspace(-2, 2, n)  u, x in enumerate(xvalues):     v, y in enumerate(yvalues):         z = 0         c = complex(x, y)         in range(maxiter):             z = z*z + c             if abs(z) > 2.0:                 m[v, u] = 1                 break  show(m) 

here's output image:

---

of course, whenever find iterating on numpy array indices, that's sign you're doing wrong. main point of using numpy can perform operations on whole arrays @ once internally iterating on them @ c speed; above code might using plain python lists instead of numpy arrays.

here's version gets numpy of looping. uses more ram, it's 2.5 times faster previous version, , it's shorter.

this code uses several 2d arrays. c contains complex seed numbers, carry out core mandelbrot calculation in z, initialsed zeros. mask boolean array controls mandelbrot calculation needs performed. items set true, , on each iteration true items in mask correspond items in z have escaped mandelbrot set set false.

to test if point has escaped use z.real**2 + z.imag**2 > 4.0 rather abs(z) > 2.0, save little bit of time because avoids expensive square root calculation , abs function call.

we can use final value of mask plot mandelbrot set, make points in mandelbrot set black need invert values, can 1 - mask.

import numpy np pil import image  def show(data):     img = image.frombytes('1', data.shape[::-1], np.packbits(data, 1))     img.show()     img.save('mset.png')  n = 100 maxiter = 100  = np.linspace(-2, 2, n) c = + 1.j * a[:, none] z = np.zeros((n, n), np.complex128) mask = np.ones((n, n), np.bool)  in range(maxiter):     mask[mask] = z[mask].real**2 + z[mask].imag**2 < 4.0     z[mask] = z[mask]**2 + c[mask]  show(1 - mask)