Mini-Lecture 5-06 - The attractor in a chaotic model.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"

def rectangle(x1, y1, x2, y2, **options):
    return polygon([(x1,y1), (x2,y1), (x2,y2), (x1,y2)], **options)

Learning goal:

• Know that even when a model has chaotic behavior, there is an attractor in the state space, and know what this new kind of attractor is called.

Definition: An attractor is a set of points* $A$ in the state space with the following properties:

• A trajectory that starts within $A$ stays within $A$, and

• If a trajectory is perturbed a small distance away from $A$, it will move back towards $A$.

* The “set of points” here can be a single point, a curve, or even (in higher dimensional state spaces) a surface or larger region of the state space.

The two types of attractors we've encountered so far are

1. A stable equilibrium point (or single-point attractor): equilibrium behavior

2. A limit cycle attractor: oscillatory behavior

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

def hastings_interactive():
    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
    delta_t = 0.1
    t_range = srange(0, 10000, delta_t)
    initial_state = (1, 1, 1)
    attractor = desolve_odeint(vectorfield, initial_state, t_range, state_vars)[int(2000/delta_t):]

    t_range = srange(0, 4000, delta_t)
    ics = {
        (3.2, 0.9,  7): "orange", 
        (3.2, 0.4,  7): "green", 
        (0.2, 0.1,  8): "blue", 
        (1.7, 0.6, 11): "darkgray", 
        (0.5, 0.1, 11): "red", 
        (1.7, 0.5,  8.5): "fuchsia", 
    }

    fig = plotly.graph_objects.FigureWidget()
    fig.layout.showlegend = False
    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.xaxis.range = (0,  3.2)
    fig.layout.scene.yaxis.range = (0,  1.2)
    fig.layout.scene.zaxis.range = (5, 11.0)
    fig.layout.scene.aspectmode = "cube"
    fig.add_scatter3d(x=attractor[:,0], y=attractor[:,1], z=attractor[:,2], 
                      mode="lines", line_color="purple", opacity=0.5)
    attractor = fig.data[-1]
    trajectories = []
    for initial_state, color in ics.items():
        solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
        fig.add_scatter3d(x=solution[:,0], y=solution[:,1], z=solution[:,2], 
                          mode="lines", line_color=color)
        trajectories.append(fig.data[-1])

    choices = srange(1, len(ics) + 1) + ["All", "None"]
    @interact(which=slider(choices, default=1, label="Which trajectory:"), 
              show_attractor=checkbox(False, label="Show attractor"))
    def update(which, show_attractor):
        with fig.batch_update():
            attractor.visible = show_attractor
            for n, trajectory in enumerate(trajectories):
                trajectory.visible = (which == "All" or which == n + 1)

    return fig

hastings_interactive()

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

def rossler_interactive():
    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
    t_range = srange(0, 1000, delta_t)
    initial_state = (5, 5, 1)
    attractor = desolve_odeint(vectorfield, initial_state, t_range, state_vars)[int(100/delta_t):]

    t_range = srange(0, 200, delta_t)
    ics = {
        (-25,   0,  40): "orange", 
        (  0, -25,  30): "green", 
        (-15, -15,  45): "blue", 
        ( 25,  25,   1): "darkgray", 
        (0.1, 0.1,  48): "red", 
        (0.1, 0.1,   0): "fuchsia", 
    }

    fig = plotly.graph_objects.FigureWidget()
    fig.layout.showlegend = False
    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.xaxis.range = (-27, 31)
    fig.layout.scene.yaxis.range = (-28, 26)
    fig.layout.scene.zaxis.range = (  0, 65)
    fig.layout.scene.aspectmode = "cube"
    fig.add_scatter3d(x=attractor[:,0], y=attractor[:,1], z=attractor[:,2], 
                      mode="lines", line_color="purple", opacity=0.5)
    attractor = fig.data[-1]
    trajectories = []
    for initial_state, color in ics.items():
        solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
        fig.add_scatter3d(x=solution[:,0], y=solution[:,1], z=solution[:,2], 
                          mode="lines", line_color=color)
        trajectories.append(fig.data[-1])

    choices = srange(1, len(ics) + 1) + ["All", "None"]
    @interact(which=slider(choices, default=1, label="Which trajectory:"), 
              show_attractor=checkbox(False, label="Show attractor"))
    def update(which, show_attractor):
        with fig.batch_update():
            attractor.visible = show_attractor
            for n, trajectory in enumerate(trajectories):
                trajectory.visible = (which == "All" or which == n + 1)

    return fig

rossler_interactive()

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

def lorenz_interactive():
    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, 500, delta_t)
    initial_state = (0.01, 0.01, 0.01)
    attractor = desolve_odeint(vectorfield, initial_state, t_range, state_vars)[int(100/delta_t):]

    t_range = srange(0, 100, delta_t)
    ics = {
        ( 20, -10,  10): "orange", 
        ( 20,  10, -10): "green", 
        ( 40, -20,  25): "blue", 
        ( 40,  20, -25): "darkgray", 
        (  0,  20,  25): "red", 
        ( 25,   8,   8): "fuchsia", 
    }

    fig = plotly.graph_objects.FigureWidget()
    fig.layout.showlegend = False
    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.xaxis.range = (  0, 53)
    fig.layout.scene.yaxis.range = (-20, 20)
    fig.layout.scene.zaxis.range = (-28, 28)
    fig.layout.scene.aspectmode = "cube"
    fig.add_scatter3d(x=attractor[:,0], y=attractor[:,1], z=attractor[:,2], 
                      mode="lines", line_color="purple", opacity=0.5)
    attractor = fig.data[-1]
    trajectories = []
    for initial_state, color in ics.items():
        solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
        fig.add_scatter3d(x=solution[:,0], y=solution[:,1], z=solution[:,2], 
                          mode="lines", line_color=color)
        trajectories.append(fig.data[-1])

    choices = srange(1, len(ics) + 1) + ["All", "None"]
    @interact(which=slider(choices, default=1, label="Which trajectory:"), 
              show_attractor=checkbox(False, label="Show attractor"))
    def update(which, show_attractor):
        with fig.batch_update():
            attractor.visible = show_attractor
            for n, trajectory in enumerate(trajectories):
                trajectory.visible = (which == "All" or which == n + 1)

    return fig

lorenz_interactive()

The three-species food chain model, revisited

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 = 4000
ymax = 12.5
delta_t = 0.1
t_range = srange(0, tmax, delta_t)
initial_state = (0.5, 0.1, 11) # The red one from the original interactive
solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)

fig = plotly.graph_objects.Figure()
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[:,0], y=solution[:,1], z=solution[:,2], 
                  mode="lines", line_color="red")
fig.show()

Time series of the trajectory from the previous slide:

def hastings_timeseries():
    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 = 1000
    ymax = 12.5
    delta_t = 0.1
    t_range = srange(0, tmax, delta_t)
    initial_state = (0.5, 0.1, 11) # The red one from the previous interactive
    solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
    solution = np.insert(solution, 0, t_range, axis=1)

    @interact(show_transient=checkbox(False, label="Show transient and steady state"))
    def update(show_transient):
        t_crit = 180
        p  = list_plot(solution[:,(0,1)], plotjoined=True, color="darkgreen", legend_label="$G$ (grass)")
        p += list_plot(solution[:,(0,2)], plotjoined=True, color="darkgray",  legend_label="$S$ (sheep)")
        p += list_plot(solution[:,(0,3)], plotjoined=True, color="black",     legend_label="$W$ (wolves)")
        if show_transient:
            p += rectangle(0, 0, t_crit, ymax, color="red", fill=True, alpha=0.2)
            p += rectangle(t_crit, 0, tmax, ymax, color="green", fill=True, alpha=0.2)
            p += text("transient", (0.5*t_crit, 0.95*ymax), color="darkred", fontsize=16)
            p += text("steady state →", (0.8*t_crit + 0.2*tmax, 0.95*ymax), 
                      color="darkgreen", fontsize=16)
        p.show(ymax=ymax, axes_labels=("$t$", "Populations"), aspect_ratio=0.5*tmax/ymax, figsize=7.5)

hastings_timeseries()

Conclusions:

1. Even in the presence of chaotic behavior, there is an attractor in the state space. It is called a strange attractor, or just a chaotic attractor.

2. This means that chaotic behavior, just like stable oscillations or stable equilibrium behavior, is a type of steady state behavior.