[SOLVED] MTH3039 2024 Computational Nonlinear Dynamics Coursework 2 Matlab

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.