Assignment Task
Analysis Of A Coupled Phase-Oscillator Model
Introduction This coursework is based on a variant of the Kuramoto model for a population of N globally coupled limit cycle oscillators that you investigated in Workshop 3. The phase dynamics of the jth oscillator in the variant model are given by
20240706053431AM-2090518050-1502546089.PNG
where ωj = 2π/Tj is the natural frequency of the jth oscillator and K is the coupling strength. Recall that we assume that the ωj s are drawn from a symmetric, unimodal probability density function g(ω).
As part of this coursework, we will be simulating (for finite N), the transition from asynchrony to synchrony that occurs with increasing K in eqns. (1), and comparing this to the corresponding transition in the Kuramoto model. As in Workshop 3, we will use the Kuramoto order parameter
20240706053825AM-951776847-786896288.PNG
in which r(t) and Φ(t) measure the collective amplitude and collective phase of the oscillator population. Recall that
20240706053431AM-2074565898-772040526.PNG
where 〈·〉 denotes average over the oscillator population (i.e. over 1 ≤ k ≤ N).
Question 1
In the case N = 2, system (1) can be written in the form
θ˙1 = ω1 + ε sin(2(θ1 − θ2 )),
θ˙2 = ω2 + ε sin(2(θ2 − θ1 )),
where ε is the coupling strength
By considering the relative phase φ = θ1 − θ2, show that if |ω1 − ω2 | < 2>
Question 2
In this question, you will use MATLAB function var_kuramotomod_sim.m provided below to simulate synchronous and asynchronous states in the oscillator population. Set the following:
N = 100
ωj to be normally distributed with mean frequency ω = 0.2618 (so that the mean period is circadian) and standard deviation σω = 0.025
θj (0) to be normally distributed with mean phase 0 and standard deviation σθ = 1
∆t = 0.0
tMAX = 144 Hint: Use the MATLAB function randn to generate normally distributed numbers.
(i) Solve the model with the above settings for coupling strengths K = 0.005 and K = 0.5. Compute and plot the collective amplitude r(t) as a function of time, for both coupling strengths.
(ii) Compute and plot the deviation from the collective mean field, {φj (t) = θj (t) − Φ(t) : 1 ≤ j ≤ N}, for both coupling strengths.
(iii) Plot the initial phase distribution cos θj (0) , sin θj (0) : 1 ≤ j ≤ N and compare with the following:
(a) The phase distribution for K = 0.005 at t = tMAX.
(b) The phase distribution for K = 0.5 at t = tMAX
Question 3
As your simulations for question 2 should show, for a sufficiently large coupling strength K, the oscillator population splits into two synchronized clusters: one with collective phase θ¯(t) and the other with collective phase θ¯(t) + π.
Assume that there are n oscillators in the cluster with phase θ¯(t). Show that the asymptotic collective amplitude, r(∞) = limt→∞ r(t), is given by the following
20240706053431AM-1950110976-1771365910.PNG
Question 4
In this question, you will simulate the change in the collective amplitude r(∞) that occurs as the coupling strength K increases through a critical value Kc , at which the two synchronised clusters form. Set the following:
N = 1000
ωj to be normally distributed with mean frequency ω = 0.2618 (so that the mean period is circadian) and standard deviation σω = 0.025
j (0) to be normally distributed with mean phase 0 and standard deviation σθ = 1
∆t = 0.01
tMAX = 720
(i) For a set of uniformly spaced coupling strengths in the range 0 ≤ K ≤ 0.15, solve eqns. (1) over [0, tMAX ] for each K value, and then estimate the corresponding asymptotic collective amplitude rK(∞) by averaging r(t) over the interval [tMAX − tAV G, tMAX ] to smooth out fluctuations due to the finite size of the oscillation population, where tAV G should be chosen appropriately.
(ii) Plot your estimates of rK(∞) against K. How does the simulated transition compare with the corresponding transition in the Kuramoto model?
(iii) By considering the system of equations obtained from (1) by setting φj = 2θj, calculate the value of Kc in the continuum limit N → ∞, and compare with your simulated transition
Bifuracation Analysis Of The FitzHugh–Nagumo (FHN) model.
Introduction
The FitzHugh–Nagumo model (FHN) is a prototypical excitable system (e.g., a neuron). It is an example of a relaxation oscillator. The FHN model is given by a pair of non-linear differential equations:
The FHN model is a simplified version of the Hodgkin–Huxley model of a spiking neuron. For a = 0 and b = 0 the FHN model simplifies to the Van der Pol oscillator. You will investigate the dynamics of this model in terms of parameters a and b. The other parameters are fixed at c = 2, I = 0. If implemented correctly the right-hand side of (1) should evaluate to (x˙, ˙y) = (4.40691065671982,−1.93915758875542) for a = 6, b = 1 and (x, y) = (p 2,p 3).
Question 1:
One-Parameter Investigation With Brute-Force Methods
(a) For each value of a and b given below, use a brute-force approach to summarise the dynamics across a range of parameter values. Summarise your observations in terms of the observed dynamics (number and types of coexisting solutions). List the types of bifurcation that you can deduce from the plots.
For a = −0.2, a = 0 and a = 0.2, find stable and unstable equilibria and stable periodic orbits (PO) of (1) while varying b ∈ [0.1, 3].
Please identify all unique fixed points (stable and unstable). To this end you need to solve system (1) for zeros (x˙ = 0, ˙y = 0); e.g., by using fsolve, solve, vpasolve or fzero with different initial guesses, you might also consider using nullclines and unique or uniquetol.
Please identify stable equilibria and stable POs. To this end you need to integrate system (1) in time, e.g., by using ode45, and analysing differences between minimum and maximum values reached by a simulation after any transient behaviour has passed). [submit 3 plots, one for each a]
Repeat the analysis for b = 1.03 and b = 1.33, while varying a ∈ [−0.6, 0.6]. [submit 2 plots, one for each b]
(b) Pick (a, b) parameter values at which 2 stable equilibria coexists with a stable periodic orbit. Make a phase-plane plot of (1) for the selected (a, b) values. The plot should include:
two stable equilibria and one unstable equilibrium,
stable PO,
examples of trajectories that converge to each of the solutions,
nullclines, – marked and labelled basins of attraction of each of the solutions.
Question 2
Two-Parameter Investigation With Brute-Force Methods
Summarise the dynamics over the parameter ranges (a, b) ∈ [−0.6,0.6] × [0.1,3] in terms of:
(a) number of the co-existing equilibria (both stable and unstable)
(b) number of the co-existing stable/ attracting solutions (equilibria and POs)
Question 3
Bifurcation Analysis Using Numerical Continuation
(a) Track branches of equilibria varying b (for fixed a = 0) using COCO or other numerical continuation package. Identify all co-dimension 1 bifurcations in your bifurcation diagram. For any Hopf bifurcations branch off and follow POs and identify their bifurcations. Compare the results of the continuation (in the same plot) with the brute-force computations from Q1.
(b) Track any co-dimension 1 bifurcations of equilibria in the (a, b)-plane and identify the locations of Bogdanov-Takens and Cusp co-dimension 2 bifurcations. Compare these results (in the same plot) with the brute-force computations from Q2.