%pylab notebook

Populating the interactive namespace from numpy and matplotlib

2- A simple model

2.1 Motivation for the Logistic Equation

$$x_{t+1} = F(x_t)$$ $$x_{t+1} = a x_t (1-x_t)$$

def f(x,a):
    return a*x*(1-x)

2.2 Explorations of $1<a<3$

Numerics for $a=2$

x=np.array([0.81131431,0.13443,0.26575])
listx=[]
listt=[]
for t in range(30):
    listx.append(x)
    listt.append(t)
    x=f(x,2)
#a=2,10.0/3,3.82843
plot(listt,listx,'-o');

Fixed points </red>

$$x_{n+1} = f(x_n) = x_n \text{ for } x_n = x^{(1)} \Rightarrow f(x^{(1)}) = x^{(1)} \Rightarrow x^{(1)}= 1-1/a$$

2.3 Numerical Exploration for $a=3.2$

listx=[]
listt=[]
x=np.array([0.13,0.23,0.93])
for t in range(30):
    listx.append(x)
    listt.append(t)
    x=f(x,10./3)
#a=2,10.0/3,3.82843
plot(listt,listx,'-o');

Periodic orbit (period 2) analysis

</red>

$$x^{(2)} = f(f(x^{(2)})) \Rightarrow$$ $$x_\pm^{(2)} = \frac{1}{2a}\left((a+1)\pm\sqrt{a^2-2a-3}\right)$$

Two periodic points of the period 2 orbit: $x_- = f(x_+)$ and $x_+=f(x_-)$

def x2(a):
    return (a+1-np.sqrt(a*a-2*a-3))/(2*a)

f(x2(3.2),3.2)

0.7994554904673701

f(f(x2(3.2),3.2),3.2)

0.5130445095326298

2.4 Stability Analysis

A periodic orbit of period $n$ of the map $F$ is linearly stable if

$$|\lambda| \equiv |\prod_{i=1}^{n} F'(x_i) | < 1,$$

where $x_i$ is the $i$-th point of the orbit and $F'(x) \equiv \frac{dF}{dx}$. If $|\lambda| >1$ the orbit is unstable.

Stability analysis </red>

Fixed point $x^{(1)}$: $$\lambda = 2- a \Rightarrow |\lambda|< 1 \text{ for } 1< a <3$$ Periodic orbit of period 2, $x^{(2)}$: $$\lambda_2 = -a^2+2a+4 \Rightarrow |\lambda| < 1 \text{ for } 3< a < 1+\sqrt{6} \approx 3.45$$ </font>

2.5 Dependence on $a$ (parameter analysis)

aa=np.arange(1,4,.001)
a=np.arange(3,4,.001)
xlabel("a")
ylabel("x")
plot(aa,(1-1/aa),'--',color="blue")
plot(a,((1/(2*a))*(a+1+np.sqrt(a**2-2*a-3))),'--',color="red")
plot(a,((1/(2*a))*(a+1-np.sqrt(a**2-2*a-3))),'--',color="red")

[<matplotlib.lines.Line2D at 0x7ff4bd34fac8>]

Numerical Experiment

# Empty lists to collect results
lista=[]
listx=[]
#Loop for parameter a in a given range
for a in np.arange(1.,4.0,0.001):
    x=np.arange(0.001,1,0.01) #Initial conditions, uniformly distributed
    for t in range(500): # Iterate the map many times
        x=f(x,a)
    lista.append(np.full(len(x),a)) #List of values of a, convenient to plot
    listx.append(x) #At the end, store output of x in listx

#Plotting the results of simulations above
plot(lista,listx,'.',color="black",markersize=.15);

#Plotting the results of simulations above
xlabel("a")
ylabel("x")
plot(lista,listx,'.',color="black",markersize=.1);

2.6 Gaphical method

a=4
#a=4.0
x=np.arange(0,1,.001)
xlabel("$x_n$")
ylabel("$x_{n+1}$")
plot(x,f(x,a))
#plot(x,f(f(x,a),a))
#plot(x,f(f(f(x,a),a),a))
plot(x,x,"--",color="black");

# Birth of the period 2 orbit: period-doubling bifurcation
x=np.arange(0,1,.001)
a=2.9
plot(x,f(f(x,a),a),color="blue")
a=3.2
plot(x,f(f(x,a),a),color="red")
plot(x,x,"--",color="black")

[<matplotlib.lines.Line2D at 0x7ff4bec44668>]

# Birth of the period 3 orbit: saddle-node bifurcation
x=np.arange(0,1,.001)
a=3.82
plot(x,f(f(f(x,a),a),a),color="blue")
a=3.84
plot(x,f(f(f(x,a),a),a),color="red")
plot(x,x,"--",color="black")

[<matplotlib.lines.Line2D at 0x7ff4bc8987f0>]