--- 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:Assignment8)= # Assignment 8: Feigenbaum Constant Revisited :::{margin} Note: our notation is slightly different from that in {cite}`Schroeder:1991`. We use $r_{n}$ and $R_{n}$ to represent the bifurcation and superstable points for period-$n$. {cite}`Schroeder:1991` uses these notations for period $2^n$. ::: Let us revisit the computation of the Feigenbaum constant. From {ref}`sec:Assignment6` you should have been able to at least estimate this from the location of the bifurcation points $r_{n}$ where the period doubles to $2^n$, but computing these is a little tricky. {cite}`Schroeder:1991` suggests a better approach in Chapter 12 section **Self Similarity in the Logistic Parabola** using the values of $r=R_{2^n}$ that give **superstable** orbits of period $2^n$. ## Stability Consider a fixed-point $x$ such that $f(x) = x$. Perturbing this we have \begin{gather*} f(x+\epsilon) \approx f(x) + \epsilon f'(x) = x + \epsilon f'(x). \end{gather*} Thus, the fixed-point is attractive if $\abs{f'(x)} < 1$, repulsive if $\abs{f'(x)} > 1$, and indeterminant if $\abs{f'(x)} = 1$. A **stable** fixed point will **attract** nearby points. A fixed-point is **superstable** if $f'(x) = 0$. This notion can be applied to orbits of period $n$, i.e. fixed points of $x = f^{(n)}(x)$. ## Question 1: Find $R_1$ Find parameter values $R_{1}$ that give rise to a superstable fixed-point of period $1$. Demonstrate numerically that this point is indeed exhibits very fast convergence. ```{code-cell} R1 = 1.5 # WRONG! Put in the correct value ``` ```{code-cell} %matplotlib inline import numpy as np, matplotlib.pyplot as plt from math_583.plotting import CobWeb plt.rcParams['figure.figsize'] = (5, 3) def f(x, r, d=0): """Logistic growth function (or derivative)""" if d == 0: return r*x*(1-x) elif d == 1: return r*(1-2*x) elif d == 2: return -2*r else: return 0 x = np.linspace(0, 1) fig, axs = plt.subplots(1, 2, figsize=(8, 2)) ax = axs[0] ax.plot(x, x, '-k') ax.plot(x, f(x, r=R1), label=f"$r={R1}$") ax.legend() ax.set(xlabel="$x$", ylabel="$f_r(x)$"); CobWeb(f).cobweb(r=R1, ax=axs[1]); ``` ```{code-cell} # Maybe this will help. It should converge in about 7 steps, with extremely # rapid convergence once it gets close if you have the correct value for R1 r = R1 xs = [0, 0.2] while abs(np.diff(xs[-2:])) > 1e-12 and len(xs) < 100: x = xs[-1] xs.append(r*x*(1-x)) print(x) ``` :::{solution} We look for solutions to the equations \begin{gather*} x = f(x) = rx(1-x), \qquad 0 = f'(x) = r(1-2x),\\ x = \frac{1}{2}, \qquad R_1 = \frac{1}{1-x} = 2. \end{gather*} ::: ## Solving for $R_n$ You should have found that the superstability condition $f'(x) = 0$ gives $x_0 = 1/2$. Argue that $x_0 = 1/2$ must be in the superstable orbits of higher period. I.e., let \begin{gather*} x_1 = f(x_0), \qquad x_2 = f(x_1), \qquad \text{etc.} \end{gather*} Argue that, if you find a period $n$ fixed point such that $x_n = x_0$, then the condition that the orbit is superstable \begin{gather*} \frac{\mathrm{d}{f^{(n)}}(x_0)}{\mathrm{d} x_0} = 0 \end{gather*} implies that the point $1/2$ is on the orbit. We will start here $x_0 = 1/2$. ```{code-cell} # Play with this, or make your own plot, and convince yourself that # the derivative of f^(n)(1/2) is always zero. try: from ipywidgets import interact except ImportError: # On ReadTheDocs, we disable interact def interact(**kw): def wrapper(f): f() return wrapper @interact(r=(1.0, 4.0), n=(1, 20)) def draw(r=1.0, n=1): x = x0 = np.linspace(0, 1, 1000) for _n in range(n): x = f(x, r) plt.plot(x0, x) ``` This gives us a strategy. We simply look for roots of the objective function as a function of $r$. \begin{gather*} f^{(n)}_{r}(\tfrac{1}{2}) - \tfrac{1}{2} = 0. \end{gather*} ```{code-cell} def objective(r, n): """Return f^n(1/2) - 1/2.""" x = 0.5 for _n in range(n): x = r*x*(1-x) return x - 1/2 @interact(n=(1, 20), a=(1.0, 4.0, 0.001), b=(1.0, 4.0, 0.001)) def draw(n=1, a=1.0, b=4.0): fig, ax = plt.subplots(figsize=(20,3)) r = np.linspace(a, b, 1000) ax.plot(r, objective(r, n)) ax.plot(r, 0*r, 'k') ax.set(xlabel="r") ``` ## A Systematic Approach Consider the equation $f^{(2)}{}'(x) = 0$ first to find a point on the orbit, then compute $r$. :::{solution} We are now looking for solutions to the equations \begin{gather*} x = f^{(2)}(x) = f(f(x)), \qquad 0 = f^{(2)}{}'(x) = f'(f(x))f'(x). \end{gather*} The second equation still has the solution $x=1/2$ as before, so this point must be on the orbit. This point iterates to: \begin{gather*} x_0 = \frac{1}{2}, \qquad x_1 = f(x_0) = -\frac{r}{4}, \qquad x_2 = f(x_1) = rx_1(1-x_1) = -\frac{r}{4}, \qquad \end{gather*} This same reasoning holds for any value of $n$. We then have \begin{gather*} x = r(rx(x-1))(x-1), \qquad 0.5 = r^2(0.5)^3 \end{gather*} ::: You saw from the previous assignment that the logistic map appeared to be chaotic for $r=4$: \begin{gather*} x \mapsto 4x(1-x). \end{gather*} [Staisław Ulam][] and [John von Neumann][] suggested using this as a way of generating pseudo-random numbers. Your task is to explore how well this works and try to improve the naïve implementation. ## References * Section **A Case of Complete Chaos** in Chapter 12 of {cite}`Schroeder:1991`. * [Section 5.1: **The Multiple Reduction Copy Machine Metaphor** of Peitgen et al. "Chaos and Fractals: New Frontiers of Science"](https://ebookcentral.proquest.com/lib/wsu/reader.action?docID=3085512&ppg=252) (Online acces through WSU Library.) * [Section 6.4: **Random Number Generator Pitfall** of Peitgen et al. Chaos and Fractals: New Frontiers of Science](https://ebookcentral.proquest.com/lib/wsu/reader.action?docID=3085512&ppg=354) (Online acces through WSU Library.) * Wikipedia: [Logistic Map: Solutions when $r=4$](https://en.wikipedia.org/wiki/Logistic_map#Solution_when_r_=_4). [Staisław Ulam]: [John von Neumann]: [Sierpinski gasket]: [PCG64]: