---
jupytext:
  formats: ipynb,md:myst,py:light
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.13.8
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
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
```

(sec:Mandelbrot)=
Mandelbrot Set
==============
:::{margin}
If you have not already, please review {ref}`sec:Julia`.
:::
The Mandelbrot set $M$ is the set of all points $c$ such that the iteration $f: z
\mapsto z^2+c$ remains bounded starting from $z_0 = 0$.  It turns out that corresponding
Julia sets $J_c$ are connected and have finite measure.

We can approximate the Mandelbrot set with the same escape-time algorithm as for the
Julia sets, but instead of fixing $c$, we now let $c \in \mathbb{C}$ vary over the
plane, and fix the initial point $z_0 = 0+0\I$.  In this sense, $M$ describes what
happens to the origin of $K_c$, the filled-in Julia set: if $c\in M$, then $0 \in K_c$
and vice versa.  It turns out that if $c \in M$, then $J_c$ and $K_c$ are connected.

```{code-cell}
from math_583.mandelbrot import escape_time

def get_mandelbrot(c=-0.75, r=1.5, Nx=256*4, Ny=None, maxiter=100, z0=0, ax=None):
    """Draw M centered on `c` with radius `c+-r`."""
    x = np.linspace(c.real-r, c.real+r, Nx)
    if Ny is None:
        Ny = Nx
    ry = Ny*r/Nx
    y = np.linspace(c.imag-ry, c.imag+ry, Ny)
    z = x[:, None] + 1j*y[None, :]
    res = escape_time(c=z, z0=z0)
    return x, y, res

x, y, M = get_mandelbrot()
fig, ax = plt.subplots()
ax.pcolormesh(x, y, np.log(1+M).T, cmap='magma')
ax.set(aspect=1, xlabel=r"$x=\Re c$", ylabel="$y=\Im c$");

phi = np.linspace(0, 2*np.pi, 1000)
q = np.exp(1j*phi)
c1 = q*(2-q)/4
c2 = q/4 - 1

# From {cite}`Giarrusso:1995`
w = np.arcsinh((88-27*q)/80/np.sqrt(5))/3
c3s = [-7/4 - 20/9*(np.sinh(w + 2j*k*np.pi/3) - 1/4/np.sqrt(5))**2
       for k in [0, 1, 2]]
ax.plot(c1.real, c1.imag, '-', lw=0.5)
ax.plot(c2.real, c2.imag, '--', lw=0.5)
for c3 in c3s:
    ax.plot(c3.real, c3.imag, ':k', lw=0.5)
```

## A Surprise

In 1991, while exploring the point $c=-3/4$ at the neck between the main cardioid and
the largest bulb, David Boll stumbled on the following surprise.  He was trying to see
if this neck was infinitely thin, and so exploring the escape-times for points $c =
-3/4 + \I \epsilon$ along the imaginary axis when he found this series:

```{code-cell}
eps = 1/10**np.arange(8)

[escape_time(c=-3/4+eps*1j, R=4.0, maxiter=1000000000) for eps in eps]
```

He also approached the butt of the main cardioid $c=1/4 + \epsilon$ where he found this:

```{code-cell}
eps = 1/10**np.arange(13)

[escape_time(c=1/4+eps, R=4.0, maxiter=1000000000) for eps in eps]
```

This second example is simpler since the iteration is purely real (and was know
significantly earlier as a property of tangent bifurcations.  See {cite}`Peitgen:2004`
and references therein.

:::::{admonition} Question: Does a similar phenomenon happen at $c=(-1 + 3\sqrt{3}\I)/8$?

(Incomplete)

This point is on the not so trivial: if one simply parameterizes
\begin{gather*}
  c = \frac{-1 + 3\sqrt{3}\I}{8} + \epsilon
\end{gather*}
then the escape time fluctuates rapidly $c$ passes in and out of $M$.  Instead, we must move
along a path midway between the main cardioid $c_1$ and tertiary bulb $c_3$.  The main
cardioid $c_1$ and secondary bulb $c_2$ have simple forms (see {ref}`sec:Julia`),
but the tertiary bud requires solving a cubic {cite}`Giarrusso:1995`.

Suggestion: use the results of {cite}`Giarrusso:1995` (or derive your own) to find a
second-order approximation to both the main cardioid `c_1` and the tertiary bulb `c_3`
at $c$, and then approach long the midpoint between these.

\begin{gather*}
  c_1 = \frac{e^{\I \theta}(2-e^{\I\theta})}{4}, \qquad
  c_3 = ???
\end{gather*}

To find $c_3$ one can look for the solutions to the following two equations,
additionally factoring out the period 1 solutions:
\begin{gather*}
  f^{(3)}(z) = z = ((z^2 + c_3)^2 + c_3)^2 + c_3,\\
  \abs{f^{(3)}_{,z}(z)} = \abs{8z(z^2 + c_3)((z^2 + c_3)^2 + c_3)} = 1.
\end{gather*}
:::::

```{code-cell}
s = 2*np.pi/3
q = np.exp(1j*s)
c = q*(2-q)/4
fig, axs = plt.subplots(1, 2, figsize=(8, 4))

ax = axs[0]
x, y, M = get_mandelbrot(c, r=1)
ax.pcolormesh(x, y, np.log(1+M).T, cmap='magma')
ax.plot([c.real], [c.imag], 'xg')
ax.set(aspect=1, xlabel=r"$x=\Re c$", ylabel="$y=\Im c$");

ax = axs[1]
x, y, M = get_mandelbrot(c, r=0.1, maxiter=2000, Nx=1024, Ny=256)
ax.pcolormesh(x, y, np.log(1+M).T, cmap='magma')
ax.plot([c.real], [c.imag], 'xg')
ax.axhline([c.imag], ls=':', c='g')
ax.set(aspect=1, xlabel=r"$x=\Re c$", ylabel="$y=\Im c$");
```


## Question:
### References

* {cite}`Peitgen:2004` discusses this ([see 1st ed. through WSU
  Library](https://ebookcentral.proquest.com/lib/wsu/reader.action?docID=3085512&ppg=892)).
  
* [Pi and fractal sets: The Mandelbrot set - Dave Boll - Gerald
   Edgar](http://www.pi314.net/eng/mandelbrot.php).

[Numba]: <https://numba.pydata.org/>