Mini-Lecture 5-07 - Patterns in chaos.ipynb

Kernel: SageMath 9.3

import numpy as np
import plotly.graph_objects

mydefault = plotly.graph_objects.layout.Template()
mydefault.layout.xaxis.showgrid = False
mydefault.layout.yaxis.showgrid = False
mydefault.layout.xaxis.showline = True
mydefault.layout.yaxis.showline = True
mydefault.layout.yaxis.linewidth = 2
mydefault.layout.xaxis.ticks = "outside"
mydefault.layout.yaxis.ticks = "outside"
mydefault.layout.hovermode = False
mydefault.layout.scene.hovermode = False
mydefault.layout.xaxis.showspikes = False
mydefault.layout.yaxis.showspikes = False
mydefault.layout.scene.xaxis.showspikes = False
mydefault.layout.scene.yaxis.showspikes = False
mydefault.layout.scene.zaxis.showspikes = False
plotly.io.templates["mydefault"] = mydefault
plotly.io.templates.default = "mydefault"

Learning goals:

• Be able to explain why there are patterns in the long-term behavior of chaotic dynamics.

• Be able to explain the relationship between this fact, and the fact that we “can't predict the long-term future” in a system with chaotic behavior.

Example: The three-species food chain model (the Hastings model)

Let $G = \text{density of grass}$, ParseError: KaTeX parse error: Expected 'EOF', got '#' at position 11: S = \text{#̲ of sheep (in t…, ParseError: KaTeX parse error: Expected 'EOF', got '#' at position 11: W = \text{#̲ of wolves}.

Assumptions:

• In the absence of sheep, the grass grows logistically.

• Each sheep eats grass at a rate that increases with amount of grass, up to a maximum/saturation level*.

• The birth rate of the sheep is proportional to how much grass they consume.

• The sheep die natural deaths at a constant per-capita rate.

• Each wolf consumes sheep at a rate that increases with the number of sheep, up to a maximum/saturation level*.

• The birth rate of wolves is proportional to how much food (sheep) they consume.

• The wolves die at a constant per-capita rate.

* Use a (non-sigmoid) saturating function for this. In ecology, this type of function is often called
a “Holling type II” function.

The differential equations for this model:

$$\begin{cases} G' = r G \cdot \big( 1 - \frac{G}{k} \big) - m_1 \frac{G}{h_1 + G} \cdot S \\ S' = b_1 \cdot \big( m_1 \frac{G}{h_1 + G} \cdot S \big) - d_1 S - m_2 \frac{S}{h_2 + S} \cdot W \\ W' = b_2 \cdot \big( m_2 \frac{S}{h_2 + S} \cdot W \big) - d_2 W \end{cases}$$

Lots of parameters!

We will use $r = 1$, $k = 3$, $m_1 = 2.5$, $h_1 = 1$, $b_1 = \frac{2}{3}$, $d_1 = 0.4$, $m_2 = 0.1$, $h_2 = 1$, $b_2 = 0.5$, $d_2 = 0.01$.

The differential equations for this model:

$$\begin{cases} G' = G \cdot \big( 1 - \frac{G}{3} \big) - 2.5 \frac{G}{1 + G} \cdot S \\ S' = \frac{2}{3} \big( 2.5 \frac{G}{1 + G} \cdot S \big) - 0.4 S - 0.1 \frac{S}{1 + S} \cdot W \\ W' = 0.5 \big( 0.5 \frac{S}{1 + S} \cdot W \big) - 0.01 W \end{cases}$$

state_vars = list(var("G, S, W"))
system = (
    G*(1 - G/3) - 2.5*G/(1 + G)*S, 
    2/3*2.5*G/(1 + G)*S - 0.4*S - 0.1*S/(1 + S)*W, 
    0.5*0.1*S/(1 + S)*W - 0.01*W, 
)
vectorfield(G, S, W) = system
tmax = 5000
delta_t = 0.1
t_range = srange(0, tmax, delta_t)
initial_state = (1, 1, 1)
solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
solution = np.insert(solution, 0, t_range, axis=1)[int(2000/delta_t):]

end = int(1000/delta_t)
fig = plotly.graph_objects.Figure()
fig.layout.title = "A trajectory of the three-species food chain model"
fig.layout.scene.xaxis.title.text = "G (grass)"
fig.layout.scene.yaxis.title.text = "S (sheep)"
fig.layout.scene.zaxis.title.text = "W (wolves)"
fig.layout.scene.aspectmode = "cube"
fig.add_scatter3d(x=solution[:end,1], y=solution[:end,2], z=solution[:end,3], 
                  mode="lines", line_color="purple")
fig.show()

fig = plotly.graph_objects.Figure()
fig.layout.title = "Time series of the three-species food chain model"
fig.layout.xaxis.title.text = "t"
fig.layout.yaxis.title.text = "Populations"
fig.layout.xaxis.range = (2000, 3000)
fig.layout.yaxis.range = (0, 11)
fig.layout.dragmode = "pan"
fig.add_scatter(x=solution[:,0], y=solution[:,1], mode="lines", line_color="darkgreen", name="G (grass)")
fig.add_scatter(x=solution[:,0], y=solution[:,2], mode="lines", line_color="darkgray",  name="S (sheep)")
fig.add_scatter(x=solution[:,0], y=solution[:,3], mode="lines", line_color="black",     name="W (wolves)")
fig.show()

Example: The Romeo and Juliet and Tybalt model (a.k.a. the Rössler model)

state_vars = list(var("R, J, T"))
system = (
    -J - T,
    R + 0.1*J,
    0.1 - 18*T + R*T,
)
vectorfield(R, J, T) = system
delta_t = 0.01
tmax = 300
t_range = srange(0, tmax, delta_t)
initial_state = (5, 5, 1)
solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
solution = np.insert(solution, 0, t_range, axis=1)[int(100/delta_t):]

fig = plotly.graph_objects.Figure()
fig.layout.title = "A trajectory of the Romeo and Juliet and Tybalt model"
fig.layout.scene.xaxis.title.text = "R (Romeo)"
fig.layout.scene.yaxis.title.text = "J (Juliet)"
fig.layout.scene.zaxis.title.text = "T (Tybalt)"
fig.layout.scene.aspectmode = "cube"
fig.add_scatter3d(x=solution[:,1], y=solution[:,2], z=solution[:,3], 
                  mode="lines", line_color="purple")
fig.show()

fig = plotly.graph_objects.Figure()
fig.layout.title = "Time series of the Romeo and Juliet and Tybalt model"
fig.layout.xaxis.title.text = "t"
fig.layout.yaxis.title.text = "Love/Hate"
fig.layout.xaxis.range = (100, 130)
fig.layout.dragmode = "pan"
fig.add_scatter(x=solution[:,0], y=solution[:,1], mode="lines", line_color="blue",  name="R (Romeo's love)")
fig.add_scatter(x=solution[:,0], y=solution[:,2], mode="lines", line_color="red",   name="J (Juliet's love)")
fig.add_scatter(x=solution[:,0], y=solution[:,3], mode="lines", line_color="black", name="T (Tybalt's hatred)")
fig.show()

Example: The Romeo, Juliet and Juliet's nurse model (a.k.a. the Lorenz model)

state_vars = list(var("R, J, N"))
system = (
    J*N - 8/3*R,
    10*(N - J),
    28*J - N - R*J,
)
vectorfield(R, J, N) = system
delta_t = 0.01
t_range = srange(0, 200, delta_t)
initial_state = (0.01, 0.01, 0.01)
solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
solution = np.insert(solution, 0, t_range, axis=1)[int(100/delta_t):]

fig = plotly.graph_objects.Figure()
fig.layout.title = "A trajectory of the model of Romeo, Juliet, and Juliet's Nurse"
fig.layout.scene.xaxis.title.text = "R (Romeo)"
fig.layout.scene.yaxis.title.text = "J (Juliet)"
fig.layout.scene.zaxis.title.text = "N (Nurse)"
fig.layout.scene.aspectmode = "cube"
fig.add_scatter3d(x=solution[:,1], y=solution[:,2], z=solution[:,3], 
                  mode="lines", line_color="purple")
fig.show()

fig = plotly.graph_objects.Figure()
fig.layout.title = "Time series of the model of Romeo, Juliet, and Juliet's Nurse"
fig.layout.xaxis.title.text = "t"
fig.layout.yaxis.title.text = "Happiness"
fig.layout.xaxis.range = (100, 115)
fig.layout.dragmode = "pan"
fig.add_scatter(x=solution[:,0], y=solution[:,1], mode="lines", line_color="blue", name="R (Romeo)")
fig.add_scatter(x=solution[:,0], y=solution[:,2], mode="lines", line_color="red",  name="J (Juliet)")
fig.add_scatter(x=solution[:,0], y=solution[:,3], mode="lines", line_color="gold", name="N (The nurse)")
fig.show()

Conclusions:

• The existence of a strange attractor (or chaotic attractor) in the state space means that, despite the fact that the behavior is chaotic, we can still find patterns or trends in the behavior, even far into the future.

• This does not contradict what we said a few days ago about sensitive dependence on initial conditions: that property implies that we can't make accurate predictions of the exact state of the system at a specific time far in the future. But there is an important difference between predicting

  1. the exact values of the state variables at specific times, and

  2. patterns of behavior.