%pylab notebook
def f(x,a):
return a*x*(1-x)
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');
$$ 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$$
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');
</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)
f(f(x2(3.2),3.2),3.2)
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.
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>
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")
# 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);
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")
# 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")