|
| 1 | +# Bifurcation diagram of two populations of heterogeneous Kuramoto-Sakaguchi oscillators |
| 2 | + |
| 3 | +These are commands for the Python interface of `auto-07p`. |
| 4 | + |
| 5 | +At the end of this file we provide instructions to install the software, althought this might not work depending on your system configuration. See the official documentation. |
| 6 | + |
| 7 | +Notice that $\mu$ in the manuscript corresponds to `p` in these codes. |
| 8 | + |
| 9 | +## Provided files |
| 10 | + |
| 11 | +The software requires at least two files to run, which we provide and should not be modified: |
| 12 | + |
| 13 | +- `oa2.f90` provides the system equations and the corresponding Jacobian. We also modified the options in the file to provide the maxima and minima of limit-cycles. |
| 14 | +- `c.oa2` provides the default constants for `auto-07`. These can be overwritten within the program, thus no need to modify this file. |
| 15 | + |
| 16 | +## Bifurcations along $K$ for fixed $p=0.9$ |
| 17 | + |
| 18 | +### Fixed points: |
| 19 | + |
| 20 | +First, we use Euler method to find a fixed point for $K=7$: |
| 21 | + |
| 22 | +``` |
| 23 | +init=run('oa2',ICP=[1,2,3],IPS=-2,NMX=100000,PAR={'K': 7.0, 'p' : 0.9, 'alpha' : 1.2}) |
| 24 | +``` |
| 25 | +The fixed point corresponds to an assymetric (chimera) state, as $R_a\neqR_b$. |
| 26 | +We can continue this solution increasing $K$: |
| 27 | + |
| 28 | +``` |
| 29 | +ic=init(201) |
| 30 | +asym=run(ic,IPS=1,NMX=10000,ICP=[1,2]) |
| 31 | +``` |
| 32 | +We see that the continuation stops at $K=200$, because we specified so in the `c.oa2` file. |
| 33 | +Also, a Hopf bifurcation (`HB`) has been detected at $K\approx 7.37$. We will investigate this later. |
| 34 | +For the moment, we save it on a new variable: |
| 35 | + |
| 36 | +``` |
| 37 | +hopf = asym('HB1') |
| 38 | +``` |
| 39 | + |
| 40 | +The Python interface of `auto-07p` is just Python. Thus we can plot the previous results using `matplotlib`: |
| 41 | + |
| 42 | +``` |
| 43 | +import matplotlib.pyplot as plt |
| 44 | +plt.plot(asym['K'],asym['Ra'],asym['K'],asym['Rb']) |
| 45 | +plt.ylabel('Ra,Rb') |
| 46 | +plt.xlabel('K') |
| 47 | +plt.show() |
| 48 | +``` |
| 49 | +*(Depending on your version of numpy, this might produce errors due to deprecation of `np.bool`. |
| 50 | +If this happens, a quick workaround is to do `import numpy as np; np.bool=np.bool_;`)* |
| 51 | +We can now run backwards in $K$ by specifying `DS="-"`: |
| 52 | + |
| 53 | +``` |
| 54 | +new_branch = run(ic,IPS=1,NMX=10000,ICP=[1,2],DS="-") |
| 55 | +``` |
| 56 | + |
| 57 | +A bifurcation at $K\approx 6.66$ is detected (and then Auto turns backwards again). |
| 58 | +From this bifurcation, a new solution branch is detected, and auto computes the solution directly. |
| 59 | + |
| 60 | +The two branches can be accessed as: |
| 61 | +``` |
| 62 | +asym = new_branch[0] # The chimera state for the full $K$-range. |
| 63 | +sym = new_branch[1] # This corresponds to the homogeneous solution of the system. |
| 64 | +``` |
| 65 | + |
| 66 | +Let's visualize the results: |
| 67 | + |
| 68 | +``` |
| 69 | +plt.plot(asym['K'],asym['Ra'],'black') |
| 70 | +plt.plot(asym['K'],asym['Rb'],'black') |
| 71 | +plt.plot(sym['K'],sym['Ra'],'black') |
| 72 | +plt.ylabel('Ra,Rb') |
| 73 | +plt.xlabel('K') |
| 74 | +plt.show() |
| 75 | +``` |
| 76 | + |
| 77 | +We see that the homogeneous branch has been continued to negative values of $R$! |
| 78 | +In order to obtain a physically meaningfull solution, and obtain the entire branch we rerun |
| 79 | +the continuation starting from the Kuramoto synchronization transition that Auto has detected: |
| 80 | +The way Auto works, we cannot initialize the new simulation at `sym('LP1')`, we have to |
| 81 | +specify the original solution instead: |
| 82 | + |
| 83 | +``` |
| 84 | +sol = run(new_branch('LP2'),DS="-") |
| 85 | +``` |
| 86 | + |
| 87 | +Notice we had to change direction again with `DS="-"` (since we were going backwards). |
| 88 | + |
| 89 | +Again, auto computes authomatically additional solutions from pitchfork bifurcations (we can avoid this turning off the detection of new branches). |
| 90 | +From these two branches, we are only interested on the first, as we already have the second: |
| 91 | + |
| 92 | +``` |
| 93 | +sym=sol[0] |
| 94 | +plt.plot(asym['K'],asym['Ra'],'black') |
| 95 | +plt.plot(asym['K'],asym['Rb'],'black') |
| 96 | +plt.plot(sym['K'],sym['Ra'],'black') |
| 97 | +plt.ylabel('Ra,Rb') |
| 98 | +plt.xlabel('K') |
| 99 | +plt.show() |
| 100 | +``` |
| 101 | + |
| 102 | +We observe a new pitchfork (`BP`) of the homogeneous state at $K\approx 15.67$. |
| 103 | +Auto did not check this branch. We can continue it with: |
| 104 | + |
| 105 | +``` |
| 106 | +another_branch=run(sol('BP2'),ISW=-1,STOP=['BP1']) |
| 107 | +``` |
| 108 | + |
| 109 | +Notice the `ISW=-1` to force a branc switch, otherwise auto computes the solution we already know. Also we use `STOP=['BP2']` to avoid recomputing solutions already known for us. |
| 110 | + |
| 111 | +This continuation provides two new branches. One is a branch of (unstable) assymetric fixed points that vanish on a subcritical pitchfork bifurcation to the antiphase state, which |
| 112 | +is a new solution for us: |
| 113 | + |
| 114 | +``` |
| 115 | +asym2 = another_branch[0] |
| 116 | +antiphase = another_branch[1]; |
| 117 | +``` |
| 118 | + |
| 119 | +Let's visualize this: |
| 120 | + |
| 121 | +``` |
| 122 | +plt.plot(asym['K'],asym['Ra'],'black') |
| 123 | +plt.plot(asym['K'],asym['Rb'],'black') |
| 124 | +plt.plot(sym['K'],sym['Ra'],'black') |
| 125 | +plt.plot(asym2['K'],asym2['Ra'],'gray') |
| 126 | +plt.plot(asym2['K'],asym2['Rb'],'gray') |
| 127 | +plt.plot(antiphase['K'],antiphase['Ra'],'blue') |
| 128 | +plt.ylabel('Ra,Rb') |
| 129 | +plt.xlabel('K') |
| 130 | +plt.show() |
| 131 | +``` |
| 132 | + |
| 133 | +We should continue the antiphase solution downstream and this would have provided all the relevant fixed points in the system for this value of $p$. For instance, we can do |
| 134 | + |
| 135 | +``` |
| 136 | +antiphase = run(another_branch('UZ1'),DS="-",ISW=1,STOP=[]) |
| 137 | +``` |
| 138 | + |
| 139 | +The plotting commands used previously should show now the full antiphase state. |
| 140 | +Notice that here auto breaks when $R=0$ (as it should). |
| 141 | +Also notice a Hopf bifurcation from the antiphase state. We save it for later use: |
| 142 | + |
| 143 | +``` |
| 144 | +hopf2 = antiphase('HB1') |
| 145 | +``` |
| 146 | + |
| 147 | +### Limit-cycles |
| 148 | + |
| 149 | +Now, let's turn our attention back to the limit-cycle solutions |
| 150 | +emerging from the Hopf of the chimera states. This Hopf bifurcation is at |
| 151 | + |
| 152 | +``` |
| 153 | +hopf['K'] |
| 154 | +``` |
| 155 | + |
| 156 | +Auto can authomatically continue the resulting limit-cycles. |
| 157 | +We just have to specify that we are interested in periodic orbits with `IPS=2`. |
| 158 | +We also turn on the detection of bifurcations from periodic orbits with `ISP=2`: |
| 159 | + |
| 160 | +``` |
| 161 | +lc=run(hopf,IPS=2,ISP=2,ICP=[1,11,3,4,5],NMX=50000,ISW=1,DSMAX=0.01,NTST=200,NCOL=7, STOP=['BP1']) |
| 162 | +``` |
| 163 | + |
| 164 | +Now, the simulation halts at $K\approx 7.88$ without detecting any bifurcation. |
| 165 | +However we can see that the period of the oscillation is quite large, indicating |
| 166 | +a possible homoclinic bifurcation. |
| 167 | + |
| 168 | +A closer inspection shows that the period orbit is colliding with the homogeneous state. |
| 169 | +This can also be seen here if we plot the maxima and minima of the limit cycle, |
| 170 | +together with the fixed points: |
| 171 | + |
| 172 | +``` |
| 173 | +plt.plot(asym['K'],asym['Ra'],'black') |
| 174 | +plt.plot(asym['K'],asym['Rb'],'black') |
| 175 | +plt.plot(sym['K'],sym['Ra'],'black') |
| 176 | +plt.plot(asym2['K'],asym2['Ra'],'gray') |
| 177 | +plt.plot(asym2['K'],asym2['Rb'],'gray') |
| 178 | +plt.plot(antiphase['K'],antiphase['Ra'],'blue') |
| 179 | +plt.plot(lc['K'],lc['MAX Ra'],'r', lc['K'],lc['MIN Ra'],'r',lc['K'],lc['MAX Rb'],'r',lc['K'],lc['MIN Rb'],'r') |
| 180 | +plt.ylabel('Ra,Rb') |
| 181 | +plt.xlabel('K') |
| 182 | +plt.show() |
| 183 | +``` |
| 184 | + |
| 185 | +Let's go now to the Hopf bifurcation of the antiphase state, and see if it can complete the picture we have here: |
| 186 | + |
| 187 | +``` |
| 188 | +lc2=run(hopf2,IPS=2,ISP=2,ICP=[1,11,3,4,5],NMX=50000,ISW=1,DSMAX=0.01,NTST=200,NCOL=7, STOP=['BP1']) |
| 189 | +``` |
| 190 | + |
| 191 | +So, finally, we have that: |
| 192 | + |
| 193 | +``` |
| 194 | +plt.plot(asym['K'],asym['Ra'],'black') |
| 195 | +plt.plot(asym['K'],asym['Rb'],'black') |
| 196 | +plt.plot(sym['K'],sym['Ra'],'black') |
| 197 | +plt.plot(asym2['K'],asym2['Ra'],'gray') |
| 198 | +plt.plot(asym2['K'],asym2['Rb'],'gray') |
| 199 | +plt.plot(antiphase['K'],antiphase['Ra'],'blue') |
| 200 | +plt.plot(lc['K'],lc['MAX Ra'],'r', lc['K'],lc['MIN Ra'],'r',lc['K'],lc['MAX Rb'],'r',lc['K'],lc['MIN Rb'],'r') |
| 201 | +plt.plot(lc2['K'],lc2['MAX Ra'],'r', lc2['K'],lc2['MIN Ra'],'r',lc2['K'],lc2['MAX Rb'],'r',lc2['K'],lc2['MIN Rb'],'r') |
| 202 | +plt.ylabel('Ra,Rb') |
| 203 | +plt.xlabel('K') |
| 204 | +plt.show() |
| 205 | +``` |
| 206 | + |
| 207 | +We see the two limit-cycle branches are about to join, but the continuation from the antiphase solution breaks down when the orbit touches $R=0$. To obtain the trajectory in this narrow space one can perform simulations of the system. In spite of this (numerical) constrain, a detailed analysis shows that the limit-cycles from both Hopf solutions join at a Double Homoclinic bifurcation (see manuscript for more details). |
| 208 | + |
| 209 | +All these solutions can be exported either using Python or directly with the auto commands, e.g., `save(lc,"lc")`. The computed data contains more information, such as the stability of each solution and the corresponding eigenvalues, see the official documentation for more information. |
| 210 | + |
| 211 | +Auto generates several auxiliary files that might not be needed. |
| 212 | +Before closing, do not forget to `clean()`! |
| 213 | + |
| 214 | + |
| 215 | +## Installing `auto-07p` |
| 216 | + |
| 217 | +In order to install `auto-07p` in Linux from the official repository you can use: |
| 218 | + |
| 219 | +``` |
| 220 | +mkdir auto-07p |
| 221 | +git clone https://github.com/auto-07p/auto-07p auto-07p |
| 222 | +cd auto-07p |
| 223 | +./configure |
| 224 | +make |
| 225 | +make install |
| 226 | +``` |
| 227 | + |
| 228 | +Depending on your system, there might be some conflicts. |
| 229 | +I suggest installing a minimal version without the provided plotting tools, |
| 230 | +as they require some dependencies that are outdated or conflict with current packages. |
| 231 | +To do so, diable them in the configuration step of the previous instructions: |
| 232 | + |
| 233 | +``` |
| 234 | +./configure --enable-plaut04=no --enable-plaut=no --enable-plaut-qt=no |
| 235 | +``` |
| 236 | + |
| 237 | +Some other conflicts might appear anyhow, please see the official documentation. |
| 238 | + |
| 239 | +If everything goes according to plan, typing `auto` in a terminal should start the interface. |
| 240 | + |
| 241 | + |
| 242 | + |
| 243 | + |
| 244 | + |
| 245 | + |
0 commit comments