<jupyter_start><jupyter_text>Cobwebbing<NAME>Created 4/20/21Based on a Mathematica class demonstration.See a standalone version of the main cobwebbing widget (shown [below](Cobweb-Diagram-for-Logistic-Growth)) [here](./cobwebbing-standalone.ipynb).[Main Project Page](.././index.ipynb)<jupyter_code># Initialization code

%matplotlib widget
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np

# Define parameters
COBWEB_MAX = 20 # number of cobweb iterations to generate

# Define functions

def lmap(x, lim=1.0, r=1.0, h=0.0):
    """Discrete logistic map with absolute harvesting.
    
    Positional arguments:
    x - input value
    
    Keyword arguments:
    lim (1.0) - population limit
    r (1.0) - intrinsic growth rate
    """
    
    return x + r*x*(1 - (x/lim)) - h

def tmap(x, lb=0.5, ub=1.0, r=1.0):
    """Discrete population growth with upper and lower stable population bounds.
    
    Positional arguments:
    x - input value
    
    Keyword arguments:
    lb (0.5) - lower population limit
    ub (1.0) - upper population limit
    r (1.0) - intrinsic growth rate
    """
    
    return x + r*x*(x-lb)*(ub-x)

def mmap(x, r=0.0):
    """Discrete map meant to have an adjustable slope at one of its equlibria.
    
    Positional arguments:
    x - input value
    
    Keyword arguments:
    r (0.0) - slope at intermediate equilibrium point (<= 1.5)
    """
    
    a = 4 - 4*r
    return (a*x*x - 1.5*a*x + (1.0 + 0.5*a))*x

def cobweb_update(x0, cwx, cwy, lim=1.0, r=1.0, mode=0):
    """Updates the global cobweb lists. Lists are edited in-place.
    
    Positional arguments:
    x0 - initial population value
    cwx - reference to a list of cobweb plot x-values
    cwy - reference to a list of cobweb plot y-values
    
    Keyword arguments:
    lim (1.0) - population limit
    r (1.0) - intrinsic growth rate
    mode (0) - 0 for logistic with harvesting, 1 for bounded population, 2 for slope map
    """
    
    # Generate cobweb coordinates
    cwx[0] = x0
    cwy[0] = 0.0
    for i in range(0, 2*COBWEB_MAX, 2):
        cwx[i+1] = cwx[i]
        if mode == 0:
            cwy[i+1] = max(lmap(cwx[i], r=r), 0.0)
        elif mode == 1:
            cwy[i+1] = max(tmap(cwx[i], r=r), 0.0)
        elif mode == 2:
            cwy[i+1] = max(mmap(cwx[i], r=r), 0.0)
        cwx[i+2] = cwy[i+1]
        cwy[i+2] = cwy[i+1]
    cwx[-1] = cwx[-2]
    if mode == 0:
        cwy[-1] = max(lmap(cwx[-1], r=r), 0.0)
    elif mode == 1:
        cwy[-1] = max(tmap(cwx[-1], r=r), 0.0)
    elif mode == 2:
        cwy[-1] = max(mmap(cwx[-1], r=r), 0.0)

# Generate x- and n-values
x = np.linspace(0, 1.5, 101)
nval = [np.floor((n+1)/2) for n in range(2*COBWEB_MAX+2)]<jupyter_output><empty_output><jupyter_text>Difference EquationsDifference equations are the discrete analogs of differential equations, and are used to describe discrete dynamical systems (systems that change over the course of discrete time steps). They are often used in mathematical biology for modeling populations that reproduce in discrete annual cycles rather than continuously year-round.A typical continuous population model would take the form of an autonomous differential equation\begin{align*}\dot{x}(t) = f(x(t))\end{align*}which gives the instantaneous rate of change in the population $x$ at time $t$, $\dot{x}(t) := \frac{d}{dt} x(t)$, as some function $f$ of the current population, $x(t)$.The discrete equivalent takes the form of a difference equation\begin{align*}x_{n+1} = f(x_n)\end{align*}which gives the population $x$ during the next ($n+1$) time period, $x_{n+1}$, as some function $f$ of the population during the current ($n$) time period, $x_n$. The total change in population between time steps $n$ and $n+1$ is denoted $\Delta x_n$, and is defined as\begin{align*}\Delta x_n := x_{n+1} - x_n\end{align*}In this way $\Delta x_n$ is like the discrete version of $\dot{x}(t)$, and the difference can be used for some of the same simple qualitative descriptions as the derivative can. For example, if $\Delta x_n$ is positive then population is gained between time periods $n$ and $n+1$, if $\Delta x_n$ is negative then population is lost, and if $\Delta x_n$ is zero then the population remains constant.  Equilibrium AnalysisOne of the most important questions to be answered about a dynamical systems model is what the expected long-term behavior is like. In the case of a population model we might want to know whether the population will crash to zero (i.e. extinction), or explode without bound, or settle at a finite value. This behavior might in turn depend on some of the parameters of the model, in which case we might want to know how different combinations of parameters affect the long-term behavior of the system. These kinds of _qualitative_ rather than _quantitative_ questions are one of the big differences between a dynamical systems course and a differential equations course. While differential equations courses tend to focus almost entirely on analytically or numerically solving systems in order to be able to write the solution $x(t)$ in closed form, dynamical systems courses tend to focus on questions that can be answered without actually needing to explicitly solve the system.One of the most useful tools in studying dynamical systems is _equilibrium analysis_. An _equilibrium_ of a dynamical system is a solution that does not change once reached, meaning that if the system were to begin at an equilibrium solution, it would remain there forever. It is possible for an equilibrium to be _stable_, meaning that nearby solutions tend to move towards it, or _unstable_, meaning that nearby solutions tend to move away from it. Other, more complicated behaviors (like saddle points, limit cycles, or periodic orbits) are also possible, but one of the main things we usually want to know about a dynamical system is what the equilibria are, and whether they are stable or unstable.Since the defining characteristic of an equilibrium, $x^*$, is that the system doesn't change once it is reached, we can define an equilibrium solution to be any solution for which the difference (or the derivative, in the continuous case) is zero, i.e.\begin{align*}\Delta x^* = 0\end{align*}or equivalently\begin{align*}f(x^*) = x^*\end{align*}So, given a specific model $x_{n+1} = f(x_n)$, we can search for equilibria by simply solving the equation $f(x^*) = x^*$.  Example: Discrete Logistic GrowthOne of the most well-known discrete population models is the discrete logistic map, which is the discrete analog of the continuous logistic growth model. Both provide a simple model for a population whose growth is limited by some carrying capacity (for example, due to limited space or resources).The continuous logistic growth model takes the form of the differential equation\begin{align*}\dot{x}(t) = r x(t) \left( 1 - \frac{x(t)}{K} \right)\end{align*}where $r$ is the intrinsic growth rate of the population and $K$ is the carrying capacity.To analyze this model, we can immediately see that there are exactly two equilibria, $x^* = 0$ and $x^* = K$, since both result in the growth rate $\dot{x}(t)$ becoming zero. When the population is small, $1 - \frac{x(t)}{K}$ is nearly equal to just $1$, and the equation is nearly equal to $\dot{x}(t) = r x(t)$, which is the differential equation for exponential growth. When the population is just under $K$, $1 - \frac{x(t)}{K}$ is nearly zero, and the rate of growth slows. When the population is larger than $K$, then $1 - \frac{x(t)}{K}$ is negative and $\dot{x}(t)$ becomes negative, meaning that the population shrinks. Combining all of these observations we can conclude that $0$ is an unstable equilibrium while $K$ is a stable equilibrium.<jupyter_code>def logistic(t, x0=0.25, r=1.0, lim=1.0):
    """Logistic growth function.
    
    Positional arguments:
    t - input value
    
    Keyword arguments:
    x0 (0.25) - initial population
    r (1.0) - intrinsic growth rate
    lim (1.0) - carrying capacity
    """
    
    if (x0 <= 0.0):
        return t*0.0
    else:
        return lim/(1 + ((lim-x0)/x0)*np.exp(-r*t))

# Set up plot and initialize solution coordinates
fig0, ax0 = plt.subplots()
lx0 = np.linspace(0, 10, 101)
ly0 = np.zeros_like(lx0) # logistic y-coordinates

# Draw plot lines
@widgets.interact(x0=(0.0, 1.5, 0.01))
def update0(x0=0.25):
    global ax0, ly0
    
    # Scatter plot
    ax0.clear()
    ax0.set_xlim([0, 10])
    ax0.set_ylim([0, 1.5])
    ax0.grid(False)
    ax0.set_title("Continuous Logistic Growth")
    ax0.set_xlabel("$t$")
    ax0.set_ylabel("$x(t)$")
    ax0.plot(lx0, logistic(lx0, x0), color="C0")
    ax0.plot(lx0, np.ones_like(lx0), color="black", linestyle="dashed")<jupyter_output><empty_output><jupyter_text>The above widget shows the solution curve for the continuous logistic growth model. The dashed black line indicates the carrying capacity $K$. Click and drag the slider to adjust the initial population $x_0$.The discrete analog of this model is defined by the difference equation\begin{align*}\Delta x_n = r x_n \left( 1 - \frac{x_n}{K} \right)\end{align*}Again we see that the two equilibria are $x^* = 0$ and $x^* = K$, since both result in $\Delta x_n$ becoming zero and thus no change. We can still make some of the same basic observations as in the continuous case about when the population grows and when it shrinks, but the discrete model can display significantly more complicated behavior than the continuous model is capable of, and determining the long-term behavior of any given initial solution is not as straightforward as one might expect.For example, consider the case of $r=2.5$ and $K=1$ (note that this does not necessarily mean that the population is limited to $1$ individual, since the unit of $x$ is unspecified). Then if our initial population were $x_0 = 0.6$, our next population would be $x_1 = 2.5 x_0 (1-x_0) = 2.5 (0.6)(1-0.6) = 1.2$. After that our population would be $x_2 = 2.5 x_1 (1-x_1) = 2.5 (1.2) (1-1.2) = 0.6$, which is equal to the initial population $x_0$, meaning that from this point onward the population will continue to oscillate between $0.6$ and $1.2$ forever. Complicated cycles like this are common in discrete dynamical systems, which can make them difficult to analyze.<jupyter_code># Set up plot and initialize cobweb coordinate lists
fig1, ax1 = plt.subplots()
cwx1 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy1 = np.zeros_like(cwx1) # cobweb y-coordinates
cobweb_update(0.6, cwx1, cwy1, r=2.5, mode=0)
ax1.set_ylim([0, 1.5])
ax1.grid(False)
ax1.set_title("Time Series ($r=2.5$)")
ax1.set_xlabel("$n$")
ax1.set_ylabel("$x_n$")
ax1.plot(np.append([0], nval[1:12+2:2]), np.append([0.6], cwy1[1:12+2:2]), color="C0", marker=".", markersize=10);<jupyter_output><empty_output><jupyter_text>Finding equilibria is relatively straightforward, but determining whether an equilibrium is stable or unstable is another matter. Intuitively we could try to get a rough idea of the stability by plugging a value "close to" $x^*$ into the difference equation and then evaluating it for a few iterations to see whether the solutions get closer to or farther from $x^*$, but that might require going through some lengthy calculations. Fortunately there's a quick and easy graphical method that lets us approximate these calculations visually.  CobwebbingCobwebbing is a graphical method for quickly analyzing the behavior of a discrete dynamical system given a starting value. It is primarily meant to be used with a hand-sketched plot of the difference equation $x_{n+1} = f(x_n)$, and allows the user to come up with a rough idea of how a given initial solution will evolve after many time steps using only a paper and pencil.Given a model of the form $x_{n+1} = f(x_n)$, consider what the graph $y = f(x)$ represents. For any given input value $x$, the height of the point $(x,y)$ on the graph directly above it gives us the value of $f(x)$. Then if we draw the graph of $y = f(x)$, we will be able to quickly, visually estimate the next value $x_1$ which follows our initial value $x_0$. All we need to do is locate $x_0$ on the $x$-axis and then look at the height of the curve directly above $x_0$.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs2, ax2 = plt.subplots(1, 2, figsize=(10, 4))
cwx2 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy2 = np.zeros_like(cwx2) # cobweb y-coordinates

# Generate cobweb plot points
cobweb_update(0.25, cwx2, cwy2, r=1.5, mode=0)

# Draw first part of cobweb plot
ax2[0].clear()
ax2[0].set_xlim([0, 1.5])
ax2[0].set_ylim([0, 1.5])
ax2[0].grid(False)
ax2[0].set_title("Cobweb Plot")
ax2[0].set_xlabel("$x_n$")
ax2[0].set_ylabel("$x_{n+1}$")
ax2[0].plot(x, lmap(x, r=1.5), color="C0")
ax2[0].plot(cwx2[:0+2], cwy2[:0+2], color="C1")
ax2[0].plot(cwx2[0:0+2], cwy2[0:0+2], color="red")

# Draw first part of scatter plot
ax2[1].clear()
ax2[1].set_ylim([0, 1.5])
ax2[1].grid(False)
ax2[1].set_title("Time Series")
ax2[1].set_xlabel("$n$")
ax2[1].set_ylabel("$x_n$")
ax2[1].plot(np.append([0], nval[1:0+2:2]), np.append([0.25], cwy2[1:0+2:2]), color="C0", marker=".", markersize=10)
ax2[1].plot(nval[0+1], cwy2[0+1], color="red", marker=".", markersize=10);<jupyter_output><empty_output><jupyter_text>The pair of plots above show the process of estimating the value of $x_1$ as the height of the curve $f(x_n)$ above $x_0$. The cobweb plot on the left shows this height as a red line, while the time series on the right shows a plot of the population between times $n=0$ and $n=1$.To then find the next value $x_2$ which follows $x_1$ we could repeat the same process, locating $x_1$ on the $x$-axis and finding the height of the curve directly above it, but there's a faster way.After having found the point $(x_0,f(x_0))$ on the plot, the current $y$-coordinate is $f(x_0) = x_1$. We wish to transfer this $y$-coordinate into an $x$-coordinate, and the quickest way to accomplish that visually is to move horizontally to the $y=x$ line, to the point $(x_1,x_1)$.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs3, ax3 = plt.subplots(1, 2, figsize=(10, 4))

# Draw first part of cobweb plot
ax3[0].clear()
ax3[0].set_xlim([0, 1.5])
ax3[0].set_ylim([0, 1.5])
ax3[0].grid(False)
ax3[0].set_title("Cobweb Plot")
ax3[0].set_xlabel("$x_n$")
ax3[0].set_ylabel("$x_{n+1}$")
ax3[0].plot(x, lmap(x, r=1.5), color="C0")
ax3[0].plot(x, x, color="black")
ax3[0].plot(cwx2[:1+2], cwy2[:1+2], color="C1")
ax3[0].plot(cwx2[1:1+2], cwy2[1:1+2], color="red")

# Draw first part of scatter plot
ax3[1].clear()
ax3[1].set_ylim([0, 1.5])
ax3[1].grid(False)
ax3[1].set_title("Time Series")
ax3[1].set_xlabel("$n$")
ax3[1].set_ylabel("$x_n$")
ax3[1].plot(np.append([0], nval[1:1+2:2]), np.append([0.25], cwy2[1:1+2:2]), color="C0", marker=".", markersize=10)
ax3[1].plot(nval[1+1], cwy2[1+1], color="red", marker=".", markersize=10);<jupyter_output><empty_output><jupyter_text>The red horizontal line in the cobweb diagram above runs from $(x_0,x_1)$ to $(x_1,x_1)$, where it intersects the black $y=x$ line. Now our current $x$-coordinate is $x_1$, and we can again move vertically to intersect the curve at the point $(x_1,f(x_1))$ to find $x_2 = f(x_1)$.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs4, ax4 = plt.subplots(1, 2, figsize=(10, 4))

# Draw first part of cobweb plot
ax4[0].clear()
ax4[0].set_xlim([0, 1.5])
ax4[0].set_ylim([0, 1.5])
ax4[0].grid(False)
ax4[0].set_title("Cobweb Plot")
ax4[0].set_xlabel("$x_n$")
ax4[0].set_ylabel("$x_{n+1}$")
ax4[0].plot(x, lmap(x, r=1.5), color="C0")
ax4[0].plot(x, x, color="black")
ax4[0].plot(cwx2[:2+2], cwy2[:2+2], color="C1")
ax4[0].plot(cwx2[2:2+2], cwy2[2:2+2], color="red")

# Draw first part of scatter plot
ax4[1].clear()
ax4[1].set_ylim([0, 1.5])
ax4[1].grid(False)
ax4[1].set_title("Time Series")
ax4[1].set_xlabel("$n$")
ax4[1].set_ylabel("$x_n$")
ax4[1].plot(np.append([0], nval[1:2+2:2]), np.append([0.25], cwy2[1:2+2:2]), color="C0", marker=".", markersize=10)
ax4[1].plot(nval[2+1], cwy2[2+1], color="red", marker=".", markersize=10);<jupyter_output><empty_output><jupyter_text>The pair of plots above show the cobweb diagram and time series after estimating $x_2$ from $x_1$. The new vertical red line runs from $(x_1,x_1)$ to $(x_1,x_2)$. The time series includes the single new population value for $n=2$.This process can be repeated as many times as needed, alternating moving vertically (which may require moving either up or down) to intersect the curve and then moving horizontally (which may require moving either left or right) to intersect the $y=x$ line. The sequence of $x$-coordinates corresponds to the sequence of values $x_0,x_1,x_2,\dots$ attained by the dynamical system.The $y=x$ line also has some special significance beyond simply allowing the current $y$-coordinate to be converted into an $x$-coordinate. On a plot of $x_{n+1}$ versus $x_n$, the $y=x$ line indicates all of the points for which $x_{n+1} = x_n$. Any such points which intersect the curve $y = f(x)$ are solutions for which $f(x_n) = x_n$, which is exactly the definition of an equilibrium, therefore the points for which the curve intersects the $y=x$ line are the equilibrium solutions.  Cobweb Diagram for Logistic GrowthAs discussed earlier, the discrete logistic growth model can display unusual behavior for certain combinations of parameters, and cobwebbing can help to explain why. In particular, let's explore the case for which we fix the carrying capacity $K = 1$ and the initial population $x_0 = 0.25$ and change only the growth rate $r$.For relatively small growth rates the model behaves more or less like the continuous version, with the population increasing from its initial value and stabilizing at the carrying capacity over time.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs5, ax5 = plt.subplots(1, 2, figsize=(10, 4))
cwx5 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy5 = np.zeros_like(cwx5) # cobweb y-coordinates

# Draw plot lines
@widgets.interact(r=(0.25, 1.25, 0.05))
def update5(r=0.5):
    global ax5, cwx5, cwy5
    
    # Cobweb plot
    ax5[0].clear()
    ax5[0].set_xlim([0, 1.5])
    ax5[0].set_ylim([0, 1.5])
    ax5[0].grid(False)
    ax5[0].set_title("Cobweb Plot")
    ax5[0].set_xlabel("$x_n$")
    ax5[0].set_ylabel("$x_{n+1}$")
    cobweb_update(0.25, cwx5, cwy5, r=r, mode=0)
    ax5[0].plot(x, lmap(x, r=r), color="C0")
    ax5[0].plot(x, x, color="black")
    ax5[0].plot(cwx5[:2*COBWEB_MAX+1], cwy5[:2*COBWEB_MAX+1], color="C1")
    
    # Scatter plot
    ax5[1].clear()
    ax5[1].set_ylim([0, 1.5])
    ax5[1].grid(False)
    ax5[1].set_title("Time Series")
    ax5[1].set_xlabel("$n$")
    ax5[1].set_ylabel("$x_n$")
    ax5[1].plot(np.append([0], nval[1:2*COBWEB_MAX+1:2]), np.append([0.25], cwy5[1:2*COBWEB_MAX+1:2]), color="C0", marker=".", markersize=10)<jupyter_output><empty_output><jupyter_text>Click and drag the slider above to change the growth rate $r$. As expected, a larger growth rate causes the population to increase from its initial value to the carrying capacity sooner and in fewer time steps.What might not be expected is that larger growth rates can actually cause the population to "overshoot" the carrying capacity, after which (according to the logistic growth equation) the population must then decrease, causing the population to oscillate around the carrying capacity.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs6, ax6 = plt.subplots(1, 2, figsize=(10, 4))
cwx6 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy6 = np.zeros_like(cwx6) # cobweb y-coordinates

# Draw plot lines
@widgets.interact(r=(1.25, 2.25, 0.05))
def update6(r=1.5):
    global ax6, cwx6, cwy6
    
    # Cobweb plot
    ax6[0].clear()
    ax6[0].set_xlim([0, 1.5])
    ax6[0].set_ylim([0, 1.5])
    ax6[0].grid(False)
    ax6[0].set_title("Cobweb Plot")
    ax6[0].set_xlabel("$x_n$")
    ax6[0].set_ylabel("$x_{n+1}$")
    cobweb_update(0.25, cwx6, cwy6, r=r, mode=0)
    ax6[0].plot(x, lmap(x, r=r), color="C0")
    ax6[0].plot(x, x, color="black")
    ax6[0].plot(cwx6[:2*COBWEB_MAX+1], cwy6[:2*COBWEB_MAX+1], color="C1")
    
    # Scatter plot
    ax6[1].clear()
    ax6[1].set_ylim([0, 1.5])
    ax6[1].grid(False)
    ax6[1].set_title("Time Series")
    ax6[1].set_xlabel("$n$")
    ax6[1].set_ylabel("$x_n$")
    ax6[1].plot(np.append([0], nval[1:2*COBWEB_MAX+1:2]), np.append([0.25], cwy6[1:2*COBWEB_MAX+1:2]), color="C0", marker=".", markersize=10)<jupyter_output><empty_output><jupyter_text>When $r$ is still relatively low these oscillations dampen out and eventually give way to a stable population, but increasing $r$ eventually leads to periodic behavior similar to the alternating solution of $0.6, 1.2, 0.6, 1.2, \dots$ discussed above.Increasing $r$ by still more leads to even more complicated behavior. Eventually we reach 4-periodic solutions, and then 8-periodic solutions, and then even higher periods.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs7, ax7 = plt.subplots(1, 2, figsize=(10, 4))
cwx7 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy7 = np.zeros_like(cwx7) # cobweb y-coordinates

# Draw plot lines
@widgets.interact(r=(2.25, 3.0, 0.01))
def update7(r=2.5):
    global ax7, cwx7, cwy7
    
    # Cobweb plot
    ax7[0].clear()
    ax7[0].set_xlim([0, 1.5])
    ax7[0].set_ylim([0, 1.5])
    ax7[0].grid(False)
    ax7[0].set_title("Cobweb Plot")
    ax7[0].set_xlabel("$x_n$")
    ax7[0].set_ylabel("$x_{n+1}$")
    cobweb_update(0.25, cwx7, cwy7, r=r, mode=0)
    ax7[0].plot(x, lmap(x, r=r), color="C0")
    ax7[0].plot(x, x, color="black")
    ax7[0].plot(cwx7[:2*COBWEB_MAX+1], cwy7[:2*COBWEB_MAX+1], color="C1")
    
    # Scatter plot
    ax7[1].clear()
    ax7[1].set_ylim([0, 1.5])
    ax7[1].grid(False)
    ax7[1].set_title("Time Series")
    ax7[1].set_xlabel("$n$")
    ax7[1].set_ylabel("$x_n$")
    ax7[1].plot(np.append([0], nval[1:2*COBWEB_MAX+1:2]), np.append([0.25], cwy7[1:2*COBWEB_MAX+1:2]), color="C0", marker=".", markersize=10)<jupyter_output><empty_output><jupyter_text>Eventually when $r$ becomes large enough the behavior of the system becomes impossible to describe in terms of simple periodicity and we enter a region of chaos.Chaotic behavior often arises in deceptively simple dynamical systems models. Like the simple systems seen earlier, chaotic systems are still deterministic in the sense that they produce exactly the same output when given exactly the same input, but they are extremely sensitive to minor changes in initial conditions.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs8, ax8 = plt.subplots(1, 2, figsize=(10, 4))
cwx8 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy8 = np.zeros_like(cwx8) # cobweb y-coordinates

# Draw plot lines
@widgets.interact(step=(0, 2*COBWEB_MAX, 1), r=(0.5, 3.0, 0.01), x0=(0.0, 1.25, 0.01))
def update8(step=4, r=1.5, x0=0.25):
    global ax8, cwx8, cwy8
    
    # Cobweb plot
    ax8[0].clear()
    ax8[0].set_xlim([0, 1.5])
    ax8[0].set_ylim([0, 1.5])
    ax8[0].grid(False)
    ax8[0].set_title("Cobweb Plot")
    ax8[0].set_xlabel("$x_n$")
    ax8[0].set_ylabel("$x_{n+1}$")
    cobweb_update(x0, cwx8, cwy8, r=r, mode=0)
    ax8[0].plot(x, lmap(x, r=r), color="C0")
    ax8[0].plot(x, x, color="black")
    ax8[0].plot(cwx8[:step+2], cwy8[:step+2], color="C1")
    ax8[0].plot(cwx8[step:step+2], cwy8[step:step+2], color="red")
    
    # Scatter plot
    ax8[1].clear()
    ax8[1].set_ylim([0, 1.5])
    ax8[1].grid(False)
    ax8[1].set_title("Time Series")
    ax8[1].set_xlabel("$n$")
    ax8[1].set_ylabel("$x_n$")
    ax8[1].plot(np.append([0], nval[1:step+2:2]), np.append([x0], cwy8[1:step+2:2]), color="C0", marker=".", markersize=10)
    ax8[1].plot(nval[step+1], cwy8[step+1], color="red", marker=".", markersize=10)<jupyter_output><empty_output><jupyter_text>Click and drag the sliders above to experiment with different growth rates $r$ and initial populations $x_0$, and observe how the population changes as the time steps advance.  Cobweb Diagram with Multiple EquilibriaSince cobweb diagrams can be used to estimate the behavior of a system around some equilibrium value, they can also be a valuable tool for classifying different equilibria. Solutions should generally approach stable equilibria and move away from unstable ones, and so if a system has multiple equilibria, a quick way to classify them is to draw several cobweb diagrams for different choices of initial solution.For example, consider a population model for a species with both an upper bound $U$ (representing, for example, resource limitations) _and_ a lower bound $L$ (representing, for example, a minimum breeding population) which takes the form\begin{align*}\Delta x_n = r x_n (L - x_n) (U - x_n)\end{align*}This model has three equilibria: $x^* = 0$, $x^* = L$, and $x^* = U$, since setting $x_n$ to any of these values results in $\Delta x_n = 0$. It's not immediately obvious which are stable and which are unstable, but a cobweb diagram can help.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs9, ax9 = plt.subplots(1, 2, figsize=(10, 4))
cwx9 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy9 = np.zeros_like(cwx9) # cobweb y-coordinates

# Draw plot lines
@widgets.interact(x0=(0.0, 1.25, 0.01))
def update9(x0=0.65):
    global ax9, cwx9, cwy9
    
    # Cobweb plot
    ax9[0].clear()
    ax9[0].set_xlim([0, 1.5])
    ax9[0].set_ylim([0, 1.5])
    ax9[0].grid(False)
    ax9[0].set_title("Cobweb Plot")
    ax9[0].set_xlabel("$x_n$")
    ax9[0].set_ylabel("$x_{n+1}$")
    cobweb_update(x0, cwx9, cwy9, r=1.5, mode=1)
    ax9[0].plot(x, tmap(x, r=1.5), color="C0")
    ax9[0].plot(x, x, color="black")
    ax9[0].plot(cwx9[:2*COBWEB_MAX+1], cwy9[:2*COBWEB_MAX+1], color="C1")
    
    # Scatter plot
    ax9[1].clear()
    ax9[1].set_ylim([0, 1.5])
    ax9[1].grid(False)
    ax9[1].set_title("Time Series")
    ax9[1].set_xlabel("$n$")
    ax9[1].set_ylabel("$x_n$")
    ax9[1].plot(np.append([0], nval[1:2*COBWEB_MAX+1:2]), np.append([x0], cwy9[1:2*COBWEB_MAX+1:2]), color="C0", marker=".", markersize=10)<jupyter_output><empty_output><jupyter_text>Click and drag to change the initial population $x_0$.In this case we see that selecting an initial population just below or just above $U$ causes the solution to converge to $U$, making $x^* = U$ a stable equilibrium. Selecting an initial population on either side of $L$ causes the solution to diverge from $L$, making $x^* = L$ unstable. Between $0$ and $L$ the population decays to zero, and so $x^* = 0$ is also stable in this model, unlike in the logistic growth model, where it was unstable.<jupyter_code># Set up side-by-side plots and initialize cobweb coordinate lists
figs10, ax10 = plt.subplots(1, 2, figsize=(10, 4))
cwx10 = np.zeros(2*COBWEB_MAX+2) # cobweb x-coordinates
cwy10 = np.zeros_like(cwx10) # cobweb y-coordinates

# Draw plot lines
@widgets.interact(step=(0, 2*COBWEB_MAX, 1), r=(0.5, 4.0, 0.01), x0=(0.0, 1.25, 0.01))
def update10(step=4, r=1.5, x0=0.65):
    global ax10, cwx10, cwy10
    
    # Cobweb plot
    ax10[0].clear()
    ax10[0].set_xlim([0, 1.5])
    ax10[0].set_ylim([0, 1.5])
    ax10[0].grid(False)
    ax10[0].set_title("Cobweb Plot")
    ax10[0].set_xlabel("$x_n$")
    ax10[0].set_ylabel("$x_{n+1}$")
    cobweb_update(x0, cwx10, cwy10, r=r, mode=1)
    ax10[0].plot(x, tmap(x, r=r), color="C0")
    ax10[0].plot(x, x, color="black")
    ax10[0].plot(cwx10[:step+2], cwy10[:step+2], color="C1")
    ax10[0].plot(cwx10[step:step+2], cwy10[step:step+2], color="red")
    
    # Scatter plot
    ax10[1].clear()
    ax10[1].set_ylim([0, 1.5])
    ax10[1].grid(False)
    ax10[1].set_title("Time Series")
    ax10[1].set_xlabel("$n$")
    ax10[1].set_ylabel("$x_n$")
    ax10[1].plot(np.append([0], nval[1:step+2:2]), np.append([x0], cwy10[1:step+2:2]), color="C0", marker=".", markersize=10)
    ax10[1].plot(nval[step+1], cwy10[step+1], color="red", marker=".", markersize=10)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14100474739601784, 'fraction_numerical': 0.04701339190816977, 'mean_word_length': 4.432448036951501, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '408433', 'license': 'MIT', 'n_tokens_mistral': 10482, 'n_tokens_neox': 10065, 'n_words': 3289, 'path': 'calc-diffeq-analysis/cobwebbing.ipynb', 'repo_name': 'adam-rumpf/jupyter-class-demonstrations', 'size': '28230'}