TP 4 : Modules mathématiques, partie 1 : Numpy et Matplotlib¶
Correction.
Exercice (Numpy et complexes) Ecrire une fonction
tab_complexe(N)qui renvoie un tableau Numpy de taille NxN, dont les éléments sont de typenp.complex, et tel que siA= tab_complexe(N)alorsA[x,y]est $x + i y$, où ici $i^2 = -1$. Sachez que pour créer des nombres complexes en Python, on peut utiliser la syntaxe (un peu curieuse)I= 1jpour queIsoit égal au complexe $i$, et ensuite fairex + y*I.Après avoir fait une première version, essayer les commandes :
L= np.arange(10)
np.stack([L,L,L])
np.stack([L,L,L], axis= 1)Utiliser alors ces commandes pour faire une version beaucoup plus rapide de
tab_complexe(N).
# voici une première version:
def tab_complexe(N):
A= np.zeros(shape=(N,N), dtype= np.complex)
I= 1j
for x in range(0,N):
for y in range(0,N):
A[x,y]= x + I*y
return A
tab_complexe(5)
Essayons les commandes en question :
L= np.arange(5, dtype= np.complex)
Y= np.stack([L,L,L,L,L])
Y
X= np.stack([L,L,L,L,L], axis= 1)
X
On comprend donc qu'il suffit de faire :
X+ 1j*Y
def tab_rapide(N):
L= np.arange(N, dtype= np.complex)
LL= [L for i in range(N)]
X= np.stack(LL, axis= 1)
Y= np.stack(LL)
return X + 1j*Y
%time tab_complexe(4000)
%time tab_rapide(4000)
Exercice (booléens Numpy 1) Ecrire une fonction
f(A)qui prend en paramètre un tableau NumpyA(n'importe quelle taille, n'importe quelle nombre de dimensions) et qui renvoie un tableauBqui est une copie deAdans laquelle les nombres $<0$ ont été remplacés par $0$ et les nombres $>1$ ont été remplacés par $1$. Vous pouvez faireA.copy()pour obtenir une copie deA.
def f(A):
B= A.copy()
B[ A < 0 ]= 0
B[ A > 1 ]= 1
return B
Exercice (booléens Numpy 2) Si
Aest un tableau Numpy, expliquez ce que donne la commande
A[ A<0 ]= - A[ A<0 ]et pourquoi.
Réponse: les coefficients de $A$ sont changés, et en fait, remplacés par leurs valeurs absolues.
A= np.array([[-4, 12], [0.5, 0.7]])
np.abs(A)
A[ A<0 ]= - A[ A<0 ]
A
Exercice (courbes de Lissajou) Pour diverses valeurs entières de $p$ et $q$, tracer la courbe paramétrée donnée par
$$ x(t) = \sin(pt), \qquad y(t) = \sin(qt) . $$
Il s'agit d'utiliser la commande
plt.plotpour tracer une ligne brisée qui sera une bonne approximation. On note qu'il suffit de prendre $t$ entre $0$ et $2\pi$.Au début, essayer de fixer $p=1, 2$ ou $3$, puis augmenter $q$.
p= 2
q= 3
N= 200
t= np.linspace(0, 7, N) # 7 = un peu plus de 2pi
x= np.sin( p*t )
y= np.sin( q*t )
plt.plot(x,y)
Exercice (triangle de Sierpinksi et nombres binomiaux) Faire d'abord une fonction
binomial(n,k)qui renvoie le nombre binomial$$ \frac {n!} {k!(n-k)!}$$
(On peut reprendre le code pour la factorielle du TP précédent.) Attention, pour $k > n$ on veut que la fonction renvoie $0$.
Puis, faire un tableau Numpy
Acarré, de taille $N\times N$ où $N$ est un entier que vous choisissez, tel queA[n,k] = binomial(n,k) % 2. Puis, affichezAcomme une image. Augmentez $N$.Enfin, essayez de remplacer le modulo 2 par un modulo $p$ pour un nombre premier $p$ que vous choisissez (ça marche un peu moins bien, mais on repère quand même des phénomènes).
def factorielle(n):
if n == 0:
return 1
else:
return n*factorielle(n-1)
def binomial(n,k):
if k > n or k < 0:
return 0
else:
return factorielle(n) // ( factorielle(k)*factorielle(n-k) )
factorielle(10)
N= 2**9
A= np.array([ [binomial(n,k)%2 for k in range(N)] for n in range(N) ])
# ou alors A= np.zeros(shape=(N,N)) puis on remplit A[n,k]= binomial(n,k)%2 ...
plt.imshow(A, cmap= "hot")
Exercice (ensemble de Julia) Soit $c$ un nombre complexe. L'ensemble de Julia au point $c$ est
$$ J_c= \{ w \in \mathbb{C} ~~|~~ \text{la suite définie par}~ z_0 = w, z_{n+1} = z_n^2 + c ~\text{reste bornée} \}.$$
Tracer l'ensemble de Julia pour diverses valeurs de $c$, en vous inspirant du code ci-dessus. Une remarque heuristique (qu'on peut partiellement justifier mathématiquement) est que si on prend une valeur $c$ telle que l'ensemble de Mandelbrot est "intéressant" (ou joli) autour de $c$, alors $J_c$ sera également intéressant. Là encore, faites partager vos plus belles images à la classe.
def temps_generalise(wx, wy,cx, cy, R, nmax):
"""comme la fonction temps, mais on considère la suite
avec z0= wx + i*wy au lieu de z0= 0. Donc
temps(cx, cy, R, nmax) == temps_generalise(0,0, cx, cy, R, nmax)
"""
n= 0
x= wx
y= wy
xx= x*x
yy= y*y
while xx + yy < R and n < nmax:
x, y = xx - yy + cx, 2*x*y + cy
xx= x*x
yy= y*y
n= n+1
return n
def julia(cx, cy, x, y, a, res= 100, R=10, nmax= 100):
# rappel (x,y) est le centre de l'image!!
haut= y+ a
gauche= x-a
bas= y-a # bas et droite ne sont pas utilisés, mais c'est pour
droite= x+a # vous aider à comprendre
dx= float( 2*a / res )
dy= dx # c'est plus clair d'avoir dy et dx je trouve, mais c'est inutile bien sûr
julia= np.zeros(shape= (res, res), dtype= np.int64)
for i in range(res):
for j in range(res):
julia[i,j]= temps_generalise(gauche+j*dx, haut-i*dy, cx, cy, R, nmax)
return julia
J= julia(-0.8, 0.156, 0, 0, 1,res= 600, nmax= 1000)
plt.imshow(J, cmap= 'gnuplot2')
Exercice (défi Mandelbrot) Utilisez Numpy habilement pour accélérer le calcul de l'ensemble de Mandelbrot (en évitant toutes les boucles for).
# I= le nombre complexe i, celui-là même dont le carré est -1
I= np.complex256(1j)
def mandelbrot_malin(centre, a, res= 1000, R= 10, nb_iterations= 100):
# on prépare un tableau numpy avec des nombres complexes,
# presque comme dans l'exo tab_complexe. La seule différence,
# c'est qu'on tient compte du fait que l'on regarde une portion
# carrée du plan centrée sur "centre", de côté 2*a, découpée en
# autant de pixels que la résolution l'exige.
all_x= np.linspace(centre.real-a, centre.real+a, num= res, dtype= np.complex256)
all_y= np.linspace(centre.imag+a, centre.imag-a, num= res, dtype= np.complex256)
X= np.stack([ all_x for _ in range(res) ])
Y= np.stack([ all_y for _ in range(res) ], axis= 1)
nombres_complexes= X+I*Y
# on aura besoin d'une copie qu'on ne touchera pas
valeurs_initiales= nombres_complexes.copy()
# un tableau numpy dans lequel on mettra l'image
image= np.zeros(shape= nombres_complexes.shape, dtype= np.int)
# un tableau numpy avec des True/False, qui nous dit quels sont les
# points pour lesquels le calcul n'est pas encore terminé
a_calculer= np.ones(shape= nombres_complexes.shape, dtype= np.bool)
# ici True= pas encore fini
# au début, on met tout sur "True". Et pour ça, on créé un tableau avec des 1
# en précisant que le type est "np.bool" ; Numpy a la convention que 1 = True
# dans cette situation.
# l'algorithme principal
for n in range(nb_iterations):
z= nombres_complexes[a_calculer] # tableau 1D extrait de "nombres_complexes"
c= valeurs_initiales[a_calculer] # ceci pour que c ait la même taille que z
# la formule principale pour calculer la suite est z_(n+1) = z_n^2 + c
# ce qui donne:
z= z**2 + c
# c'est précisément ici que numpy accélère énormément les choses!
pas_fini= np.abs(z) < R # tableau 1_D avec des True/False
# on met à jour:
nombres_complexes[a_calculer]= z # et hop on remet dans 'numbers', juste au bon endroit
# donc z_n a été remplacé par z_n+1, en gros
# la mise à jour la plus fine est sans doute:
a_calculer[a_calculer]= pas_fini
# la ligne précédente laisse les False contenus dans "a_calculer" sur False,
# et pour les True, ça les remplace par la valeur correspondante de "pas_fini"
# c'est pas tout ça, mais il faudrait penser à faire une image
masque= np.array(a_calculer, dtype= np.int) # convertit True en 1, False en 0
image= image + masque # on ajoute 1 si on a fait un autre tour de boucle...
return image
%time image= mandelbrot_malin(0 + 0.*I, 1.5, res= 4000)
plt.imshow(image, cmap= 'hot')
Commentaires sur la vitesse : pour une image d'une résolution de 1000x1000, on passe de 6.69 secondes avec la version précédente à 1.89 secondes ici (sur le même ordi), c'est donc 3.5 fois plus rapide. Pour une résolution de 4000x4000, on passe de 105 secondes à 28s, c'est 3.75 fois plus rapide.
En augmentant la résolution, l'écart grandit certainement, mais lentement.
J'avais promis une accélération bien plus spectaculaire que ça, mais ce n'est pas avec Numpy qu'on peut y arriver, en fait.