Computational Nonlinear Dynamics
(MTH3039, 2024) Coursework 2
Advanced concepts of stability in dynamical systems
This coursework is worth 60% of the credits for this unit. The maximum number of marks for this coursework is 100, separated into into thematic task groups.
Deadline for this coursework is 13th of December (Friday), 11:59am.
Submitting your work
Submission is done via ELE. You need to submit a .zip file containing:
• A Jupyter notebook CW2 .ipynb which contains the solutions for all tasks. The notebook is run so that all code is accompanied with its numeric or visualized output, and text-based explanations and comments are appended where necessary. You must structure the notebook in sections corresponding to the individual tasks in this coursework (to add a section create a markdown cell and type # Task 1). The notebook must contain the line include("CW2.jl") in its first code block.
• A Julia source code file CW2.jl that contains function definitions that you use in the above Jupyter notebook. The file CW2.jl must contain definitions for all named functions mentioned in this coursework, e.g., a function named fg_attractors or ppha_mapper, etc.. Once CW2.jl is run, all named functions mentioned in this coursework must be accessible.
• (Optional) Any other Julia source code files that you include() in the above file, in case you want to organize your codebase into several files. You are allowed to put all of your Julia source code into the single CW2.jl file if you wish so.
• (Optional) Any binary files that you saved to disk while working on the coursework, see the discussion below on JLD2.
Marking criteria for code
Marking criteria for submitted Jupyter notebook and source code:
• 100%: Code output (typically a figure) is fully correct and visually clear: different things are colored differently or are plotted with different markers or linestyles, etc. Code runs and creates the same figure as saved in the Jupyter notebook. Code runs within 10 minutes or less per Task. Code is readable.
• 75-99%: Code output is mostly correct but it has minor mistakes that can be cor- rected with minor work. Or, code takes unnecessarily long to run even if fully correct.
• 40-74%: Code output is not correct, but there is some code that is related with the given task, and the code runs. The code can be corrected with extensive work.
• 0-39%: There is no output from the code and/or the code does not run.
Point deductions on the quality of the code are the same as with coursework 1 with the exceptions:
• You do not have to submit a PDF version of your Jupyter notebook. But remember that all figure output must be in the submitted Jupyter notebook without re-running it!
• Poorly structured code will not incur any penalties unless it also affects the correct- ness or computational performance of the output. Remember that all individually named functions throughout this coursework must exist by themselves irrespectively of source code structure.
Restrictions on used packages
Throughout this coursework you have to use parts of the DynamicalSystems.jl library. You can create dynamical systems (via the CoupledODEs function). You can use the trajectory function at will, or the interactive trajectory GUI. You can use the query- ing and altering interface for DynamicalSystem such as the functions step! , reinit! , or set parameter! (in general, all functions described in the documentation of the DynamicalSystem type). Other components of DynamicalSystems.jl are forbid- den from your solution code, but you may use them outside your submitted solution to check whether your solutions are correct.
For plotting please use only the Makie package and its derivatives such as CairoMakie. You may use JLD2 to save data to binary files that you can load in later sections of the notebook. This allows you to perform. an extensive computation (such as computing the basins of attraction), then save it to disk, and then resume your work later by loading this file in. You must also include such binary files in your submission if you chose this approach. All packages belonging to the Julia standard library (those installed with Julia, such as LinearAlgebra or Statistics) may be used as well.
No other external packages are necessary for this coursework. If any are used, they need to be justified extensively in the submitted coursework, otherwise the resulting grade will be penalized appropriately. If you intend to submit coursework using such other external packages, we recommend that you ask your lecturer’s approval before doing so.
Multistability [30 points]
1. Consider the following dynamical system:
where A = 2.0551, B = −2.6, C = 0.4, D = 1, E = 0.4. This dynamical system represents complicated predator-prey dynamics between predator population y and prey population x.
The system has three attractors at the given parameters. Find them, not with any formal method, but just by eye, by plotting trajectories of different initial conditions that end up at different attractors. For your initial conditions (x,y) pick random values, with x being between 0.001 to 1, and y between 0.001 to 0.05. Use time sampling of ∆t = 0.1 and evolve initial conditions for a total of T = 200. To show the attractors simply provide a plot with at least three trajectories, each going to each attractor.
2. Estimate and plot the basins of attraction (BoA). To do so, initialize a grid of initial conditions for x ranging from 0.001 to 1, and for y ranging from 0.001 to 0.05. Use a total of 100 values in both ranges, which gives you 100×100 = 10,000 initial conditions. Each initial condition is evolved for time T = 500 and then given a color according to the attractor it converged to. The BoA are a heatmap of the initial conditions grid color-coded by the attractors.
Classify each initial condition (x(0), y(0)) according to the final point (x(T), y(T)) of its trajectory and the following thresholds:
• (x(0), y(0)) converged to attractor 1 if y(T) ≤ ythr and x(T) < xthr ,
• (x(0), y(0)) converged to attractor 2 if y(T) ≤ ythr and x(T) ≥ xthr ,
• (x(0), y(0)) converged to attractor 3 if y(T) > ythr.
I suggest the values xthr = 0.5 and ythr = 10−3 . Convince yourself that the result does not depend strongly on the threshold values.
3. Perform. the same exercise for shorter and shorter T. Find a value of T below which the resulting basin-of-attraction plot start to deviate visibly from the plot for T = 500 (this means that T was not large enough for trajectories to converge to their respective attractors).
4. Create your own general purpose “featurize and group” function for mapping initial conditions into attractors. Name your function fg simple(ds, featurizer, ics; Ttr = 400, Dt = 1, T = 100, r = 0 .1). The input arguments are an instance of a DynamicalSystem (in your case this is a CoupledODEs), a featurizing function (see below), and an iterable of initial conditions ics. Inside the function, each initial condition generates a trajectory according to the keywords: first evolve for transient time Ttr, then store the trajectory X from Ttr to Ttr + T by sampling with time step Dt. Cast each trajectory X into a feature-vector using the featurizer input, which is itself a single-argument function that takes as an input a trajectory (vector of state space points) and outputs a feature (a single vector).
For grouping the features into individual groups, create a function
simple_grouping(features, r) that you call in your main function using the r keyword. simple_grouping implements the simple grouping algorithm (via group centers) shown in the slides of Week 7 with threshold r and returns a vector of the group-labels of the features. fg_simple returns this same vector as the attractor labels.
Apply your function to a dynamical system representing the predator-prey model of (1) using the same initial conditions as in the previous Task. For this and the next exercises use the following featurizer:
using Statistics: mean
function featurizer(X) # X is the trajectory, a vector of points
x = mean(u[1] for u in X) # average of 1st dimension
y = mean(u[2] for u in X)/0 .05 # scaled average of 2nd dimension return [x, y]
end
Create the same BoA plot using the new labels you estimated.
5. So far your function fg simple is finding groups that represent attractors. Extend it into a new function fg attractors such that the new version also returns the actual attractors in state space. To do this, you need to keep track of which initial conditions generated the features that became group “centers” . For each initial condition that became a group center, re-create the trajectory for it (with same Ttr, Dt, T), and store it in memory. These trajectories are the different attractors, and these should be returned as a second output of your (now modified) function. Apply this extended function to the example system of the previous task and confirm that you obtain the same attractors you found in Task 1. Scatterplot these attractors on top of the BoA plot which you have created in the previous Task.
6. Extend your function fg attractors to fg attractors fractions, such that the extended function also returns a third output: the BoA fractions of each attractor. Your function ppha attractors fractions should return three outputs: the labels of the input initial conditions (vector of integers), the attractors (vector of state space sets), and the fractions of the BoAs.
7. Apply your generic function fg_attractors_fractions to the following dynamical system
with parameters F = 6.9, G = 1.3, a = 0.24, b = 4.07. Use a fixed set of initial conditions given by
yg = xg = zg = range(-1, 3; length = 11)
ics = [[x, y, z] for x in xg for y in yg for z in zg]
Plot a scatter-plot of the system attractors with different colors in 3D space. Add a plot label to each attractor plot that is their respective state space basin fractions (or alternatively, show the state space fractions as the title of the axis, or to the right of the axis, or simply print it in your Jupyter notebook output).
Hint: for this dynamical system you need to come up with your own featurizerfunc- tion, you can’t re-use the one from the previous Tasks. Create one by examining trajectories of this dynamical system. You also need to decide a value for the thresh- old r as it is based on the featurizer. You can set the Ttr keyword to 100 for this Task.
Nonlocal Stability [20 points]
8. This task continues from Task 5 (or 6), where you have calculated the BoA of the predator prey model with Eqs. (1) and plotted them along with the attractors. Use the approximated BoA to calculate the minimal critical shock (MCS) for the two fixed point attractors of the system using the technique discussed in the slides of Week 8 (that utilizes the BoA). Include these two vectors as arrows into the plot of the basins from Task 5.
To calculate distance between points in the state space of this dynamical system, you need to use a weighted Euclidean distance. Use the function
9. The dynamical system (1) has a limit cycle attractor (at the quoted default pa- rameters). Calculate the minimal critical shock for 100 points along the limit cycle trajectory sampled with sampling time δt = Π/100 where Π is the approximate period of the limit cycle estimated approximately from a timeseries plot of any of the variables. For the initial condition to generate the limit cycle trajectory, use any point on the limit cycle. Report the magnitude and direction of the average minimal critical shock. Compare this magnitude with the average of the magnitudes of the individual minimal critical shocks for each point on the limit cycle. Explain the difference between the two.
10. Now estimate the MCS without relying on the BoA you have already estimated. First, create a new function ppha_mapper(ds, u) that given the dynamical sys- tem ds representing Eq. (1), and an initial condition u it returns an integer 1, 2, or 3, corresponding to one of the three attractors u converged to. To map individual initial conditions to attractors use the same algorithm as in Task 2. Then, create a new function ppha_mcs_random(ds, u; n = 10_000) that utilizes ppha_mapper(ds, u). Given an initial condition u, this new function returns the minimal critical shock corresponding to u using the random search algorithm dis- cussed in Week 8 lecture slides. Keyword n is how many random points v to try. Generate random points v in the same ranges for x,y in which you have originally es- timated the basins. Apply this function to the fixed point attractors of Eqs. (1) and show the MCSs you obtain, which are approximately the same as with the previous Task for large n.
11. Estimate the final state sensitivity averaged over the whole BoA that you estimated in Task 5 (or 6) following the algorithm outlined in the slides of Week 8. The BoA has a 100×100 size. First use a 10×10-sized box to estimate in it the probabilities pij , and then use a 5×5-sized box. Report the value of the average final state sensitivity in both cases and explain the similarity or difference between the two numbers.
Global continuation [30 points]
12. Consider the tristable predator prey model of Eq. (1) at the quoted parameters. Perform. a global continuation for this system while varying C from 0.3 to 0.6 with steps of 0.01. To do this, need tore-use the function fg attractors fractions(ds, featurizer, ics; kw . . .) with the same featurizer as in Task 6. Use the generic continuation algorithm described in the slides of Week 9 and ignore the “seeding” component. For matching use the simple matching by centroid distance described in the lecture slides. For all parameter values use the same fixed set of initial conditions to sample the state space:
values = 31
xg = range(0 .01, 1 .1; length = values)
yg = range(0 .001, 0 .1; length = values)
ics = [[x, y] for x in xg for y in yg]
Visualize your continuation by plotting three quantities versus the C parameter: 1) the BoA fractions, 2) the mean of the x coordinate of the attractors, 3) the standard deviation of the y coordinate of the attractors (use a different Axis for each of the quantities, don’t plot all of them in the same axis).
Hint: If your code takes too long to run (10+ minutes), you can decrease values to 21.
13. In this range C ∈ [0.3, 0.6] the system undergoes three bifurcations: scenarios where either the type, or number of attractors, change fundamentally. Without perform- ing linear stability analysis, but simply by looking at the plot you generated in the previous task (or equivalently by analyzing the attractor output), report the ap- proximate values Ci, i = 1, 2, 3 of when the bifurcations occurred, as well as what happened during the bifurcation. Hint: you do not need to attach formal names to the bifurcations, such as “Hopf” or “Saddle-node” . Just descriptions “attractor of type X appears/dissapears” .
14. Consider the parameter curve given by: C(θ) = 0.4 + 0.05cos(θ), E(θ) = 0.4 + 0.1sin(θ),θ ∈ [0, 2π]. This defines an ellipsis in the C-E parameter plane. Perform. a global continuation, by re-using the code of task 12, but now versus this parameter curve as θ varies from 0 → 2π (use 100 discrete values in this interval).
15. Since this parameter curve is closed (its end is also its beginning), it is guaranteed that the attractors we have found at θ = 2π are exactly the same as when we started the continuation at θ = 0. However, with the simple matching scheme we are re-using from Task 12, the attractors will not be matched correctly because some attractors that existed when θ = 0 have disappeared as θ > 0, and then reappeared when θ → 2π . Enhance the matching part of your implementation so that it also considers vanished attractors when matching, as is discussed in the lecture slides.
16. In this global continuation over the closed parameter curve, there is one more new local bifurcation that takes place which we did not encounter during the continuation of task 12. Specify what this bifurcation is, without using linear stability analysis, but simply looking at the continuation curves (mean of x coordinate and standard deviation of the y coordinate). Formally label the bifurcation from one of the local bifurcations you have learned so far (e.g., “Saddle-Node”, “Pitchfork”, “Hopf”,etc.).
Tipping points (14 points)
17. A simple model that can showcase rate-dependent tipping is Stommel’s box model for Atlantic thermohaline circulation given by
Here T,S are dimensionless temperature and salinity differences between the polar and equatorial ocean basins respectively, and η,α,β are parameters. Keep α = 1,β = 0.3 and produce a bifurcation diagram for T vs. η from 2 to 4. To produce the bifurcation diagram you need to utilize code that you wrote in coursework 1.
18. Continue from the above task and assume that the system starts with η = 2.5. If we are increasing η infinitesimally slowly, what kind of tipping will the system showcase? Find the critical value ηc at which this type of tipping will happen.
19. Continue from the above task and do a simulation where η is increasing from 2.5 to 3.3 with a finite constant rate δη during the time evolution. Always start your simulation with the system being at the only fixed point at η 1 = 2.5. For each δη find at which η value the system is tipping from the initial stable branch to the second one, by crossing the unstable fixed point branch. There is a critical rate δηc below which the system does not tip at all; report this critical rate.
Hint: to solve this task you need to re-define your RHS equations so that they depend on time. In particular, now η is no longer a parameter of the system, but η0 , its starting value is. η is instead a function of time that needs to be evaluated inside the RHS function.
Chaos (6 points)
20. Consider the following dynamical system
with a = 5.2, b = 0.1. Create a function
max_lyap_exp(ds, u; n=10_000, dt=1.0, d0=1e-6) which calculates the maxi- mum Lyapunov exponent of a dynamical system ds representing Eq. (5) and for the given initial condition u (a 3-element vector), utilizing the algorithm described in the slides of Week 11 lecture (rescaling of a test trajectory). The keyword n quantifies how many times to evolve the system for a time period of dt and then to rescale the test trajectory to the reference trajectory. Keyword d0 is both the initial as well as the rescaled distance of the test trajectory from the reference one.
Apply your max_lyap_exp(ds, u) to the following two initial conditions: u1 = [4, 5, 0], u2 = [-4, -5, 0]. As this dynamical system is bistable, these two initial conditions lead to different attractors. Categorize these attractors into one of [fixed point, limit cycle, chaotic] solely based on the maximum Lyapunov exponent and without looking at plots of the attractors (you must justify your choice!). Then visualize the attractors (by plotting trajectories of the two given initial conditions) and confirm your classification.
21. Apply your max_lyap_exp function also to the system of Eqs. (2), but now instead of the quoted parameters use the parameters F,G,a,b = 6.886, 1.347, 0.255, 4.0]. Use the following three initial conditions:
• u1 = [0.1, 0.1, 0.1]
• u2 = [-0.1, 0.5, 0]
• u3 = [-1.5, 1.2, 1.5]
As with the previous task, classify each initial condition’s attractor to one of [fixed point, limit cycle, chaotic] without actually visualizing the attractors. You can optionally confirm your results with an accompanying plot.