0

I'm trying to wrap my head around the perturbation theory paper by K.I. Martin. I've made a toy Python script to try it out using float64 as 'high precision' for the reference and calculating the rest of the pixels in float32. Unfortunately, the code I've written is doing what it should be. I have a few questions about the implementation.

My understanding is that the formula $X_{n+1} = X^2_n + X_0$ is to be iterated first in high precision. During iteration, I keep track of X_n over time. While doing this I also keep track of $A_{n+1}=2X_n A_n + 1$, $B_{n+1}=2X_n B_n + A^2_n$ and $C_{n+1}=2X_n C_n + 2A_n B_n$.

On this math stackexchange it is suggested that $A_0$, $B_0$ and $C_0$ should be initialised with 1, 0 and 0 respectively. However, these values are multiplied with $X$ which is a complex number. Should, for example, $A_0$ be initialised as 1+1i?

Then I iterate over all the pixels in low-precision. For each pixel Y, I calculate $\delta = \Delta_0$ using $\Delta_0 = Y_0 - X_0$. Then using this delta I recalculate Y, update $\Delta$ using $\Delta_n = A_n \delta + B_n \delta^2 + Cn \delta^3$ and the n calculate the new $Y_n$ using the new delta and the stored $X_n$. Then I check the new Y for escape as normal.

Here's the code I've written:

import numpy as np
import matplotlib.pyplot as plt
import copy

Complex multiplication


def cMul(a, b):
    return np.array([a[0]b[0] -  a[1]b[1], a[0]b[1] + a[1]  b[0]])


Complex addition


def cAdd(a, b):
    return np.array([a[0] + b[0], a[1] + b[1]])


Define a function to check if a point belongs to the Mandelbrot set


def mandelbrot(c, max_iter):
    Y = np.array([0.,0.])
    for n in range(max_iter):
        if np.sqrt(np.sum(np.square(Y))) > 2:
            return n
        Y = cAdd(cMul(Y,Y), c)


return max_iter




Define the range and resolution of the plot


x_min, x_max = -2, 2
y_min, y_max = -2, 2
resolution = 1000


Generate the grid of complex numbers


x = np.linspace(x_min, x_max, resolution, dtype=np.float32)
y = np.linspace(y_min, y_max, resolution, dtype=np.float32)
X, Y = np.meshgrid(x, y)
c = X + 1j * Y


Set the maximum number of iterations


max_iter = 100


Compute the Mandelbrot set


mandelbrot_set = np.zeros_like(c, dtype=np.float32)
for i in range(resolution):
    for j in range(resolution):
        x_coord = X[i,j]
        y_coord = Y[i,j]
        mandelbrot_set[i, j] = mandelbrot(np.array([x_coord,y_coord]), max_iter)


Plot the Mandelbrot set


plt.figure(figsize=(10, 8))
plt.imshow(mandelbrot_set, extent=(x_min, x_max, y_min, y_max), cmap='hot', origin='lower')
plt.colorbar(label='Iteration count')
plt.title('Mandelbrot Set')
plt.xlabel('Real axis')
plt.ylabel('Imaginary axis')
plt.show()


Keep track of X, A, B and C in high precision


def mandelbrot_history(max_iter, Xn, An, Bn, Cn, dtype):
    z = np.array([0.,0.], dtype=dtype)
    for n in range(max_iter):
        if np.sqrt(np.sum(np.square(z))) > 2:
            return n 
        z = cAdd(cMul(z,z), Xn[0])


    An1 = 2*cMul(Xn[-1],An[-1]) + 1
    Bn1 = 2*cMul(Xn[-1],Bn[-1]) + cMul(An[-1],An[-1])
    Cn1 = 2*cMul(Xn[-1],Cn[-1]) + 2*cMul(An[-1],Bn[-1])
    Xn.append(z)

    An.append(An1)
    Bn.append(Bn1)
    Cn.append(Cn1)

return max_iter




def perturb(delta0, Xn, An, Bn, Cn, max_iter):
    delta = delta0
    for n in range(max_iter):
        # obtain Y_n by adding X_n to delta_n
        Y = Xn[n] + delta
        if np.sqrt(np.sum(np.square(Y))) > 2:
            return n
        # calculate next delta (n+1)

        delta = cAdd(cAdd(cMul(An[n + 1], delta0), cMul(Bn[n + 1], cMul(delta0, delta0))), cMul(Cn[n + 1], cMul(delta0, cMul(delta0, delta0))))
    return n


init comes from https://math.stackexchange.com/questions/939270/perturbation-of-mandelbrot-set-fractal Neal Lawton comment


Xn = [np.array([0.,0.], dtype=np.float64)]
An = [np.array([1.,1.], dtype=np.float64)]
Bn = [np.array([0.,0.], dtype=np.float64)]
Cn = [np.array([0.,0.], dtype=np.float64)]


Obtain history of Xn, An, Bn, Cn in float 64/high precision


mandelbrot_history(max_iter, Xn, An, Bn, Cn, np.float64)


mandelbrot_set_orig = copy.deepcopy(mandelbrot_set)
for i in range(resolution):
    for j in range(resolution):
        x_coord = X[i,j]
        y_coord = Y[i,j]
        # calculate delta0
        delta0 = np.array([x_coord,y_coord],dtype=np.float32) - Xn[0]
        mandelbrot_set[i, j] = perturb(delta0, Xn, An, Bn, Cn, max_iter)


Plot the Mandelbrot set


plt.figure(figsize=(10, 8))
plt.imshow(mandelbrot_set, extent=(x_min, x_max, y_min, y_max), cmap='hot', origin='lower')
plt.colorbar(label='Iteration count')
plt.title('Perturbed Mandelbrot Set')
plt.xlabel('Real axis')
plt.ylabel('Imaginary axis')
plt.show()
Matige KunstIntelligentie
  • 1

0 Answers0