I wrote some code to render the Mandelbrot set with continuous coloring, but the bulbs on the period-2 blub are not colored correctly. regions that should not escape are colored as though the escape after 100 iterations, which happens to be my iteration limit.
What I tried

1. Increasing resolution and max iteration count
2. Increasing precision (np.float32 -> np.float64). Unfortunately, I could not finish execution because of the time complexity.

import numpy as np
import matplotlib.pyplot as plt

def mandelbrot(cmax, width, height, maxiter):
    real = np.linspace(-2, 0.5, width, dtype=np.float32)
    imag = np.linspace(0, cmax.imag, height // 2 + 1, dtype=np.float32)
    c_real, c_imag = np.meshgrid(real, imag)
    
    output = np.zeros((height // 2 + 1, width), dtype=np.uint16)
    z_real = np.zeros_like(c_real)
    z_imag = np.zeros_like(c_imag)
    mask = np.ones_like(c_real, dtype=bool)

    p = np.sqrt((c_real - 0.25)**2 + c_imag**2)
    cardioid = c_real < p - 2 * p**2 + 0.25
    period2_bulb = (c_real + 1)**2 + c_imag**2 < 0.0625
    mask[cardioid | period2_bulb] = False
    output[~mask] = maxiter - 1

    epsilon = 1e-10
    for i in range(maxiter):
        
        zr2 = z_real[mask] ** 2
        zi2 = z_imag[mask] ** 2
        
        z_imag_new = 2 * z_real[mask] * z_imag[mask] + c_imag[mask]
        z_real[mask] = zr2 - zi2 + c_real[mask]
        z_imag[mask] = z_imag_new
        
        diverged = zr2 + zi2 >= 4.0
        abs_z = zr2 + zi2 + epsilon  
        output[mask] = i + 1 - np.log(np.log(abs_z)) / np.log(2)
        mask[mask] = ~diverged

    output[output == maxiter - 1] = 0

    full_output = np.vstack([np.flipud(output), np.flipud(output[-2::-1, :])])
    return full_output


width, height, maxiter = 12000, 9000, 100
cmax = complex(0, 1)
mandelbrot_set = mandelbrot(cmax, width, height, maxiter)

plt.figure(figsize=(16, 12))
plt.imshow(mandelbrot_set, extent=[-2, 0.5, -1, 1], cmap="inferno", interpolation="bilinear")
plt.colorbar(label="Iterations to Divergence")
plt.title("Mandelbrot Set")
plt.xlabel("Re(c)")
plt.ylabel("Im(c)")
plt.show()

I wrote some code to render the Mandelbrot set with continuous coloring, but the bulbs on the period-2 blub are not colored correctly. regions that should not escape are colored as though the escape after 100 iterations, which happens to be my iteration limit.
What I tried

1. Increasing resolution and max iteration count
2. Increasing precision (np.float32 -> np.float64). Unfortunately, I could not finish execution because of the time complexity.

import numpy as np
import matplotlib.pyplot as plt

def mandelbrot(cmax, width, height, maxiter):
    real = np.linspace(-2, 0.5, width, dtype=np.float32)
    imag = np.linspace(0, cmax.imag, height // 2 + 1, dtype=np.float32)
    c_real, c_imag = np.meshgrid(real, imag)
    
    output = np.zeros((height // 2 + 1, width), dtype=np.uint16)
    z_real = np.zeros_like(c_real)
    z_imag = np.zeros_like(c_imag)
    mask = np.ones_like(c_real, dtype=bool)

    p = np.sqrt((c_real - 0.25)**2 + c_imag**2)
    cardioid = c_real < p - 2 * p**2 + 0.25
    period2_bulb = (c_real + 1)**2 + c_imag**2 < 0.0625
    mask[cardioid | period2_bulb] = False
    output[~mask] = maxiter - 1

    epsilon = 1e-10
    for i in range(maxiter):
        
        zr2 = z_real[mask] ** 2
        zi2 = z_imag[mask] ** 2
        
        z_imag_new = 2 * z_real[mask] * z_imag[mask] + c_imag[mask]
        z_real[mask] = zr2 - zi2 + c_real[mask]
        z_imag[mask] = z_imag_new
        
        diverged = zr2 + zi2 >= 4.0
        abs_z = zr2 + zi2 + epsilon  
        output[mask] = i + 1 - np.log(np.log(abs_z)) / np.log(2)
        mask[mask] = ~diverged

    output[output == maxiter - 1] = 0

    full_output = np.vstack([np.flipud(output), np.flipud(output[-2::-1, :])])
    return full_output


width, height, maxiter = 12000, 9000, 100
cmax = complex(0, 1)
mandelbrot_set = mandelbrot(cmax, width, height, maxiter)

plt.figure(figsize=(16, 12))
plt.imshow(mandelbrot_set, extent=[-2, 0.5, -1, 1], cmap="inferno", interpolation="bilinear")
plt.colorbar(label="Iterations to Divergence")
plt.title("Mandelbrot Set")
plt.xlabel("Re(c)")
plt.ylabel("Im(c)")
plt.show()

• python
• numpy
• math
• fractals
• mandelbrot

• How exactly were you expecting them to be colored? It is hard to tell why it doesn't go as you expected without knowing what you expected :D. – chrslg Commented Mar 22 at 22:22
• Matplotlib documentation actually has an example of plotting the Mandelbrot set matplotlib./stable/gallery/showcase/mandelbrot.html – RuthC Commented Mar 23 at 12:17
• 1 I find your computation quite strange. mask seems to be both a marker of what has diverged and what will never (or, to be accurate, the opposite). Which seems to me to be a complex intrication motivated by your colormap objective. But, well, that doesn't explain the problem, I think. Just, it makes it harder to think about it. And it seems to me that it would be better to 1st compute the "escape iteration", and the to deal with rendering afterward. – chrslg Commented Mar 23 at 14:21
• The other point, and that one, I believe, is what is your problem (tho it is hard to tell what is on purpose and what is not, since you still haven't described what result you expected), is that - np.log(np.log(abs_z)) / np.log(2). Since, abs_z is probably very small (otherwise it would have diverged before) almost everywhere (but on the boundary), that log(log()) is probably negative. So you end up with values greater than maxiter. And combine with my first remark, that means that after 100 iterations, everything is probably > maxiter – chrslg Commented Mar 23 at 14:34
• I mean everything, but what has diverged, and what you excluded in the first iterations outside the loop. So the whole border except bulb 1 and 2 (what you circled is not just bulb 2). The only reason why all the bulbs (but bulb-1 and main bulb 2) are not yellow, is because sometimes it is so small, that this is log(0). So -∞. So you end up with 0 in those areas (even before the output[output==maxiter-1]=0 – chrslg Commented Mar 23 at 14:42

1 Answer 1

The updating of the output and the mask inside the loop (iterations) is incorrect. We need to focus only on the points inside mask and escaped in current iteration.

def mandelbrot(cmax, width, height, maxiter):
    
    real = np.linspace(-2, 0.5, width, dtype=np.float32)
    imag = np.linspace(0, cmax.imag, height // 2 + 1, dtype=np.float32)
    c_real, c_imag = np.meshgrid(real, imag)
    
    output = np.zeros((height // 2 + 1, width), dtype=np.uint16)
    z_real = np.zeros_like(c_real)
    z_imag = np.zeros_like(c_imag)
    mask = np.ones_like(c_real, dtype=bool)

    #p = np.sqrt((c_real - 0.25)**2 + c_imag**2)
    #cardioid = c_real < p - 2 * p**2 + 0.25
    q = (c_real - 0.25)**2 + c_imag**2
    cardioid = q*(q + c_real - 1/4) <= 1/4*c_imag**2
    period2_bulb = (c_real + 1)**2 + c_imag**2 < 1/16 #0.0625
    mask[cardioid | period2_bulb] = False # these points never escape
    
    epsilon = 1e-10
    for i in range(maxiter):

        z_imag_new = 2 * z_real[mask] * z_imag[mask] + c_imag[mask]
        z_real[mask] = z_real[mask]**2 - z_imag[mask]**2 + c_real[mask]
        z_imag[mask] = z_imag_new
        
        diverged = z_real**2 + z_imag**2 >= 4.0 # escaped in this iteration
        abs_z = z_real[mask & diverged]**2 + z_imag[mask & diverged]**2 + epsilon
        output[mask & diverged] = i + 1 - np.log(np.log(abs_z)) / np.log(2) # only the points inside mask that escaped in the current iteration are need to be updated in the output
        mask[mask & diverged] = False # only the escaped points inside the mask are no longer needed to be considered
 
    #output[output == maxiter - 1] = 0
    full_output = np.vstack([np.flipud(output), np.flipud(output[-2::-1, :])])
    return full_output

Although this may not be a very efficient implementation and for more efficient implementations one could use numba jit , as suggested here.