--- 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:Solution10)= # Solution 10: Fractal Dimension of Julia Sets What is the fractal dimension of the Julia set $J_{N g(z)}$ corresponding to Newton's method for finding the roots of $g(z) = z^3 - 1$? This is the Julia set corresponding to the function \begin{gather*} f(z) = z - \frac{g(z)}{g'(z)} = \frac{2z^3 + 1}{3z^2} = \tfrac{2}{3}z + \frac{1}{3z^2}. \end{gather*} :::{margin} There might be a nice way of seeing this by showing that $z\mapsto f(z)$ contracts to the appropriate root if $z$ is close enough. For example, consider the root at 1. \begin{gather*} \abs{f'(z)} = \frac{2}{3}\frac{\abs{z^3 - 1}}{\abs{z}^3}. \end{gather*} One can probably argue that if the numerator $\delta = \abs{g(z)} = \abs{z^3 - 1}$ is sufficiently small, then the denominator $\abs{z}^3$ is sufficiently large to guarantee contraction. ::: The box-counting dimension $d$ relates the linear size $\epsilon$ of the cells $A = \epsilon^2$ to the number of cells in the Julia set $N(\epsilon) \propto \epsilon^{-d}$ so that \begin{gather*} d = -\lim_{\epsilon \rightarrow 0} \frac{\ln N}{\ln \epsilon}. \end{gather*} First we define some code to get the boxes that contain points from the Julia set. We do this as follows: 1. First we construct a grid with $N_{G} = 2^{G} + 1$ point in each direction. This allows use to recursively analyze nested subgrids of size $N_{n} = 2^{n}+1$ for $n\in \{0, 1, \cdots, N\}$. Combined with the extents $L$ of the region, we can compute $\epsilon = L/N_{N}$. *(In higher dimension, we take the geometric mean, i.e. $\epsilon = \sqrt{\d{x} \d{y}}$.)* 2. We then iterate each point on the grid with $z\mapsto f(z)$ until $\abs{g(z)} < 0.9$. *(See the discussion at the bottom for why this value works.)* 3. We assign to each point on the grid an integer $0$, $1$, or $2$, depending on which root is closest. 4. We treat these points as the boundaries of the cells, and claim that the cell contains points from $J$ is the 4 boundary points are not all attracted to the same attractor. *(Technically, we do this by taking the differences with `np.diff` and checking if any of the differences is non-zero.)* 5. We repeat this for all subgrids, counting how many cells $N(\epsilon)$ there are of "length" $\epsilon$, then fitting a line to the logarithms to extract the dimension. ```{code-cell} :tags: [hide-input] def g(z): return z**3 -1 def f(z): # Little trick using masked arrays to prevent 1/0 errors. # We fill these values with 0 so that the point at z=0 stays there rather than # diverging to np.inf where it raises floating point warning. return 2/3*z + np.ma.divide(1, 3*z**2).filled(0) roots = np.exp(2j*np.pi / 3 * np.arange(3)) assert np.allclose(g(roots), 0) def get_closest_root(z, maxiter=20, roots=roots): for n in range(maxiter): z = f(z) closest_root = np.argmin(abs(z[:, :, None] - roots[None, None, :]), axis=-1) return closest_root def get_J(attractor): """Return boolean array of cells containing points in J.""" corners = np.array([ attractor[:-1, :-1], attractor[:-1, 1:], attractor[1:, :-1], attractor[1:, 1:]]) inside = abs(np.diff(corners, axis=0)).max(axis=0) > 0 return inside def get_counts(z0, R, G=10, maxiter=30, ax=None): """Return `(epss, Ns)` box-counting the Newton fractal. See analyze for argument description. Returns ------- epss : list(float) List of the 1D cell lengths. Ns : list(int) Corresponding list of cell counts. """ N = 2**G + 1 x = np.linspace(z0.real-R, z0.real+R, N) y = np.linspace(z0.imag-R, z0.imag+R, N) z = x[:, None] + 1j*y[None, :] attractor = get_closest_root(z, maxiter=maxiter) epss = [] Ns = [] for n in 2**np.arange(0, G+1)[::-1]: # Reverse the order - large cells to small. J = get_J(attractor[::n, ::n]) x_, y_ = x[::n], y[::n] epss.append(np.sqrt(np.diff(x_).mean() * np.diff(y_).mean())) Ns.append(J.sum()) if ax is not None: ax.pcolormesh(x_, y_, J.T, alpha=0.1) if ax is not None: ax.set(aspect=1, xlabel="$x=\Re z$", ylabel="$y=\Im z$"); epss, Ns = np.array(epss), np.array(Ns) return epss, Ns def analyze(z0=0, R=1, G=10, maxiter=10, skip=5, plot=True): """Plot an extract the dimension. Arguments --------- z0 : complex Center of outer box. R : float "Radus" of the outer box, which will be a square +-R centered on z0. G : int Largest grid will have `2**G + 1` points. maxiter : int Maximum number of iterations. skip : int Number of largest subgrids to skip when computing the dimenson. plot : bool If `True`, then plot. """ if plot: fig, (ax, ax1) = plt.subplots(1, 2, figsize=(10, 5), gridspec_kw=dict(width_ratios=(2, 1))) else: ax = None epss, Ns = get_counts(z0=z0, R=R, G=G, maxiter=maxiter, ax=ax) dim, b = np.polyfit(np.log(1/epss)[skip:], np.log(Ns)[skip:], deg=1) if plot: ax1.loglog(1/epss, Ns, '+-') ax1.plot(1/epss, np.exp(np.polyval([dim, b], np.log(1/epss)))) ax1.set(xlabel="Inverse cell length $1/\epsilon$", ylabel="$N$", title=f"dim = {dim:}") return dim ``` Here we apply these to the Newton fractal for $g(z) = z^3 - 1$: ```{code-cell} %time analyze(z0=0, R=1.25) ``` It seems like the box-counting dimension is about $d = 1.4$. Let's try zooming in. ```{code-cell} %time analyze(z0=0, R=0.1) ``` Hmm. This looks strange. If we zoom in, the dimension seems to drop. After playing a bit, I found that the issues was my default choice of `maxiter`. Increasing this fixes the issue. :::{warning} You should always check that your results are converged! In this case, you should vary `maxiter` and `G`, making sure that your results do not depend on these parameters. If there is a significant dependence, then you must consider another method, or figure out the nature of this dependence so you can extrapolate to the limits of large `G` and `maxiter`. ::: ```{code-cell} %time analyze(z0=0, R=0.1, maxiter=100) ``` Zooming out has a similar issue. Increasing `maxiter` fixes this again. ```{code-cell} analyze(z0=0, R=100) %time analyze(z0=0, R=100, maxiter=100) ``` Apparently this Newton fractal is a **[multifractal][]** (see {cite}`Nauenberg:1989` for details). I expected that this meant that different regions of the Julia set might therefore have different dimensions. To check this, we zoom in off the axis: ```{code-cell} %time analyze(z0=-0.63+0.2575j, R=0.001, maxiter=100) ``` To within the precision of my analysis, I don't see a significant variation here. Have you found anything? Perhaps the notion of a multifractal is more complicated. ## Convergence Tests Here we perform some rudimentary convergence tests. ```{code-cell} args = dict(z0=0, R=1.5, plot=False) %time [analyze(maxiter=maxiter, **args) for maxiter in [10, 20, 30]] ``` ```{code-cell} args['maxiter'] = 50 %time [analyze(G=G, **args) for G in [9, 10, 11]] ``` ## Speed Improvements By plotting some contours, we show that if $\abs{g(z)} < 0.9$, then $z$ lies in the basin of attraction for the nearest root. This allows us to loosen our effective cutoff criterion. ```{code-cell} x = np.linspace(-1, 1.5, 256) y = np.linspace(-1.2, 1.2, 256) z = x[:, None] + 1j*y[None, :] attractor = get_closest_root(z) fig, ax = plt.subplots() ax.pcolormesh(x, y, attractor.T) cs = ax.contour(x, y, abs(g(z)).T, levels=[0.6, 0.8, 0.9, 1], colors='w') ax.clabel(cs, cs.levels, inline=True) ax.set(aspect=1, title="Contours of $g(z)$ on basins of attraction"); ``` Here we code a vectorized version... **Incomplete**. ## References [Newton's method]: [list of fractals by Hausdorff dimension]: [Multifractal]: