Diagramme de Feigenbaum

Diagramme de Feigenbaum

Vidéo de Veritasium : This equation will change how you see the world

import numpy as np
import matplotlib.pyplot as plt

## Question 1
"""
Pour tout x entre 0 et 1, 0 <= µx(1-x) <= µ/4 <= 1

On conclut par récurrence
"""

## Suite logistique

def logistique(mu,x0,n):
    """
    Renvoie la liste contenant les n premieres valeurs de la suite logistique de parametre mu et de premier terme x0
    """
    res = [x0]
    for i in range(n-1):
        res.append(mu*res[-1]*(1-res[-1]))
    return res

## Representation graphique d'une suite

def trace_suite(f,x0,n,N=1000):
    plt.clf()
    x = np.linspace(0,1,N)
    y = f(x)
    plt.plot(x,y) #On trace f
    plt.plot(x,x) #On trace la premiere bissectrice
    suite = [x0]
    for i in range(n-1):
        suite.append(f(suite[-1]))    
    abs = []
    ord = [x0]
    for i in range(n-1):
        abs.append(suite[i])
        abs.append(suite[i])
        ord.append(suite[i+1])
        ord.append(suite[i+1])
    ord.pop(-1) #Pour que abs et ord aient la meme taille
    plt.plot(abs,ord) #On trace la suite
    plt.draw(); plt.show()

# trace_suite(lambda x: 3*x*(1-x),.9,100)

r = np.linspace(1,4,200)
for mu in r:
    plt.clf()
    plt.title('µ='+str(mu))
    trace_suite(lambda x: mu*x*(1-x),.9,100)
    plt.pause(.01)

## Diagramme de bifurcation

# plt.plot(x,y,linestyle='',marker=',')

def bifurcation(x0=.9,rmin=2,rmax=4,N=1000):
    mus = np.linspace(rmin,rmax,N)
    abs = []
    ord = []
    plt.clf()
    for mu in mus:
        x = x0
        for i in range(500):
            if i>400:
                abs.append(mu)
                ord.append(x)
            x = mu*x*(1-x)
    plt.plot(abs,ord,marker=',',linestyle='')
    plt.draw(); plt.show()

## Diagrammes de phase

def phase(mu,x0,n):
    plt.clf()
    suite = logistique(mu,x0,n)
    abs = suite[:-1]
    ord = suite[1:]
    plt.axis([.2,.8,.8,1])
    plt.plot(abs,ord,marker=',',linestyle='')
    plt.draw(); plt.show()

def phases(x0,n):
    plt.clf()
    r = np.linspace(3.6,4,100)
    for mu in r:
        suite = logistique(mu,x0,n)
        plt.axis([0.2,.8,.8,1])
        plt.plot(suite[:-1],suite[1:],marker=',',linestyle='',color=plt.cm.rainbow((mu-3.5)/.5))
    plt.draw(); plt.show()