displaying Mandelbrot set in python using matplotlib.pyplot and numpy

Question

I am trying to get a plot of a Mandelbrot set and having trouble plotting the expected plot.

As I understand, the Mandelbrot set is made up of values c, which would converge if are iterated through the following equation z = z**2 + c. I used the initial value of z = 0.

Initially, I was getting a straight line. I look for solutions online to see where I went wrong. Using the following link in particular, I attempted to improve my code:

https://scipy-lectures.org/intro/numpy/auto_examples/plot_mandelbrot.html

Here is my improved code. I don't really understand the reason of using np.newaxis and why I am plotting the final z values that converge. Am I misunderstanding the definition of the Mandelbrot set?

# initial values 
loop = 50 # number of interations
div = 600 # divisions
# all possible values of c
c = np.linspace(-2,2,div)[:,np.newaxis] + 1j*np.linspace(-2,2,div)[np.newaxis,:] 
z = 0 
for n in range(0,loop):
      z = z**2 + c

plt.rcParams['figure.figsize'] = [12, 7.5]
z = z[abs(z) < 2] # removing z values that diverge 
plt.scatter(z.real, z.imag, color = "black" ) # plotting points
plt.xlabel("Real")
plt.ylabel("i (imaginary)")
plt.xlim(-2,2)
plt.ylim(-1.5,1.5)
plt.savefig("plot.png")
plt.show()

and got the following image, which looks closer to the Mandelbrot set than anything I got so far. But it looks more of a starfish with scattered dots around it. Image

For reference, here is my initial code before improvement:

# initial values 
loop = 50
div = 50
clist = np.linspace(-2,2,div) + 1j*np.linspace(-1.5,1.5,div) # range of c values 
all_results = []

for c in clist: # for each value of c
    z = 0 # starting point
    for a in range(0,loop): 
        negative = 0 # unstable

        z = z**2 + c 

        if np.abs(z) > 2: 
            negative +=1
        if negative > 2: 
            break

    if negative == 0:
        all_results.append([c,"blue"]) #converging
    else:
        all_results.append([c,"black"]) # not converging

Answer 1

Alternatively, with another small change to the code in the question, one can use the values of z to colorize the plot. One can store the value of n where the absolute value of the series becomes larger than 2 (meaning it diverges), and color the points outside the Mandelbrot set with it:

import pylab as plt
import numpy as np
# initial values 
loop = 50 # number of interations
div = 600 # divisions
# all possible values of c
c = np.linspace(-2,2,div)[:,np.newaxis] + 1j*np.linspace(-2,2,div)[np.newaxis,:] 
# array of ones of same dimensions as c
ones = np.ones(np.shape(c), np.int)
# Array that will hold colors for plot, initial value set here will be
# the color of the points in the mandelbrot set, i.e. where the series
# converges.
# For the code below to work, this initial value must at least be 'loop'.
# Here it is loop + 5
color = ones * loop + 5
z = 0
for n in range(0,loop):
      z = z**2 + c
      diverged = np.abs(z)>2
      # Store value of n at which series was detected to diverge.
      # The later the series is detected to diverge, the higher
      # the 'color' value.
      color[diverged] = np.minimum(color[diverged], ones[diverged]*n)

plt.rcParams['figure.figsize'] = [12, 7.5]
# contour plot with real and imaginary parts of c as axes
# and colored according to 'color'
plt.contourf(c.real, c.imag, color)
plt.xlabel("Real($c$)")
plt.ylabel("Imag($c$)")
plt.xlim(-2,2)
plt.ylim(-1.5,1.5)
plt.savefig("plot.png")
plt.show()

Answer 2

The plot doesn't look correct, because in the code in the question z (i.e. the iterated variable) is plotted. Iterating z = z*z + c, the Mandelbrot set is given by those real, imaginary part pairs of c, for which the series doesn't diverge. Hence the small change to the code as shown below gives the correct Mandelbrot plot:

import pylab as plt
import numpy as np
# initial values 
loop = 50 # number of interations
div = 600 # divisions
# all possible values of c
c = np.linspace(-2,2,div)[:,np.newaxis] + 1j*np.linspace(-2,2,div)[np.newaxis,:] 
z = 0 
for n in range(0,loop):
      z = z**2 + c

plt.rcParams['figure.figsize'] = [12, 7.5]
p = c[abs(z) < 2] # removing c values for which z has diverged 
plt.scatter(p.real, p.imag, color = "black" ) # plotting points
plt.xlabel("Real")
plt.ylabel("i (imaginary)")
plt.xlim(-2,2)
plt.ylim(-1.5,1.5)
plt.savefig("plot.png")
plt.show()