Mini-Lecture 2 - Attractors, and limit cycle attractors.ipynb

Kernel: SageMath 9.8

from itertools import product
import numpy as np


def find_zeros(field, *ranges, **options):
    """Numerically approximate all zeros of a vector field within a box

    Each 'range' should have the form (x, xmin, xmax), where 'x' is one of the
    variables appearing in the vector field, and 'xmin' and 'xmax' are the
    bounds on that variables.

    Options:
        'intervals' - How many subintervals to use in each dimension (default
            20). If there are zeros that are very close together, you may need
            to increase this to find them. But be warned that the running time
            is roughly this to the n power, where n is the number of variables.
        'tolerance' - How close to approximate roots, roughly (default 1e-4)
        'maxiter' - Maximum number of iterations of Newton's Method to run for
            any one (potential) solution. This defaults to -2*log(tolerance).
            There usually isn't much harm in increasing this, unless there are
            many false hits.
        'round' - Round the components of the solutions to this many decimal
            places (default None, meaning do not round them at all)
    """

    # Initialization
    intervals = options.get("intervals", 20)
    tolerance = options.get("tolerance", 1e-5)
    maxiter = options.get("maxiter", int(round(-2*log(tolerance))))
    roundto = options.get("round", None)
    n = len(ranges)
    if len(field) != n:
        raise ValueError("Dimension of vector field is {}, but {} ranges "
                         "given".format(len(field), n))

    mins  = [xmin - (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]
    maxes = [xmax + (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]
    deltas = [(xmax - xmin) / intervals for xmin, xmax in zip(mins, maxes)]
    powers = [1 << i for i in range(n)]
    J = jacobian(field, [x for x, xmin, xmax in ranges])
    def dist(v, w):
        return sqrt(sum(((a - b) / d)**2 for a, b, d in zip(v, w, deltas)))

    # Set up the array of positive/negative signs of the vector field
    signs = np.zeros((intervals + 1,) * n, dtype=int)
    sranges = [srange(m, m + d*(intervals + 0.5), d) for m, d in zip(mins, deltas)]
    for vertex, index in zip(product(*sranges), product(range(intervals + 1), repeat=n)):
        v = field(*vertex)
        signs[index] = sum(powers[i] for i in range(n) if v[i] > 0)

    # Now search through that array for potential solutions
    solutions = []
    mask = int((1 << n) - 1)
    for index in product(range(intervals), repeat=n):
        indexpowers = list(zip(index, powers))
        all0s = 0
        all1s = mask
        for k in range(mask + 1):
            newindex = [i + (1 if k & p else 0) for i, p in indexpowers]
            vertex = signs[tuple(newindex)]
            all0s |= vertex
            all1s &= vertex
        if all0s & ~all1s == mask:
            # Now do Newton's method!
            v = vector(m + d*(i + 0.5) for i, m, d in zip(index, mins, deltas))
            for i in range(maxiter):
                previous_v = v
                v = v - J(*v).solve_right(field(*v))
                if dist(v, previous_v) < tolerance:
                    break
            else:
                warn("{} iterations reached without convergence for solution "
                     "{}".format(maxiter, v), RuntimeWarning)
            if solutions and min(dist(v, w) for w in solutions) < 2*tolerance:
                continue
            if not all(xmin <= x <= xmax for x, xmin, xmax in zip(v, mins, maxes)):
                continue
            solutions.append(vector(RDF, v))

    # A convenience: round to some number of decimal places, if requested
    if roundto is not None:
        solutions = [vector(round(x, roundto) for x in v) for v in solutions]
    return solutions


def plot_phase_portrait(field, x_range, y_range, init_conds, **options):
    """Plot a vector field, together with equilibrium points and trajectories

    'x_range' should have the form (x, xmin, xmax), and likewise for 'y_range'.
    'init_conds' should be a dictionary whose keys are the points that are to
        be used as initial conditions for trajectories. For each such initial
        condition, the value in the dictionary should be the range of t values
        over which to plot the trajectory starting at that initial point, such
        as srange(0, 20, 0.01) for t from 0 to 20 with a step size of 0.01. If
        you don't want to plot any trajectories, just use {} for this.

    Options:
        'field_color' - The color for vector field arrows (default "limegreen")
        'nullcline_colors' - A pair of colors for the two nullclines (default
            ("lightskyblue", "lightsalmon") )
        'eqpt_color' - The color for equilibrium points (default "darkmagenta")
        'trajectory_color' - The color for trajectories (default "red")
        Specifying None (or any False value) for any of the above will cause
        it to not be plotted at all. (Note that for the nullclines, this must
        be a tuple (None, None).)
    All other options are as for the various 'plot*' commands, and are passed
    directly on to them.
    """

    # Initialization
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    field_color = options.pop("field_color", "limegreen")
    xnullcline_color, ynullcline_color = options.pop("nullcline_colors",
            ("lightskyblue", "lightsalmon"))
    eqpt_color = options.pop("eqpt_color", "darkmagenta")
    trajectory_color = options.pop("trajectory_color", "red")
    options.setdefault("frame", False)
    options.setdefault("axes", True)
    options.setdefault("aspect_ratio", "auto")
    options.setdefault("xmin", xmin)
    options.setdefault("xmax", xmax)
    options.setdefault("ymin", ymin)
    options.setdefault("ymax", ymax)
    p = Graphics()

    # Plot the vector field itself
    if field_color:
        options["color"] = field_color
        p += plot_vector_field(field, x_range, y_range, **options)

    # Plot the "vertical" nullcline
    if xnullcline_color:
        options["color"] = xnullcline_color
        p += implicit_plot(field[0] == 0, x_range, y_range, linewidth=2, **options)

    # Plot the "horizontal" nullcline
    if ynullcline_color:
        options["color"] = ynullcline_color
        p += implicit_plot(field[1] == 0, x_range, y_range, linewidth=2, **options)

    # Find and plot the equilibrium points
    if eqpt_color:
        options["color"] = eqpt_color
        equilibria = find_zeros(field, x_range, y_range)
        p += points(equilibria, size=50, **options)

    # Plot trajectories, if any initial conditions were given
    options["color"] = trajectory_color
    options["thickness"] = 2
    for init_cond, t_range in init_conds.items():
        solution = desolve_odeint(field, init_cond, t_range, [x, y])
        p += list_plot(solution, plotjoined=True, **options)

    return p

We've seen oscillatory behavior before

Recall the Lotka–Volterra predator-prey model, one of the first models we studied in LS 30A.

Specifically, think back to the first version of this we ever looked at: sharks and tuna.

$$\begin{cases} T' = 0.3T - 0.02TS \\ S' = 0.01TS - 0.1S \end{cases}$$

state_vars = list(var("T, S"))  # T = number of tuna, S = number of sharks

system = (                      # The "system" of differential equations: 
    0.3*T - 0.02*T*S,           #     T' = 0.3*T - 0.02*T*S
    0.01*T*S - 0.1*S,           #     S' = 0.01*T*S - 0.1*S
)
vectorfield(T, S) = system

initial_state = (20, 5)
t_range = srange(0, 200, 0.1)

solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
solution = np.insert(solution, 0, t_range, axis=1)
tuna_plot = list_plot(solution[:,(0,1)], plotjoined=True, ymin=0, 
    color="lightblue", thickness=2, legend_label="$T$ (Tuna)")
sharks_plot = list_plot(solution[:,(0,2)], plotjoined=True, ymin=0, 
    color="gray", thickness=2, legend_label="$S$ (Sharks)")
p = tuna_plot + sharks_plot
p.show(axes_labels=("$t$", "Populations"), ymax=50, figsize=5)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [1], line 12
      9 initial_state = (Integer(20), Integer(5))
     10 t_range = srange(Integer(0), Integer(200), RealNumber('0.1'))
---> 12 solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
     13 solution = np.insert(solution, Integer(0), t_range, axis=Integer(1))
     14 tuna_plot = list_plot(solution[:,(Integer(0),Integer(1))], plotjoined=True, ymin=Integer(0), 
     15     color="lightblue", thickness=Integer(2), legend_label="$T$ (Tuna)")
File /ext/sage/9.8/src/sage/calculus/desolvers.py:1716, in desolve_odeint(des, ics, times, dvars, ivar, compute_jac, args, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg)
   1713     des = [des]
   1715 if ivar is None:
-> 1716     all_vars = set().union(*[de.variables() for de in des])
   1717     ivars = all_vars - set(dvars)
   1719     if len(ivars)==1:
File /ext/sage/9.8/src/sage/calculus/desolvers.py:1716, in <listcomp>(.0)
   1713     des = [des]
   1715 if ivar is None:
-> 1716     all_vars = set().union(*[de.variables() for de in des])
   1717     ivars = all_vars - set(dvars)
   1719     if len(ivars)==1:
File /ext/sage/9.8/src/sage/structure/element.pyx:494, in sage.structure.element.Element.__getattr__()
    492         AttributeError: 'LeftZeroSemigroup_with_category.element_class' object has no attribute 'blah_blah'
    493     """
--> 494     return self.getattr_from_category(name)
    495 
    496 cdef getattr_from_category(self, name):
File /ext/sage/9.8/src/sage/structure/element.pyx:507, in sage.structure.element.Element.getattr_from_category()
    505     else:
    506         cls = P._abstract_element_class
--> 507     return getattr_from_other_class(self, cls, name)
    508 
    509 def __dir__(self):
File /ext/sage/9.8/src/sage/cpython/getattr.pyx:356, in sage.cpython.getattr.getattr_from_other_class()
    354     dummy_error_message.cls = type(self)
    355     dummy_error_message.name = name
--> 356     raise AttributeError(dummy_error_message)
    357 cdef PyObject* attr = instance_getattr(cls, name)
    358 if attr is NULL:
AttributeError: 'FreeModule_ambient_field_with_category.element_class' object has no attribute 'variables'

trajectory_plot = list_plot(solution[:500,1:], plotjoined=True, color="red", thickness=2)
trajectory_plot.show(axes_labels=("$T$ (Tuna)", "$S$ (Sharks)"), xmin=0, ymin=0, xmax=90, ymax=90)

But remember, there isn't just one trajectory! There are many!

ics = [(20,5), (20,10), (15,15), (30,3), (60,3)]
p = Graphics()
for initial_state in ics:
    solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)
    p += list_plot(solution[:800], plotjoined=True, color="red", thickness=2)
p.show(axes_labels=("$T$ (Tuna)", "$S$ (Sharks)"), xmin=0, ymin=0, xmax=90, ymax=90)

Thought experiment: Suppose at some point, we kill off 75% of the sharks.

What will happen?

initial_state = (20, 5)         # Start with 20 tuna and 5 sharks again

t_range = srange(0, 80, 0.1)   # Simulate this model from time t = 0 to t = 120
solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)

new_initial_state = solution[-1] * (1, 0.25) # Kill off 75% of the sharks!

new_t_range = srange(t_range[-1] + 0.1, 200, 0.1)
new_solution = desolve_odeint(vectorfield, new_initial_state, new_t_range, state_vars)

solution = np.insert(solution, 0, t_range, axis=1)
new_solution = np.insert(new_solution, 0, new_t_range, axis=1)
full_solution = np.append(solution, new_solution, axis=0)

tuna_plot = list_plot(full_solution[:,(0,1)], plotjoined=True, ymin=0, 
    color="lightblue", thickness=2, legend_label="Tuna")
sharks_plot = list_plot(full_solution[:,(0,2)], plotjoined=True, ymin=0, 
    color="gray", thickness=2, legend_label="Sharks")
p = tuna_plot + sharks_plot
p.show(axes_labels=("$t$", "Populations"), ymax=90)

state_space_plot = list_plot([(0,0)], size=0, xmin=0, xmax=90, ymin=0, ymax=90, 
                             axes_labels=("Tuna", "Sharks"))
frames = []
for i in srange(400, 1600, 5):
    frame = state_space_plot + list_plot(full_solution[:i,(1,2)], plotjoined=True, 
                                         color="red", thickness=1.5)
    frame += point(full_solution[i,(1,2)], size=40, color="black")
    frame += text("t = {:.1f}".format(full_solution[i,0]), (80, 80), color="black")
    frames.append(frame)
a = animate(frames)
#a.save("Sharks and Tuna 1.gif", delay=5)

Thought experiment, continued: Now suppose we kill off 75% of the sharks at a different time in the cycle.

Will we get the same behavior?

initial_state = (20, 5)         # Start with 20 tuna and 5 sharks again

t_range = srange(0, 97, 0.1)    # Simulate this model from time t = 0 to t = 97
solution = desolve_odeint(vectorfield, initial_state, t_range, state_vars)

new_initial_state = solution[-1] * (1, 0.25) # Kill off 75% of the sharks!

new_t_range = srange(t_range[-1] + 0.1, 200, 0.1)
new_solution = desolve_odeint(vectorfield, new_initial_state, new_t_range, state_vars)

solution = np.insert(solution, 0, t_range, axis=1)
new_solution = np.insert(new_solution, 0, new_t_range, axis=1)
full_solution = np.append(solution, new_solution, axis=0)

tuna_plot = list_plot(full_solution[:,(0,1)], plotjoined=True, ymin=0, 
    color="lightblue", thickness=2, legend_label="Tuna")
sharks_plot = list_plot(full_solution[:,(0,2)], plotjoined=True, ymin=0, 
    color="gray", thickness=2, legend_label="Sharks")
p = tuna_plot + sharks_plot
p.show(axes_labels=("$t$", "Populations"), ymax=90)

state_space_plot = list_plot([(0,0)], size=0, xmin=0, xmax=90, ymin=0, ymax=90, 
                             axes_labels=("Tuna", "Sharks"))
frames = []
for i in srange(500, 1700, 5):
    frame = state_space_plot + list_plot(full_solution[:i,(1,2)], plotjoined=True, 
                                         color="red", thickness=1.5)
    frame += point(full_solution[i,(1,2)], size=40, color="black")
    frame += text("t = {:.1f}".format(full_solution[i,0]), (80, 80), color="black")
    frames.append(frame)
a = animate(frames)
#a.save("Sharks and Tuna 2.gif", delay=5)

The kind of behavior we just saw is extremely unusual in systems that occur in nature.

In most “real world” dynamical systems, if the system is oscillating, and something external to the system perturbs the state away from that oscillatory trajectory, the system will eventually return back to that original oscillatory trajectory.

Sound familiar?

This is the same idea as a stable equilibrium point: the system is in equilibrium, and if some external force perturbs the state away from that equilibrium point, the system will eventually return to that equilibrium point.

So we want to generalize the idea of stable equilibrium point so that it applies to oscillatory behavior, not just equilibrium behavior.

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 simplest type of attractor consists of a single point. An attractor with a single point is better known as a stable equilibrium point!

Definition: An attractor for which the set $A$ is a curve that forms a closed loop is called a limit cycle attractor.

For a limit cycle attractor, the closed loop is itself a trajectory (meaning that the system oscillates around that loop). But if the state is perturbed slightly away from that loop, the resulting trajectory will move back towards the loop.

Another way of saying the same thing is as follows: any trajectory that gets close enough to the loop will move closer and closer to it.

A better predator-prey model

This one is called the Holling–Tanner model. You'll study it in your lab this week.

The assumptions of this model, and how to come up with the differential equations, will be covered in another video.

$$\begin{aligned} N &= \text{prey population} \\ P &= \text{predator population} \end{aligned}$$

Differential equations: $$\begin{cases} N' = 0.4 N \cdot \big( 1 - \frac{N}{3000} \big) - \frac{300 N}{1000 + N} \cdot P \\ P' = 0.03 P \cdot \big( 1 - \frac{150 P}{N} \big) \end{cases}$$

# Parameters:
r1 = 0.4       # Natural per-capita growth rate of the prey population
r2 = 0.03      # Natural per-capita growth rate of the predator population
j =  150       # Size of prey population needed to support each predator
w =  300       # Maximum rate at which one predator consumes prey (per year)
d =  1000      # Prey population at which the predation rate will be half of that maximum
k =  3000      # Maximum prey population supported by the environment (carrying capacity)

state_vars = list(var("N, P"))
system = (
    r1*N*(1 - N/k) - w*P*N/(d + N),
    r2*P*(1 - j*P/N),
)
vectorfield(N, P) = system

t_range = srange(0, 2000, 0.1)
init_conds = (
    ((220, 0.8), "red"), 
    ((500, 3), "salmon"), 
    ((500, 3.5), "orange"), 
    ((1000, 0.5), "gold"), 
    ((300, 1.6), "fuchsia"), 
)
trajectory_plots = []
for ics, color in init_conds:
    solution = desolve_odeint(vectorfield, ics, t_range, state_vars)
    trajectory_plots.append(list_plot(solution, plotjoined=True, color=color))

@interact
def attractor(index=slider(srange(0, len(trajectory_plots) + 1, 1) + ["All"], default=0, label="Show trajectory:"), 
              show_attractor=checkbox(False, label="Show limit cycle attractor")):
    p = plot_vector_field(vectorfield, (N, 0, 2000), (P, 0, 4), color="limegreen")
    if index == "All":
        p += sum(trajectory_plots)
    elif index > 0:
        p += trajectory_plots[index - 1]
    if show_attractor:
        p += list_plot(solution[-1000:], plotjoined=True, color="purple", thickness=2)
    p.show(frame=False, axes_labels=("$N$ (Prey)", "$P$ (Predators)"), xmax=2000, ymax=4)

@interact(steady_state=checkbox(False, label="Show steady state bounds"))
def timeseries1(steady_state):
    ics, color = init_conds[2]
    solution = desolve_odeint(vectorfield, ics, t_range, state_vars)
    pmax = max(solution[-2000:,1])
    pmin = min(solution[-2000:,1])
    p2 = list_plot(list(zip(t_range, solution[:5000,1])), plotjoined=True, color=color)
    if steady_state:
        p2 += plot(float(pmin), (x, 0, 500), color="purple", linestyle="dashed")
        p2 += plot(float(pmax), (x, 0, 500), color="purple", linestyle="dashed")
    p2.show(ymin=0, axes_labels=("$t$", "$P$ (Predators)"), figsize=6)

@interact(steady_state=checkbox(False, label="Show steady state bounds"))
def timeseries2(steady_state):
    ics, color = init_conds[4]
    solution = desolve_odeint(vectorfield, ics, t_range, state_vars)
    pmax = max(solution[-2000:,1])
    pmin = min(solution[-2000:,1])
    p2 = list_plot(list(zip(t_range, solution[:9000,1])), plotjoined=True, color=color)
    if steady_state:
        p2 += plot(float(pmin), (x, 0, 900), color="purple", linestyle="dashed")
        p2 += plot(float(pmax), (x, 0, 900), color="purple", linestyle="dashed")
    p2.show(ymin=0, ymax=3.5, axes_labels=("$t$", "$P$ (Predators)"), figsize=6)

Definition: When the solution of a differential equation approaches a limit cycle attractor (or any attractor), the part of the solution before it “reaches” the attractor is called the transient part of the solution. The rest, after the solution has gotten so close to the attractor that it's practically indistinguishable from the attractor itself, is called the steady state part of the solution, or the asymptotic part.

As the above definition notes, this concept of transient vs steady state parts of a solution applies to any attractor, not just a limit cycle.

So, in fact, we've seen examples of this before... with the one other type of attractor we've studied.

Below are some time series of the Allee effect model: $$X' = 0.02X \cdot \big( 1 - \tfrac{X}{250} \big) \cdot \big( \tfrac{X}{40} - 1 \big)$$

@interact(initial_state=slider([150, 100, 50, 300, 350, 42, 30, 20], label="Initial state"), 
          tmax=slider([100, 200, 400], label="$t_{max}$"))
def allee_timeseries(initial_state, tmax):
    diffeq(X) = 0.02*X*(1 - X/250)*(X/40 - 1)
    df = diff(diffeq, X)
    t_values = srange(0, tmax, 0.01)
    p = point((0, initial_state), size=40, color="black")
    #p += text("Initial state: {}".format(initial_state), (tmax/2, 0.8*ymax))
    solution = desolve_odeint(diffeq, initial_state, t_values, X)
    p += list_plot(list(zip(t_values, solution)), plotjoined=True, color="blue", thickness=2)
    eq_points = [expr.rhs() for expr in solve(diffeq(X) == 0, X)]
    for eq_point in eq_points:
        stability = df(eq_point)
        color = "green" if stability < 0 else "red" if stability > 0 else "orange"
        p += list_plot(((0, eq_point), (tmax, eq_point)), plotjoined=True,
                color=color, thickness=2, linestyle="dashed")
    show(p, xmax=tmax, ymin=0, ymax=350, axes_labels=("$t$", "$X$"))

We've just seen that a solution that approaches a stable equilibrium point also has a transient part and a steady state part: the transient part is the part of the solution before the system reaches equilibrium; the steady state part is what happens after it has reached equilibrium, i.e. when it remains constant.

Often, we refer to the steady state behavior or asymptotic behavior of a solution, meaning what happens to the solution in the long run. We've now seen that steady state behavior can be either constant (a stable equilibrium point) or oscillating (a limit cycle attractor). Can you think of other types of steady state behavior a dynamical system could have?

Aside: You may have seen the word asymptotic or asymptote before, such as vertical or horizontal asymptotes of a graph in pre-calculus or calculus. In all of these contexts, the word refers to what happens as some variable approaches infinity. In the context of dynamical systems, the asymptotic behavior or asymptotic part of a solution is what happens to the behavior/solution as $t$ approaches $\infty$.