--- jupytext: formats: ipynb,md:myst text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.6 kernelspec: display_name: Python 3 (math-583) language: python name: math-583 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os from pathlib import Path import logging logging.getLogger("matplotlib").setLevel(logging.CRITICAL) logging.getLogger('numba').setLevel(logging.CRITICAL) logging.getLogger('OMP').setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt ``` (sec:Assignment3)= # Assignment 3: Julia Set The goal of this assignment is to sample points in a Julia set $J_c$ for $z^2+c$ for various values of the parameter $c$. From these samples, you should use your program to compute a fractal dimension. ## Concrete Problem Statement For those of you who like a concrete problem statement, implement the following function: ``` def julia_sample(c, N): """Return a list of `N` points from the Julia set `J_c`.""" ``` This can then be used to estimate the fractal dimension of $J_c$ using box counting for example, as demonstrated in {ref}`sec:BSMBoxCounting`. ## Exploration To get a feel for the structure of Julia sets, you can play with the following code. This uses a compiled version of the escape-time algorithm which you should implement below. For now, you can play a bit: ```{code-cell} %matplotlib inline import numpy as np, matplotlib.pyplot as plt from math_583.mandelbrot import escape_time # Try changing the ranges here x = np.linspace(-2, 2, 1024) y = np.linspace(-2, 2, 1024) z = x[:, None] + 1j*y[None, :] # Try playing with different values of c c = 0.27334+0.00742j # Try changing maxiter - especially if you get deep. %time n = escape_time(c, z, maxiter=100) fig, ax = plt.subplots() mesh = ax.pcolormesh(x, y, np.log10(n).T, shading='auto') fig.colorbar(mesh, label=r"Escape time $\log_{10} n$") ax.set(xlabel=r"$x = \Re z$", ylabel=r"$y = \Im z$", title=f"Julia set $J_c$: $c={c}$", aspect=1); ``` :::{admonition} [XaoS][] The code here works, but the UX is a bit cludgy. You might want to instead download [XaoS][] and play with it. Note: the default there is to explore the [Mandelbrot set][] $M$, which we will do next assignment. The set $M$ provides a classification for the different Julia sets. To see the corresponding Julia set $J_c$, use the `m` ([Mandelbrot/Julia switching][]) or `j` ([Fast Julia preview mode][]) keyboard shortcuts . ::: [Xaos]: [Mandelbrot/Julia switching]: [Fast Julia preview mode]: ## Coding Suggestions ### Escape-Time Algorithm Here is a simple python version of the previous algorithm. It is very slow, but you can have complete control over it. Can you think how to make the algorithm faster (not simply by compiling)? If you have any substantial improvements, we can work with you to implement those in a compiled version. ```{code-cell} def julia_escape_time(c, R=4.0, maxiter=100, x_range=(-2, 2), y_range=(-2, 2), Nx=200, Ny=200): """Return `(x, y, res)` where `res` is the number iterations to "escape". Arguments --------- R : float Bailout radius. We say that a point has "escaped" once `abs(z) > R`. maxiter : int Maximum number of iterations. We say that a point has not escaped if `abs(z) < R` for `maxiter` iterations. Such points provide an approximation of the filled-in Julia set `K_c` and will have `n == maxiter`. c : complex Parameter defining the Julia set for $z^2+c$. x_range, y_range : (float, float) Range to plot: `x_range = (x[0], x[-1])` etc. Nx, Ny : int Number of pixels in each direction. `res.shape == (Nx, Ny)`. Results ------- x, y : float arrays 1D arrays of the cell boundaries. Note that these should have `Nx` and `Ny` elements because they are cell boundaries. See the code below for how to do this. res : 2D int array This is an array of integers expressing how many iterations it takes the initial point `z=x+j*y` to escape `abs(z) > R`. Note that `res.shape = (Nx, Ny)`. """ # Some code to get you started... but feel free to explore other techniques. x = np.linspace(*x_range, Nx) # These are the pixel centers, not the boundaries. y = np.linspace(*y_range, Ny) res = np.zeros((Nx, Ny), dtype=int) # Your code here. Possibly something based on the following loop, but this will be # slow. I recommend starting with a loop, but then trying to vectorize or optimize # your code later. for nx in range(Nx): for ny in range(Ny): ####################### # Calculate res[nx, ny] # # Your code here, but this works. z = x[nx] + 1j*y[ny] for n in range(maxiter + 1): z = z**2 + c if abs(z) > R: break res[nx, ny] = n return x, y, res ``` Here is how this function could be used to generate the plot above: ```{code-cell} c = 0.27334+0.00742j args = dict(c=c, maxiter=100, x_range=(-2, 2), y_range=(-2, 2), Nx=256, Ny=256) %time x, y, res = julia_escape_time(**args) fig, ax = plt.subplots() # Shading needed when x, y give cell centers instead of boundaries # The transpose is because we are not using meshgrid (image/MATLAB compatibility). mesh = ax.pcolormesh(x, y, np.log10(res).T, shading='auto') fig.colorbar(mesh, label=r"Escape time $\log_{10} n$") ax.set(xlabel=r"$x = \Re z$", ylabel=r"$y = \Im z$", title=f"Julia set $J_c$: $c={c}$", aspect=1); ``` ### BSM - Boundary Scanning Method The previous escape-time algorithm is useful for determining the filled-in Julia sets $K_c$, but for computing the fractal dimension, we are really interested in the Julia set itself $J_c = \partial K_c$, which is the boundary. (In most cases, $\mathrm{dim} K_c = 2$ which is uninteresting.) The goal here is to determine whether or not a cell $(x, y) \in [x_0, x_1] \times [y_0, y_1]$ contains an element $z \in J_c$ of the Julia set. ```{code-cell} def get_julia(c, x_range=(-2, 2), y_range=(-2, 2), Nx=200, Ny=200): """Return `(x, y, res)` where `res` is 1 if the corresponding cell contains points in $J_c$. Arguments --------- c : complex Parameter defining the Julia set for $z^2+c$. x_range, y_range : (float, float) Range to plot: `x_range = (x[0], x[-1])` etc. Nx, Ny : int Number of pixels in each direction. `res.shape == (Nx, Ny)`. Results ------- x, y : float arrays 1D arrays of the cell boundaries. Note that these should have `Nx+1` and `Ny+1` elements because they are cell boundaries. See the code below for how to do this. res : 2D bool array This is 1 where the corresponding cell contains points of the Julia set and 0 elsewhere. Note that `res.shape = (Nx, Ny)`. """ # Some code to get you started... but feel free to explore other techniques. x = np.linspace(*x_range, Nx+1) # Add one extra point for boundaries y = np.linspace(*y_range, Ny+1) res = np.zeros((Nx, Ny), dtype=bool) # Bool arrays save a little space. # Your code here. Possibly something based on the following loop, but this will be # slow. I recommend starting with a loop, but then trying to vectorize or optimize # your code later. for nx in range(Nx): for ny in range(Ny): x0, x1 = x[nx:nx+2] # Boundary of cell. y0, y1 = y[ny:ny+2] ##################### # Calculate f[nx, ny] # # Your code here. res[nx, ny] = True return x, y, res ``` Once you have this working, you should be able to apply the [BSM Box counting][] technique. [BSM Box Counting]: