Logistic Map (Single-Species Discrete-Time Model)

4112 days ago by medlock

Logistic Model

var('R0') f(x) = R0 * x * (1 - x) show(f(x)) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(x - 1\right)} R_{0} x
\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(x - 1\right)} R_{0} x

Find Equilibria

fEqs = solve(x == f(x), x) show(fEqs) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = \frac{R_{0} - 1}{R_{0}}, x = 0\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[x = \frac{R_{0} - 1}{R_{0}}, x = 0\right]

Check Stability

show(diff(f(x), x).subs(x = 0)) show(expand(diff(f(x), x).subs(x = (R0 - 1) / R0))) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}R_{0}
\newcommand{\Bold}[1]{\mathbf{#1}}-R_{0} + 2
\newcommand{\Bold}[1]{\mathbf{#1}}R_{0}
\newcommand{\Bold}[1]{\mathbf{#1}}-R_{0} + 2
@interact def Cobweb(R0 = slider(0, 4, 0.01, 1.5), x0 = slider(0, 1, 0.01, 0.1), t = (0..100)): width = 1 X = srange(-0.1 * width, (1 + 0.1) * width, width / 100) g = f.subs(R0 = R0) Y = [g(x) for x in X] plotF = line(zip(X, Y), color = 'red') plotXY = line(zip(X, X), color = 'black') plotStart = point([x0, 0], color = 'blue', pointsize = 30) xt = x0 xtVals = [xt] X = [xt] Y = [0] for tau in range(t): fxt = g(xt) xtVals.append(fxt) X.extend([xt, fxt]) Y.extend([fxt, fxt]) xt = fxt plotXt = line(zip(X, Y), color = 'blue') plotEnd = point([X[-1], Y[-1]], color = 'blue', pointsize = 30) show(line(zip(range(0, t + 1), xtVals), marker = 'o'), axes_labels = ['$t$', '$x$'], xmin = 0) show(plotXY + plotF + plotXt + plotStart + plotEnd, axes_labels = ['$x_t$', '$x_{t+1}$']) 
       
R0 
x0 

Click to the left again to hide and once more to show the dynamic interactive window

Bifurcation Diagram

%cython cpdef Trajectory(double x0, long tEnd, long t1, double R0): cdef double x = x0 cdef long tau for tau from 1 <= tau <= t1: x = R0 * x * (1 - x) xvals = [x] for tau from t1 + 1 <= tau <= tEnd: x = R0 * x * (1 - x) xvals.append(x) return xvals cdef extern from "math.h": double log(double) cpdef Lyapunov(double x0, long tEnd, double R0): cdef double x = x0 cdef long tau S = 0 for tau from 1 <= tau <= tEnd: x = R0 * x * (1 - x) S += log(abs(R0 * (1 - 2 * x))) return S / tEnd 
@interact def BifurcationDiagram(R0Min = slider(1, 4, 0.001, 1), R0Max = slider(1, 4, 0.001, 4)): dR0 = (R0Max - R0Min) / 1000.0 xpts = [] Lpts = [] x = 0.5 for R0 in srange(R0Min, R0Max, dR0): # Iterate the model for 1100 times steps, # discarding the first 1000 xs = Trajectory(x, 1100, 1000, R0) xpts.extend([[R0, x] for x in xs]) x = xs[-1] # Update the initial condition LE = Lyapunov(0.05, 1000, R0) Lpts.append([R0, LE]) show(points(xpts, pointsize = 1), axes_labels = ['$R_0$', '$x^*$']) show(points(Lpts, pointsize = 1), axes_labels = ['$R_0$', '$\lambda$']) 
       
R0Min 
R0Max 

Click to the left again to hide and once more to show the dynamic interactive window

Sensitive Dependence on Initial Conditions

xt1 = Trajectory(0.05, 50, 0, 4) xt2 = Trajectory(0.050001, 50, 0, 4) show(line(zip(range(len(xt1)), xt1), color = 'blue', marker = 'o') + line(zip(range(len(xt2)), xt2), color = 'red', marker = 'o'), axes_labels = ['$t$', '$x$'], xmin = 0)