Double click here and enter the student numbers for all group members. Don't enter any names.

Before you start work on the project, click on this link to read the MATH0011 project instructions.

1. 23009447

Bifurcation Diagrams

This project explores the behaviour of iterative sequences, such as the one defined by $x_0 = 1/2$ and

for $n \geqslant 0$. You will see that as $\gamma$ increases the behaviour of the sequence $(x_n)$ changes in an interesting way. For $\gamma$ less than about 5, the sequence $x_n$ converges: e.g. if $\gamma=4$ then $$x_n \to 0.9729550\cdots \textrm{ as } n \to \infty.$$

When $5.1\leqslant \gamma \leqslant 8$ the sequence doesn't converge to a limit, but instead tends to cycle between two values. For example, if $\gamma=6$ then

As $\gamma$ increases the limiting behaviour becomes a cycle of 4 values, then of 8, then 16, and so on. For gamma around 9.4 stops tending to a cycle at all and becomes completely non-periodic, but there are larger values of $\gamma$ at which the sequence again tends to a cycle.

A bifurcation diagram is a diagram which summarises this behaviour. The horizontal axis shows values of $\gamma$, and the vertical axis shows values taken by the sequence after some large number of terms.

You will generate iterative sequences, plot bifurcation diagrams, and try to estimate the Feigenbaum constant, defined as the limiting ratio of differences between parameter values at which the length of the cycle changes. You will get experience in using matplotlib to produce beautiful plots and in writing code that is as efficient as possible. Approximating the Feigenbaum constant is a bit awkward, and you'll need to think carefully about the best way to do it in Python.

Use the next cell for any imports you need.

Exercise 1 - defining and iterating some functions

Write a Python function h(x, g) whose output is $gx (1- \tanh(x))$.

If $f: \mathbb{R} \to \mathbb{R}$ is a function, we define $f^{\circ n}$ to be the function $f\circ f \circ \cdots \circ f$ obtained by composing $f$ with itself $n$ times. Thus $$\begin{align*}f^{\circ 1}(x) & = f(x)\\ f^{\circ 2}(x) &= f(f(x))\\ f^{\circ 3}(x) &= f(f(f(x))) \end{align*}$$ and so on.

Write a function iterate(f, x0, n) which takes a function $f$, a number $x_0$, and a positive whole number $n$ and returns $f^{\circ n}(x_0)$,

Exercise 2 - finding the period of an iteration

The number $y_0$ is said to be $m$-periodic for the function $f$ if $m$ is the smallest positive whole number such that $y_0 = f^{\circ m}(y_0)$. We'd like to check whether $y_0$ is $m$-periodic for some $m$, but because of rounding errors we can't expect $y_0$ to be exactly equal to $f^{\circ m}(y_0)$. For this reason we'll relax our definition a bit by allowing $f^{\circ m}(y_0)$ to be within $\epsilon$ of $y_0$, where $\epsilon$ is some small number.

Write a function periodicity(y0, f, epsilon=0.0001, max_iter=1000) which takes as input a number $y_0$, a function $f$, a number $\epsilon$, and a positive whole number max_iter, and returns the smallest $m$ such that

if there is such an $m$ which is less than max_iter. If there is no such $m$ less than max_iter it should return $-1$.

I've written the first line of the function definition for you, specifying default values for epsilon and max_iter. If you can't remember how these work, review the weekly notebooks or read this link or ask a chatbot.

Here's one test case for your function. Your function h above takes two inputs, but sometimes we want to fix the value of the second input g and use the one-input function that sends a number x to h(x, g). You can create that function with def as usual, or you can write lambda x: h(x, g) which is sometimes simpler. Read more about Python lambdas at this link. If your periodicity function works correctly, then periodicity(0.3507810470182891, lambda x: h(x, 8.5)) should be 4.

Exercise 3 - plotting a bifurcation diagram

In this exercise you are going to write a function to produce bifurcation diagrams for the function $h(x, \gamma)$. The function will take the following inputs:

• g_min and g_max, the maximum and minimum values of $\gamma$ on the x-axis of your plot.

• n_values, the number of equally spaced values of $\gamma$ to use.

• preiterations, the number $M$ referred to in the description below.

• n_points, the number $L$ referred to in the description below.

I've written the first line of the function for you and suggested some default values for these quantities.

Inside the function, generate n_values equally spaced values of $\gamma$ between g_min and g_max. To make the finished diagram look good you probably need n_values to be a few thousand, but you might want to make n_values smaller to begin with to make the code run quickly.

For each of these values of $\gamma$, compute $y_0 = f^{\circ M}(1/2)$ where $f(x)=h(x, \gamma)$ and $M$ is the input preiterations of your function (experiment to see which values make the plot look best - roughly 500 seems to work well). This allows the iterations to settle down to a cycle, if they're going to. If $y_0$ is $m$-periodic - which you can check using your periodicity function with epsilon = 0.0001 say - plot all the points with coordinates

If $y_0$ isn't $m$-periodic for any $m < 1000$, plot all $L$ points with coordinates

where $L$ is the input n_points to your function. I suggest a value of $L$ of about 500 in your final diagram.

Here are some suggestions to help you create the diagram.

• np.linspace generates equally spaced points between given limits.

• You can plot points using plt.scatter. Given lists xs = [a0, a1, a2,...] of x-coordinates and ys = [b0, b1, b2, ...] of y-coordinates, plt.scatter(xs, ys) plots the points with coordinates (a0, b0), (a1, b1), .... You may need to use plt.scatter(xs, ys, s = e, marker='.') where e is some small number to get the plot looking good.

• Experiment with parameter values to get a balance between making the plot look good and making the code run quickly. Start with smaller values of n_value, preiterations, n_points so that the code runs quickly, and increase them when you're confident everything works correctly.

• Only call plt.scatter once - you should generate all of the points in your plot, then use plt.scatter to plot them. If you call plt.scatter lots of times, your code will take a long time to run.

Now call your generate_bifurcation_diagram function to make a plot.

Make a new version of generate_bifurcation_diagram that works for the sequence defined by $x_{n+1} = \gamma x_n (1-x_n)$ and use it to generate another bifurcation diagram. Good values for the maximum and minimum values of $\gamma$ are 2.5 and 4.

Exercise 4 - the Feigenbaum constant

In this final part you are going to try and estimate the Feigenbaum constant $\delta$. This is difficult to do well, and you will only get a very rough approximation by this method.

What you should see in your bifurcation diagram for $h(x, \gamma)$ is that for values of $\gamma$ up to about 5, the iterations converge to a single value. For $\gamma$ between approximately 5 and 8, the iterations eventually have period 2. For $\gamma$ between roughly 8 and 9, they have period 4. This pattern continues with the periods doubling and the intervals of $\gamma$ values getting geometrically smaller.

Let $a_1$ be the value of $\gamma$ at which period 1 changes to period 2, $a_2$ be the value at which period 2 changes to period 4, $a_3$ be the value at which period 4 changes to period 8, and so on, so that $a_1 \approx 5, a_2 \approx 8, a_3\approx 9$. The Feigenbaum constant is defined to be

Try to calculate as many of the $a_i$ as you can, and use them to approximate $\delta$. You can probably only get reliable values for the $a_i$ up to $i=5$ or 6 because of issues with the precision of Python floats. The simplest way to find the $a_i$ is

• produce a list of equally spaced $\gamma$ values between 5 and 9.5

• for each value of $\gamma$, work out the eventual periodicity of the sequence $x_{n+1} = \gamma x_n (1-\tanh(x_n))$ by iterating 1000 times starting at $x_0=1/2$ with your iterate function and then using your periodicity function. Store the period for the $i$th value of $\gamma$ in position $i$ of another list.

• use a loop to find the $\gamma$ values where the periodicity changes from 1 to 2, 2 to 4, 4 to 8, 8 to 16, and so on - these are your $a_i$.

Now calculate as many of the values $a_i$ as you can for the iteration $x_{n+1} = \gamma x_n (1-x_n)$. You should get different $a_i$ to the ones above, but the limiting value $\delta$ of the ratio of the differences of the $a_i$ is known to be the same!

Submitting your project

Check you have done the following.

0. Included all group members' student numbers at the top of this notebook.

1. Read through every exercise to check you have answered every part.

2. Carefully read and followed all of the MATH0011 project instructions.

3. Checked that all of the code in this notebook works correctly.

One group member should download the completed notebook (in CoCalc, click the File menu next to the green Save button, then click Download) and submit it on the MATH0011 Moodle page in the projects section. Please submit only one file per group.