Understanding where does the second (stochastic) attractor of the system come from.

Question

I am currently reading a paper, studying the population dynamics in 3-dimensional Lotka-Volterra model with the following interaction change descriptive system: $$ \begin{equation} \begin{cases} \dot x = ax - 0.06x^2-\Large\frac{xy}{x+10}\normalsize,\\ \dot y = -y + \Large\frac{2xy}{x+10}\normalsize-\Large\frac{0.405\cdot yz}{y+10}\normalsize,\\ \dot z = 0.038z^2 -\Large\frac{z^2}{y+20}\normalsize.\\ \end{cases} \end{equation}\tag{1} $$

Where $a$ is the variable parameter, for different values of which the model's behavior is studied.

In the paper, authors claim, that the respective non-degenerate (biologically meaningfull if $a > 0.9158$) equlibria for $(1)$ is $M_3(\bar x_3, \bar y_3, \bar z_3)$, where $\bar x_3 = \frac{a}{0.12}-5+\frac{\sqrt{(a+0.6)^2-1.52}}{0.12},~\bar y_3 = 6.31579,~\bar z_3 = -40.3 + 80.6\frac{\bar x_3}{\bar x_3 + 10}.$

Now, the problem is that further in the paper, authors consider a case, when $(1)$ posesses coexisting limit-1 cycle (red) and two-band chaotic attractor (blue).

But the case is, I understnad what the red curve is, that is just basically the set of solutions of the system $(1)$, plotted not against the time variable, but against each other $(x(t),y(t),z(t))$. I even managed to get the "red" curve by myself, by solving $(1)$ for $a=1.803$ in python, using scipy.integrate.odeint:

But then, I completely fail to understand, where does the "blue" attractor come from, and there is no any further information in the paper about where it comes from as well, rather the blue attractor is just called "attractor of a Feigenbaum's tree".

If anyone would be so kind to explain me, where the "blue" attractor here comes from, that would be great!

Thank you in advance!

Answer 1

Modern ODE solver packages have events and actions. Events test for some condition, actions change state, parameters or event the ODE function at such points. The scipy solver solve_ivp only implements events and the hard-coded actions to register the event points (always) and to terminate the integration (by a flag).

ODE function and event can be implemented as

def model(t,u,a):
    x,y,z = u
    dx = x*(a - 0.06*x - y/(x+10))
    dy = y*(-1+2*x/(x+10)-0.405*z/(y+10))
    dz = z*z*(0.038-1/(y+20))
    return dx,dy,dz

a=1.803

u0 = [29.9,0.2,15.0]

def event(t,u,a): return u[1]-24
event.direction = -1

T = 40
res = solve_ivp(model,(0,10*T),u0, args=(a,), events=event, atol=1e-6, rtol=1e-9)
res = solve_ivp(model,(0,20*T),res.y_events[0][-1], args=(a,), events=event, atol=1e-6, rtol=1e-9)

print(res.t_events[0])

As we are looking for stable attractors, it makes sense to let the solver run for some time to get close to an attractor. The period is about $T=40$, letting it run for 20 periods may seem excessive, 10 are also sufficient.

Now loop through some random initial points from $[25,30]×[0,10]×[0,20]$ and plot the stabilized trajectories, and the event locations in the x-z-plane. One finds that the events lie on a straight line, thus also draw the path of the $x$ coordinate in a successor map.

fig = plt.figure(figsize=(7,4.5), constrained_layout=True)
gs = fig.add_gridspec(2, 3)
ax = [ fig.add_subplot(gs[k, 2]) for k in [0,1]]
ax3 = fig.add_subplot(gs[:, :2])
for k in range(10):
    u0 = np.random.uniform(size=3)*np.array([5,10,20])+np.array([25,0,0])
    res = solve_ivp(model,(0,10*T),u0, args=(a,), events=event, atol=1e-6, rtol=1e-9)
    res = solve_ivp(model,(0,10*T),res.y_events[0][-1], args=(a,), events=event, atol=1e-6, rtol=1e-9)
    ax3.plot(res.y[0],res.y[2], lw=0.8)
    ax[0].plot(res.y_events[0][:,0],res.y_events[0][:,2],'-+', lw=0.5)
    ax[1].plot(res.y_events[0][:-1,0],res.y_events[0][1:,0],'-+', lw=0.5)
for ak in [ax3,*ax]: ak.grid(); 
ax3.set_xlabel("$x$"); ax3.set_ylabel("$z$")
ax[0].set_xlabel("$x_E$"); ax[0].set_ylabel("$z_E$")
ax[1].set_xlabel("$x_E[k]$"); ax[1].set_ylabel("$x_E[k+1]$")
plt.show()

In the resulting plot the stable loop is well-isolated, the point in the lower left corners, and the alternating connection of the two loops of the chaotic attractor is also visible in the second and last plot.