TP 3: Numpy, Matplotlib.¶
On commence par importer les librairies
import numpy as np
from matplotlib import pyplot as plt
On représente la fonction sur 4 périodes.
x = np.linspace(-4*np.pi, 4*np.pi, 1000)
Comme $f$ est $2\pi$-périodique, on a $f(x+2k\pi)=f(x)$ pour tout $x\in\mathbb{R}$. Pour connaitre les valeurs de $x$ par, il suffit de connaitre la valeur de $f$ au point $x$ modulo $2\pi$.
Ce qui nous donne
y = [x**2 for x in np.mod(x,2*np.pi)]
Ici np.mod(x,2*np.pi) renvoie un numpy array des valeurs de x modulo $2\pi$.
plt.plot(x, y, 'k')
plt.xlabel('x')
plt.title(r'tracer de la fct 2$\pi$ périodique définie par $f(x)=x^2$ sur $[0,2\pi[$')
plt.show()
Une autre façon est de calculer les valeurs de $f$ sur une période et de répeter cette liste auant de fois que l'on veut de période de $f$.
y = [x**2 for x in np.arange(0,2*np.pi + 0.01, 0.01)]*4
x = np.linspace(-4*np.pi, 4*np.pi, len(y)) # abscisse pour le graph représentant 4 prériodes
plt.plot(x, y, 'k')
plt.xlabel('x')
plt.title(r'tracer de la fct 2$\pi$ périodique définie par $f(x)=x^2$ sur $[0,2\pi[$')
plt.show()
x = np.arange(0, 2 + 0.01, 0.01) # les abscisses
A = [0.5, 1, 2, 3] # les différentes valeurs de a
for a in A:
plt.plot(x, x**a, label = "a = "+str(a)) # plot de chaque courbe pour les différents a
plt.xlabel('x')
plt.legend(loc = 'upper left')
plt.title(r'tracer de $x\mapsto x^a$')
plt.show()
1ère solution: fonction récursive.¶
def U(n, u0):
if n == 0:
return u0
else:
return -0.5 * U(n-1, u0) + 63
u0 = -8.41
Seq = [U(n, u0) for n in range(12)]
plt.plot(Seq, 'ko')
plt.title(r'$u_{n+1} = -0.5 u_n + 63$ avec $u_0=-8.41$ et $n\in\{0,\dots,11\}$')
plt.show()
2ème solution: fonction non récursive¶
def V(n, u0):
v = [u0]
for j in range(1, n+1):
v.append(-0.5 * v[j-1] + 63)
return v
u0 = -8.41
n = 12
Seq = V(n, u0)
plt.plot(Seq, 'ko')
plt.title(r'$u_{n+1} = -0.5 u_n + 63$ avec $u_0=-8.41$ et $n\in\{0,\dots,11\}$')
plt.show()
3ème solution: Calcul direct¶
On calcul explicitement le terme général de la suite. Il s'agit d'une suite arithmético-géométrique
$$ u_{n+1} = a u_n +b.$$
Par récurrence on montre que
$$ u_n = a^n (u_0 - c) + c$$ avec $c = b/(1-a)$
On peut alors faire le calcul direct!
def W(n,u0, a, b):
c = b/(1-a)
return a**n * (u0 - c) + c
u0 = -8.41
Seq = [W(n, u0, -0.5, 63) for n in range(12)]
plt.plot(Seq, 'ko')
plt.title(r'$u_{n+1} = -0.5 u_n + 63$ avec $u_0=-8.41$ et $n\in\{0,\dots,11\}$')
plt.show()
Exercice 4¶
Nous cherchons à approximer $\pi$ en "lançant des fléchettes aléatoires et en comptant celles qui atteignent la cible". Pour cela: on simule $n$ points dont les coordonnées sont des nombres aléatoires uniformément distribués compris entre $-1$ et $1$. On compte parmi ces points ceux qui sont dans le cercle unité. La proportion de points tombés dans le cercle unité est une approximation de l’aire du cercle unité divisée par l’aire du carré donc de $\pi/4$.¶
Dans cet exercice on va avoir besoin du sous module random de numpy.
import numpy.random as npr
1) Écrire une fonction simulation qui prend en argument un entier $n>1$, qui simule $n$ points dont les coordonnées sont des nombres aléatoires uniformément distribués compris entre $-1$ et $1$ et qui renvoie la proportion de points dans le cercle unité multiplié par $4$. Par exemple si $n=100$, et que l’on a obtenu $80$ points dans le cercle unité parmi les 100 points simulés, la proportion de points dans le cercle unité est de $0.8$, simulation(100) renverra donc $4\times 0.8$, c’est-à-dire $3.2$.¶
On rappelle que (voir TP1) que True = 1 et False = 0. On peut donc écrire la fonction. Vous pouvez aussi tester dans la console sum([True, False]) pour comprendre.
def simulation(n):
V = npr.uniform(-1, 1, (2,n))
return 4 * sum(V[0,]**2 + V[1,]**2 < 1)/n
simulation(100000)
3.1346
2) Écrire un script qui exécute 100 fois l’instruction simulation(1000) et qui calcule la moyenne des $100$ estimations de $\pi$ obtenues.¶
On genère 100 approximations de $\pi$ et on prend la moyenne.
pi_approx = np.mean([simulation(1000) for _ in range(101)])
print('approximation = %s, pour une valeur plus précise de %s' % (pi_approx, np.pi))
approximation = 3.138217821782178, pour une valeur plus précise de 3.141592653589793
Le _ après le for signifie qu'on fait un boucle sans se préoccuper de la variable de cette boucle, puisqu'elle ne joue aucun rôle ici.
3) Représenter sur un même graphique le cercle unité ainsi que les points tirés au "hasard". Vous pouvez chercher en ligne comment dessiner un cercle.¶
def Simulation_disque(n):
V = npr.uniform(-1, 1, (2,n))
return V[:, V[0,]**2 + V[1,]**2 < 1 ] # Il s'agit d'un masque
# On ne garde que les colonnes
# pour lesquelles la condition est True
n = 2000
X = Simulation_disque(n)
plt.scatter(X[0, ], X[1, ], s = 5) # affichage des points aléatoires
disk1 = plt.Circle((0, 0), 1, color='k', fill=False) # affichage du cercle
plt.gca().add_artist(disk1)
plt.axis([-1.1, 1.1, -1.1, 1.1 ])
plt.xlabel(r'$X_1$')
plt.ylabel(r'$X_2$')
plt.title('Distribution aléatoire de point uniformément distribués sur le disque unité')
plt.show()
def prod_for_loop(n, A, B):
C = np.zeros((n,n))
for j in range(n):
for l in range(n):
c = 0
for r in range(n):
c += A[j,r]*B[r,l]
C[j,l] = c
return C
n = 150
A = npr.uniform(0, 1, (n,n))
B = npr.uniform(0, 1, (n,n))
%timeit prod_for_loop(n, A, B) # produit avec 3 boucles for
987 ms ± 19.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.dot(A,B) # produit optimisé de numpy
# la différence d'exécution est importante
88.2 µs ± 6.63 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit np.matmul(A,B)
88.6 µs ± 3.2 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
C'est tout de meme possible d'utiliser les boucles for de manière raisonnable
from numba import njit
@njit
def prod_for_loop_numba(n, A, B):
C = np.zeros((n,n))
for j in range(n):
for l in range(n):
c = 0
for r in range(n):
c += A[j,r]*B[r,l]
C[j,l] = c
return C
%timeit prod_for_loop_numba(n, A, B)
3.4 ms ± 59 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
On remarque que la différence de temps d'éxecution est énorme!
R = np.linspace(0, 4, 201)
2) Pour chaque $r\in R$ représenter les points de coordonnée $(r,u_{100}),\dots,(r,u_{200})$¶
Voici une fonction (non récursive ci-dessous) qui génére la liste les point de la suite.
def suite_log(n, r, u0):
u = u0
suite = []
for j in range(1, 201):
u = r * u * (1 - u)
suite.append(u)
return suite
u0 = 0.5
n0 = 100
n1 = 200
points = [ suite_log(200, r, 0.5) for r in R ] # liste de liste des points pour les différents r
for i,r in enumerate(R):
plt.plot(r*np.ones(100), points[i][100:], 'ko', markersize = 0.1)
# affichage des point 100 à 200 pour les différents r
plt.xlabel('r')
plt.ylabel(r'$(u_{100},\dots,u_{200})$')
plt.title(r'apparition du comportement chaotique de la suite $(u_n)$ qd $r$ grandit')
plt.show()
Consulter la page https://fr.wikipedia.org/wiki/Suite_logistique pour une explication du phénomène observé.
Si vous ne comprenez pas enumerate, allez voir les exmples du deuxième cours sur AMETICE
R = np.linspace(0, 5, 11)
for i,r in enumerate(R):
print(i, ' ', r)
0 0.0 1 0.5 2 1.0 3 1.5 4 2.0 5 2.5 6 3.0 7 3.5 8 4.0 9 4.5 10 5.0