Week 6 Bifurcations and Chaos

Logistic Equation

Bunny Growth on Magical Island $$ X_{n+1} = rX_n$$

On a non-magical island $$X_{n+1} = rX_n(1-X_n)$$

In [1]:
import matplotlib.pyplot as plt
import numpy as np
import Pythonic_files as pf

%config InlineBackend.figure_format='retina'
In [16]:
plt.figure(figsize=(13,5))

r0=3.1

X_0=0.2
Steps=100
Max_x=1.1
X_old=X_0

for t in range(0,Steps):
    X_new=r0*X_old*(1-X_old)
    plt.plot(X_old,X_new,".",c='r',alpha=0.3)
    pf.Cobweb(X_old,X_new,'g',0.2)
    
    X_old=X_new
    Max_x=max(Max_x,X_new)
    
X=np.linspace(0,Max_x,100)
Y=r0*X*(1-X)
plt.plot(X,Y,c='b',alpha=0.5)

plt.plot(X,X,c='r',alpha=0.4)
plt.ylim(0,Max_x)
Out[16]:
(0.0, 1.1)
In [13]:
def log_func(x,r):
    return r*x*(1-x)

What happens when we vary $r$

In [19]:
for R in np.linspace(2,3.5,10):
    plt.figure(figsize=(13,5))

    r0=R

    X_0=0.2
    Steps=100
    Max_x=1.1
    X_old=X_0

    for t in range(0,Steps):
        X_new=log_func(X_old,r0)
        plt.plot(X_old,X_new,".",c='r',alpha=0.3)
        pf.Cobweb(X_old,X_new,'g',0.2)
    
        X_old=X_new
        Max_x=max(Max_x,X_new)
    
    X=np.linspace(0,Max_x,100)
    Y=log_func(X,r0)
    plt.plot(X,Y,c='b',alpha=0.5)

    plt.plot(X,X,c='r',alpha=0.4)
    plt.ylim(0,Max_x)
    plt.title('Parameter_value is r = '+str(R))
In [20]:
for R in np.linspace(2.95,3.6,10):
    plt.figure(figsize=(13,5))

    r0=R

    X_0=0.2
    Steps=100
    Max_x=1.1
    X_old=X_0

    for t in range(0,Steps):
        X_new=log_func(X_old,r0)
        plt.plot(X_old,X_new,".",c='r',alpha=0.3)
        pf.Cobweb(X_old,X_new,'g',0.2)
    
        X_old=X_new
        Max_x=max(Max_x,X_new)
    
    X=np.linspace(0,Max_x,100)
    Y=log_func(X,r0)
    plt.plot(X,Y,c='b',alpha=0.5)

    plt.plot(X,X,c='r',alpha=0.4)
    plt.ylim(0,Max_x)
    plt.title('Parameter_value is r = '+str(R))
In [33]:
r=3.21

X=np.linspace(0,Max_x,100)
Y=log_func(X,r)
YY=log_func(log_func(X,r),r)
plt.plot(X,Y,c='b',alpha=0.5)
plt.plot(X,YY,c='g',alpha=0.5)

plt.plot(X,X,c='r',alpha=0.4)
plt.ylim(0,Max_x)
plt.title('Parameter_value is r = '+str(r))
plt.show()
In [34]:
for R in np.linspace(2.9,3.4,10):
    plt.figure(figsize=(13,5))

    r0=R

    X_0=0.2
    Steps=100
    Max_x=1.1
    X_old=X_0

    for t in range(0,Steps):
        X_new=log_func(X_old,r0)
        plt.plot(X_old,X_new,".",c='r',alpha=0.3)
        pf.Cobweb(X_old,X_new,'g',0.2)
    
        X_old=X_new
        Max_x=max(Max_x,X_new)
    
    X=np.linspace(0,Max_x,100)
    Y=log_func(X,r0)
    plt.plot(X,Y,c='b',alpha=0.5)
    
    YY=log_func(log_func(X,r0),r0)
    plt.plot(X,YY,c='g',alpha=0.5)
    plt.plot(X,X,c='r',alpha=0.4)
    plt.ylim(0,Max_x)
    plt.title('Parameter_value is r = '+str(R))
In [35]:
for R in np.linspace(3.4,3.6,3):
    plt.figure(figsize=(13,5))

    r0=R

    X_0=0.2
    Steps=100
    Max_x=1.1
    X_old=X_0

    for t in range(0,Steps):
        X_new=log_func(X_old,r0)
        plt.plot(X_old,X_new,".",c='r',alpha=0.3)
        pf.Cobweb(X_old,X_new,'g',0.2)
    
        X_old=X_new
        Max_x=max(Max_x,X_new)
    
    X=np.linspace(0,Max_x,100)
    Y=log_func(X,r0)
    plt.plot(X,Y,c='b',alpha=0.5)
    
    YY=log_func(log_func(X,r0),r0)
    plt.plot(X,YY,c='g',alpha=0.5)
    plt.plot(X,X,c='r',alpha=0.4)
    plt.ylim(0,Max_x)
    plt.title('Parameter_value is r = '+str(R))
In [2]:
def log_multi_func(x,r,n):
    def mini_func(x,r):
        return r*x*(1-x)
    for count in range(0,n):
        x = mini_func(x,r)
    return x
In [49]:
r=3.5

X=np.linspace(0,Max_x,100)
Y=log_func(X,r)
YY=log_multi_func(X,r,2)
YYY=log_multi_func(X,r,4)

plt.plot(X,Y,c='b',alpha=0.5)
plt.plot(X,YY,c='g',alpha=0.5)
plt.plot(X,YYY,c='y',alpha=0.5)

plt.plot(X,X,c='r',alpha=0.4)
plt.ylim(0,Max_x)
plt.title('Parameter_value is r = '+str(r))
plt.show()
In [51]:
for R in np.linspace(3.5,3.6,1):
    plt.figure(figsize=(13,5))

    r0=R

    X_0=0.2
    Steps=100
    Max_x=1.1
    X_old=X_0

    for t in range(0,Steps):
        X_new=log_func(X_old,r0)
        plt.plot(X_old,X_new,".",c='r',alpha=0.3)
        pf.Cobweb(X_old,X_new,'g',0.2)
    
        X_old=X_new
        Max_x=max(Max_x,X_new)
    
    X=np.linspace(0,Max_x,1000)
    Y=log_func(X,r0)
    plt.plot(X,Y,c='b',alpha=0.5)
    
    YY=log_multi_func(X,r0,4)
    plt.plot(X,YY,c='g',alpha=0.5)
    plt.plot(X,X,c='r',alpha=0.4)
    plt.ylim(0,Max_x)
    plt.title('Parameter_value is r = '+str(R))

Lets look at these periodic solutions for different $r$ values

In [81]:
n = 1000
x0=0.2

R = np.linspace(2,3.65,n)
X = x0*np.ones(n)

YY = log_multi_func(X,R,1000)
In [88]:
for count in range(0,200):
    plt.plot(R,YY,',k',c='b',alpha=0.01)
    YY=log_func(YY,R)
In [93]:
plt.figure(figsize=(13,5))

n = 1000
x0=0.2

R = np.linspace(2,4,n)
X = x0*np.ones(n)

YY = log_multi_func(X,R,5000)


for count in range(0,600):
    plt.plot(R,YY,',k',c='black',alpha=0.1)
    YY=log_func(YY,R)
In [94]:
plt.figure(figsize=(13,5))

n = 1000
x0=0.1

R = np.linspace(2,4,n)
X = x0*np.ones(n)

YY = log_multi_func(X,R,5000)


for count in range(0,600):
    plt.plot(R,YY,',k',c='black',alpha=0.1)
    YY=log_func(YY,R)
In [97]:
rleft = 3.65
rright = 3.7

plt.figure(figsize=(13,5))

n = 1000
x0=0.1

R = np.linspace(rleft,rright,n)
X = x0*np.ones(n)

YY = log_multi_func(X,R,5000)


for count in range(0,1000):
    plt.plot(R,YY,',k',c='black',alpha=0.1)
    YY=log_func(YY,R)
In [98]:
rleft = 2.6
rright = 4

plt.figure(figsize=(13,5))

n = 1000
x0=0.1

R = np.linspace(rleft,rright,n)
X = x0*np.ones(n)

YY = log_multi_func(X,R,5000)


for count in range(0,1000):
    plt.plot(R,YY,',k',c='black',alpha=0.1)
    YY=log_func(YY,R)

Chaos and chaotic encoding

In [4]:
def log_func(x,r):
    return r*x*(1-x)
In [8]:
plt.figure(figsize=(13,5))

R=3.9
r0=R

X_0=0.2
Steps=1000
Max_x=1.1
X_old=X_0

for t in range(0,Steps):
    X_new=log_func(X_old,r0)
    plt.plot(X_old,X_new,".",c='r',alpha=0.3)
    pf.Cobweb(X_old,X_new,'g',0.2)
    
    X_old=X_new
    Max_x=max(Max_x,X_new)
    
X=np.linspace(0,Max_x,1000)
Y=log_func(X,r0)
plt.plot(X,Y,c='b',alpha=0.5)
    


plt.plot(X,X,c='r',alpha=0.4)
plt.ylim(0,Max_x)
plt.title('Parameter_value is r = '+str(R))
plt.show()
In [28]:
#plt.figure(figsize=(13,5))

R=3.9
r0=R

X_0=0.2
Steps=1000

X_old=X_0

Points=[]
for t in range(0,Steps):
    Points.append(X_old)
    X_new=log_func(X_old,r0)
    X_old=X_new
In [29]:
plt.figure(figsize=(13,5))

X = np.arange(0,Steps)


plt.plot(X[-100:],Points[-100:],'.-',alpha=0.3)
Out[29]:
[<matplotlib.lines.Line2D at 0x7fa87970c880>]
In [39]:
plt.figure(figsize=(13,5))

R=3.9
r0=R

X_0=0.2
Delta = 0.01
Steps=100000

X_old=X_0
X_old2 = X_0+Delta

Points=[]
Points2=[]

for t in range(0,Steps):
    Points.append(X_old)
    Points2.append(X_old2)
    
    X_new=log_func(X_old,r0)
    X_new2=log_func(X_old2,r0)
    
    X_old=X_new
    X_old2=X_new2




X = np.arange(0,Steps)
Window = 50

plt.plot(X[-Window:],Points[-Window:],'b',alpha=0.3)
plt.plot(X[-Window:],Points2[-Window:],'r',alpha=0.3)
Out[39]:
[<matplotlib.lines.Line2D at 0x7fa879244a00>]
In [53]:
D = np.linspace(0,16,1000)
Msg = np.cos(D)/100
plt.plot(D,Msg,'.')
plt.show()
In [48]:
plt.plot(D,Msg,',')
Out[48]:
[<matplotlib.lines.Line2D at 0x7fa898b40eb0>]
In [73]:
Startpoint = np.random.randint(1000,90000)
Startpoint
Out[73]:
4826
In [78]:
Coded_Msg = Msg + Points[Startpoint:Startpoint+len(Msg)]
In [79]:
plt.plot(D,Msg,'.')
plt.plot(D,Coded_Msg,'.')
Out[79]:
[<matplotlib.lines.Line2D at 0x7fa8ac7d0070>]
In [80]:
Decoded_Msg = Coded_Msg - Points[Startpoint:Startpoint+len(Msg)]
In [81]:
plt.plot(D,Decoded_Msg,'.')
#plt.plot(D,Coded_Msg,'.')
Out[81]:
[<matplotlib.lines.Line2D at 0x7fa8acb895e0>]

Julia Set and Mandelbrot set

$$ X_{n+1}=X_n^2 +C $$
In [102]:
def func(x,c):
    return x**2 +c

plt.figure(figsize=(13,5))

C=-1.15
c=C

X_0=0.2
Steps=50
Max_x=1.1
Min_x=-1.1
X_old=X_0

for t in range(0,Steps):
    X_new=func(X_old,c)
    plt.plot(X_old,X_new,".",c='r',alpha=0.3)
    pf.Cobweb(X_old,X_new,'g',0.2)
    Max_x=max(Max_x,X_new)
    Min_x=min(Min_x,X_new)
    X_old=X_new
    
    
X=np.linspace(Min_x,Max_x,1000)
Y=func(X,c)
plt.plot(X,Y,c='b',alpha=0.5)
    


plt.plot(X,X,c='r',alpha=0.4)
#plt.ylim(0,Max_x)
plt.title('Parameter_value is C = '+str(C))
plt.show()
In [113]:
plt.figure(figsize=(13,5))

def multi_func(x,c,n):
    def mini_func(x,c):
        return x**2+c
    for count in range(0,n):
        x = mini_func(x,c)
    return x
    
    

n =1000
x0=0

C = np.linspace(-2,0,n)
X = x0*np.ones(n)

YY = multi_func(X,C,5000)

for count in range(0,600):
    plt.plot(C,YY,',k',c='black',alpha=0.1)
    YY = func(YY,C)

Lets look at this model in the Complex plane.

In [121]:
X = np.linspace(-2,2,4)
Y = np.linspace(-2,2,4)

X,Y = np.meshgrid(X,Y)
def Z(x,y):
    return x+1j*y
In [122]:
C=Z(X,Y)
In [136]:
Z = complex(0,0)
for count in range(3):
    Z=Z**2 + C
    Points= (abs(Z)<10)
In [137]:
Points
Out[137]:
array([[False, False, False, False],
       [ True,  True,  True, False],
       [ True,  True,  True, False],
       [False, False, False, False]])
In [131]:
C
Out[131]:
array([[-2.        -2.j        , -0.66666667-2.j        ,
         0.66666667-2.j        ,  2.        -2.j        ],
       [-2.        -0.66666667j, -0.66666667-0.66666667j,
         0.66666667-0.66666667j,  2.        -0.66666667j],
       [-2.        +0.66666667j, -0.66666667+0.66666667j,
         0.66666667+0.66666667j,  2.        +0.66666667j],
       [-2.        +2.j        , -0.66666667+2.j        ,
         0.66666667+2.j        ,  2.        +2.j        ]])
In [141]:
Z = complex(0,0)
for count in range(5):
    Z=Z**2 + C
    Points= (abs(Z)<10)
    plt.imshow(Points)
In [151]:
X = np.linspace(-2,0.5,1000)
Y = np.linspace(-1.5,1.5,1000)

X,Y = np.meshgrid(X,Y)
def Z(x,y):
    return x+1j*y

C=Z(X,Y)

Z = complex(0,0)
for count in range(50):
    Z=Z**2 + C
    Points= (abs(Z)<20)
plt.imshow(Points)

<ipython-input-151-b10165bb1006>:13: RuntimeWarning: overflow encountered in absolute
  Points= (abs(Z)<20)
<ipython-input-151-b10165bb1006>:12: RuntimeWarning: overflow encountered in square
  Z=Z**2 + C
<ipython-input-151-b10165bb1006>:12: RuntimeWarning: invalid value encountered in square
  Z=Z**2 + C
<ipython-input-151-b10165bb1006>:13: RuntimeWarning: invalid value encountered in less
  Points= (abs(Z)<20)
Out[151]:
<matplotlib.image.AxesImage at 0x7fa8497e0070>
In [153]:
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        X = np.linspace(-2,0.5,1000)
        Y = np.linspace(-1.5,1.5,1000)

        X,Y = np.meshgrid(X,Y)
        def Z(x,y):
            return x+1j*y

        C=Z(X,Y)

        Z = complex(0,0)
        for count in range(100):
            Z=Z**2 + C
            Points= (abs(Z)<20)

        
        plt.imshow(Points)
In [154]:
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        ydown,yup=-1.5,1.5
        xleft,xright=-2,0.5
        
        X = np.linspace(xleft,xright,1000)
        Y = np.linspace(ydown,yup,1000)

        X,Y = np.meshgrid(X,Y)
        def Z(x,y):
            return x+1j*y

        C=Z(X,Y)

        Z = complex(0,0)
        for count in range(100):
            Z=Z**2 + C
            Points= (abs(Z)<20)

        
        plt.imshow(Points,extent=[xleft,xright,ydown,yup])
In [159]:
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        ydown,yup=-.7,-.699
        xleft,xright=-.22,-.219
        
        X = np.linspace(xleft,xright,1000)
        Y = np.linspace(ydown,yup,1000)

        X,Y = np.meshgrid(X,Y)
        def Z(x,y):
            return x+1j*y

        C=Z(X,Y)

        Z = complex(0,0)
        for count in range(100):
            Z=Z**2 + C
            Points= (abs(Z)<20)

        
        plt.imshow(Points,cmap='magma',extent=[xleft,xright,ydown,yup])
In [161]:
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        ydown,yup=-.7,-.699
        xleft,xright=-.22,-.219
        
        X = np.linspace(xleft,xright,3000)
        Y = np.linspace(ydown,yup,3000)

        X,Y = np.meshgrid(X,Y)
        def Z(x,y):
            return x+1j*y

        C=Z(X,Y)

        Z = complex(0,0)
        for count in range(100):
            Z=Z**2 + C
            Points= (abs(Z)<20)

        pic = plt.figure()
        pic.set_size_inches((12,12))
        
        plt.imshow(Points,cmap='magma',extent=[xleft,xright,ydown,yup])
In [163]:
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        ydown,yup=-.71,-.689
        xleft,xright=-.23,-.209
        
        X = np.linspace(xleft,xright,3000)
        Y = np.linspace(ydown,yup,3000)

        X,Y = np.meshgrid(X,Y)
        def Z(x,y):
            return x+1j*y

        C=Z(X,Y)

        Z = complex(0,0)
        for count in range(100):
            Z=Z**2 + C
            Points= (abs(Z)<2)

        pic = plt.figure()
        pic.set_size_inches((12,12))
        
        plt.imshow(Points,cmap='magma',extent=[xleft,xright,ydown,yup])
In [165]:
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        ydown,yup=-.71,-.689
        xleft,xright=-.23,-.209
        n=3000
        
        X = np.linspace(xleft,xright,n)
        Y = np.linspace(ydown,yup,n)

        X,Y = np.meshgrid(X,Y)
        def Z(x,y):
            return x+1j*y

        C=Z(X,Y)

        Z = complex(0,0)
        
        FinalPic = np.zeros((n,n))
        
        for count in range(100):
            Z=Z**2 + C
            Points= (abs(Z)<2)
            
            FinalPic = FinalPic + Points*np.ones((n,n))

        pic = plt.figure()
        pic.set_size_inches((12,12))
        
        plt.imshow(FinalPic,cmap='magma',extent=[xleft,xright,ydown,yup])
$$ Z_{n+1} = Z_n +C $$

Mandelbrot Set is where you vary the "C" and look at the solutions starting at 0+0i

Julia Set, it varys the starting point of Z and keeps C constant

In [168]:
#C1= -0.8 + 0.156j
C1= -0.82 + 0.156j
n=1000


ydown, yup = -1.5,1.5
xleft,xright =-2,2


X=np.linspace(xleft,xright,n)
Y=np.linspace(ydown, yup,n)

X, Y = np.meshgrid(X, Y)
def Z(x,y):
    return x + 1j * y
Z0 = np.full((n, n), C1)
Final =np.zeros((n,n))
with np.warnings.catch_warnings():
        np.warnings.simplefilter("ignore")
        Z = Z(X, Y)
        for count in range(205):
            Z=Z*Z + Z0
            Points = (abs(Z)<2)
            Final = Final +Points*np.ones((n,n))

pic = plt.figure()
pic.set_size_inches(12,12)
ax = pic.add_axes([0, 0, 1, 1], frameon=False, aspect=1)
ax.set_xticks([])
ax.set_yticks([])
          
plt.imshow(Final,cmap = 'seismic')
plt.show()
In [ ]: