alexandre de zotti draws julia sets of quadratic polynomials and exponential maps more specifically this iterates the function a fixed number of times then plots whether the absolute value of the last iterate is greater than a fixed threshold named escape radius for the exponential map this is not really an escape radius but rather a convenient way to approximate the julia set with bounded orbits the examples presented here are the cauliflower julia set see e g https en wikipedia orgwikifile juliaz22b0 25 png other examples from https en wikipedia orgwikijuliaset an exponential map julia set ambiantly homeomorphic to the examples in https www math univtoulouse frcheritatgaliigalery html and https ddd uab catpubpubmat02141493v43n102141493v43n1p27 pdf remark some overflow runtime warnings are suppressed this is because of the way the iteration loop is implemented using numpy s efficient computations overflows and infinites are replaced after each step by a large number evaluate ez c evalexponential0 0 1 0 absevalexponential1 numpy pi1 j 1e15 true absevalexponential1 j 011 j 1e15 true evalquadraticpolynomial0 2 4 evalquadraticpolynomial1 1 0 roundevalquadraticpolynomial1 j 0 imag 1 roundevalquadraticpolynomial1 j 0 real 0 create a grid of complex values of size nbpixelsnbpixels with real and imaginary parts ranging from windowsize to windowsize inclusive returns a numpy array preparegrid1 3 array1 1 j 1 0 j 1 1 j 0 1 j 0 0 j 0 1 j 1 1 j 1 0 j 1 1 j iterate the function evalfunction exactly nbiterations times the first argument of the function is a parameter which is contained in functionparams the variable z0 is an array that contains the initial values to iterate from this function returns the final iterates iteratefunctionevalquadraticpolynomial 0 3 numpy array0 1 2 shape 3 numpy rounditeratefunctionevalquadraticpolynomial 0 3 numpy array0 1 20 0j numpy rounditeratefunctionevalquadraticpolynomial 0 3 numpy array0 1 21 10j numpy rounditeratefunctionevalquadraticpolynomial 0 3 numpy array0 1 22 2560j plots of whether the absolute value of zfinal is greater than the value of escaperadius adds the functionlabel and functionparams to the title showresults 80 0 1 numpy array0 1 5 4 2 1 1 2 1 1 3 ignore some overflow and invalid value warnings ignoreoverflowwarnings evaluate e z c eval_exponential 0 0 1 0 abs eval_exponential 1 numpy pi 1 j 1e 15 true abs eval_exponential 1 j 0 1 1 j 1e 15 true eval_quadratic_polynomial 0 2 4 eval_quadratic_polynomial 1 1 0 round eval_quadratic_polynomial 1 j 0 imag 1 round eval_quadratic_polynomial 1 j 0 real 0 create a grid of complex values of size nb_pixels nb_pixels with real and imaginary parts ranging from window_size to window_size inclusive returns a numpy array prepare_grid 1 3 array 1 1 j 1 0 j 1 1 j 0 1 j 0 0 j 0 1 j 1 1 j 1 0 j 1 1 j iterate the function eval_function exactly nb_iterations times the first argument of the function is a parameter which is contained in function_params the variable z_0 is an array that contains the initial values to iterate from this function returns the final iterates iterate_function eval_quadratic_polynomial 0 3 numpy array 0 1 2 shape 3 numpy round iterate_function eval_quadratic_polynomial 0 3 numpy array 0 1 2 0 0j numpy round iterate_function eval_quadratic_polynomial 0 3 numpy array 0 1 2 1 1 0j numpy round iterate_function eval_quadratic_polynomial 0 3 numpy array 0 1 2 2 256 0j plots of whether the absolute value of z_final is greater than the value of escape_radius adds the function_label and function_params to the title show_results 80 0 1 numpy array 0 1 5 4 2 1 1 2 1 1 3 ignore some overflow and invalid value warnings ignore_overflow_warnings see file header for explanations
import warnings
from collections.abc import Callable
from typing import Any

import numpy
from matplotlib import pyplot

c_cauliflower = 0.25 + 0.0j
c_polynomial_1 = -0.4 + 0.6j
c_polynomial_2 = -0.1 + 0.651j
c_exponential = -2.0
nb_iterations = 56
window_size = 2.0
nb_pixels = 666


def eval_exponential(c_parameter: complex, z_values: numpy.ndarray) -> numpy.ndarray:
    
    return numpy.exp(z_values) + c_parameter


def eval_quadratic_polynomial(
    c_parameter: complex, z_values: numpy.ndarray
) -> numpy.ndarray:
    
    return z_values * z_values + c_parameter


def prepare_grid(window_size: float, nb_pixels: int) -> numpy.ndarray:
    
    x = numpy.linspace(-window_size, window_size, nb_pixels)
    x = x.reshape((nb_pixels, 1))
    y = numpy.linspace(-window_size, window_size, nb_pixels)
    y = y.reshape((1, nb_pixels))
    return x + 1.0j * y


def iterate_function(
    eval_function: Callable[[Any, numpy.ndarray], numpy.ndarray],
    function_params: Any,
    nb_iterations: int,
    z_0: numpy.ndarray,
    infinity: float | None = None,
) -> numpy.ndarray:
    

    z_n = z_0.astype("complex64")
    for _ in range(nb_iterations):
        z_n = eval_function(function_params, z_n)
        if infinity is not None:
            numpy.nan_to_num(z_n, copy=False, nan=infinity)
            z_n[abs(z_n) == numpy.inf] = infinity
    return z_n


def show_results(
    function_label: str,
    function_params: Any,
    escape_radius: float,
    z_final: numpy.ndarray,
) -> None:
    

    abs_z_final = (abs(z_final)).transpose()
    abs_z_final[:, :] = abs_z_final[::-1, :]
    pyplot.matshow(abs_z_final < escape_radius)
    pyplot.title(f"Julia set of ${function_label}$, $c={function_params}$")
    pyplot.show()


def ignore_overflow_warnings() -> None:
    
    warnings.filterwarnings(
        "ignore", category=RuntimeWarning, message="overflow encountered in multiply"
    )
    warnings.filterwarnings(
        "ignore",
        category=RuntimeWarning,
        message="invalid value encountered in multiply",
    )
    warnings.filterwarnings(
        "ignore", category=RuntimeWarning, message="overflow encountered in absolute"
    )
    warnings.filterwarnings(
        "ignore", category=RuntimeWarning, message="overflow encountered in exp"
    )


if __name__ == "__main__":
    z_0 = prepare_grid(window_size, nb_pixels)

    ignore_overflow_warnings()  

    nb_iterations = 24
    escape_radius = 2 * abs(c_cauliflower) + 1
    z_final = iterate_function(
        eval_quadratic_polynomial,
        c_cauliflower,
        nb_iterations,
        z_0,
        infinity=1.1 * escape_radius,
    )
    show_results("z^2+c", c_cauliflower, escape_radius, z_final)

    nb_iterations = 64
    escape_radius = 2 * abs(c_polynomial_1) + 1
    z_final = iterate_function(
        eval_quadratic_polynomial,
        c_polynomial_1,
        nb_iterations,
        z_0,
        infinity=1.1 * escape_radius,
    )
    show_results("z^2+c", c_polynomial_1, escape_radius, z_final)

    nb_iterations = 161
    escape_radius = 2 * abs(c_polynomial_2) + 1
    z_final = iterate_function(
        eval_quadratic_polynomial,
        c_polynomial_2,
        nb_iterations,
        z_0,
        infinity=1.1 * escape_radius,
    )
    show_results("z^2+c", c_polynomial_2, escape_radius, z_final)

    nb_iterations = 12
    escape_radius = 10000.0
    z_final = iterate_function(
        eval_exponential,
        c_exponential,
        nb_iterations,
        z_0 + 2,
        infinity=1.0e10,
    )
    show_results("e^z+c", c_exponential, escape_radius, z_final)