TP 6 : Numpy¶

Introduction¶

Le module numpy est l'un des plus utilisés. Son nom évoque le calcul numérique en Python, mais concrètement, il s'agit surtout de fonctionnalités liées aux tableaux.

Nous avons jusqu'ici traité les matrices comme des listes de listes. Mais il faut comprendre une chose : en Python, le fait même qu'une liste puisse contenir des éléments de n'importe quel type rend les choses très, très lentes. Lorsqu'on demande à Numpy de créer un tableau, il faut lui donner le type des éléments (par exemple "faire un tableau de 3 lignes et 4 colonnes composé de nombres entiers"), et sans rentrer dans les détails, ça change complètement la façon dont les données sont stockées en mémoire -- de manière bien plus compacte et efficace.

De plus, avec Numpy on va éviter de faire des boucles nous-mêmes, parce le module a été optimisé pour faire certaines choses très vite, comme appliquer une fonction à tous les éléments par exemple. (C'est notamment dû au fait que Numpy n'est pas lui-même écrit en Python, mais en langage C.)

Ici nous allons juste aborder les fonctionnalités les plus basiques de Numpy.

Création de tableaux¶

In [1]:
import numpy as np   # il est standard de l'appeler 'np'

Pour créer un tableau, on peut utiliser la fonction np.array :

In [2]:
A= np.array([[1,2], [3,4]])
In [3]:
A
Out[3]:
array([[1, 2],
       [3, 4]])

Nous avons laissé Numpy deviner le type des éléments. Pour vérifier :

In [4]:
A.dtype
Out[4]:
dtype('int64')

Ici int signifie "entier", et le 64 signifie "encodé sur 64 bits". En beaucoup plus clair (mais un informaticien ne vous le dira jamais comme ça), le type int64 est formé des éléments de $\mathbb{Z}/2^{64}\mathbb{Z}$. Le nombre $2^{64}$ étant énorme, on peut travailler avec ces entiers pour modéliser les éléments de $\mathbb{Z}$.

Pour choisir le type :

In [5]:
B= np.array([[1,2], [3,4]], dtype= np.int8)
print(B)

[[1 2]
 [3 4]]
In [6]:
C= np.array([[1,2], [3,4]], dtype= np.float64)  # nombre décimaux sur 64 bits
print(C)

[[1. 2.]
 [3. 4.]]

Pour faire des tableaux remplis de 0 ou de 1 :

In [7]:
np.zeros(shape= (4,9), dtype= np.int8)
Out[7]:
array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=int8)
In [8]:
np.ones(shape= (1,10), dtype= np.float64)
Out[8]:
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])

Dans ce dernier exemple, il y a une ligne et dix colonnes, mais c'est bel et bien un tableau (avec 2 dimensions). Or Numpy sait traiter n'importe quel nombre de dimensions ; voici des exemples :

In [9]:
np.array([1,2,3,4])
Out[9]:
array([1, 2, 3, 4])
In [10]:
np.ones(shape= (10,)) # notez le "10," ; un tuple avec un seul nombre, 10
Out[10]:
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
In [11]:
troisD= np.zeros(shape= (3,4,2))  # plus fort qu'une matrice : il y a trois dimensions
In [12]:
troisD  # affiche une liste de matrices
Out[12]:
array([[[0., 0.],
        [0., 0.],
        [0., 0.],
        [0., 0.]],

       [[0., 0.],
        [0., 0.],
        [0., 0.],
        [0., 0.]],

       [[0., 0.],
        [0., 0.],
        [0., 0.],
        [0., 0.]]])

Voyons comment accéder aux élements :

In [13]:
A= np.array([[1,2], [3,4]])
In [14]:
A[0,0]
Out[14]:
1
In [15]:
A[1,1]
Out[15]:
4
In [16]:
troisD[2,1,1]   # trois dimensions donc trois indices
Out[16]:
0.0

Et pour faire des modifications, sans surprise:

In [17]:
A[0,0]= 100
print(A)

[[100   2]
 [  3   4]]

Avant d'aborder l'exercice suivant, sachez que pour créer des nombres complexes en Python, on peut utiliser la syntaxe (un peu curieuse) I= 1j pour que I soit égal au complexe $i$, et ensuite faire x + y*I. Faites des essais, puis passez à l'exercice.

Exercice (Numpy et complexes) Cet exercice va paraître artificiel, mais il va servir plus loin pour créer des images intéressantes. On identifie $\mathbb{C}$ avec $\mathbb{R}^2$, et on va regarder le rectangle $[x_{min}, x_{max}] \times [y_{min}, y_{max}]$. On se donne un entier $N$ (la "résolution"), et on va prendre des valeurs dans le rectangle régulièrement espacées. Horizontalement, entre deux valeurs consécutives il y aura $dx= (x_{max} - x_{min})/N$ et verticalement il y aura $dy= (y_{max} - y_{min})/N$.

Bref, on vous demande d'écrire une fonction tab_complexe(xmin, xmax ymin, ymax, N) qui renvoie un tableau Numpy de taille (N+1)x(N+1), dont les éléments sont de type np.complex64, et tel que si A= tab_complexe(xmin, xmax, ymin, ymax, N), alors A[i,j] est (xmin + j * dx) + I*(ymax - i*dy) avec $dx$ et $dy$ comme ci-dessus, et $I= 1j$.

Note: si la formule vous paraît bizarre (ce qui n'empêche pas de faire l'exercice), rappelez vous que si $A$ est une matrice et qu'on parle de $A[i,j]$, alors $i$ est la ligne et $j$ est la colonne, et en plus on place la ligne $i+1$ en dessous de la ligne $i$ ; dans le plan, quand on parle du point $(i,j)$, c'est tout-à-fait le contraire.

Par exemple tab_complexe(0,4,0,4,4) doit renvoyer

array([[0.+4.j, 1.+4.j, 2.+4.j, 3.+4.j, 4.+4.j],
       [0.+3.j, 1.+3.j, 2.+3.j, 3.+3.j, 4.+3.j],
       [0.+2.j, 1.+2.j, 2.+2.j, 3.+2.j, 4.+2.j],
       [0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j, 4.+1.j],
       [0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j]], dtype=complex64)

Notez que l'on a bien $4i$ en haut à gauche, $4+4i$ en haut à droite, $0$ en bas à gauche et $4$ en bas à droite, comme on le voudrait.

In [18]:
def tab_complexe(xmin, xmax, ymin, ymax, N):
    A= np.zeros(shape= (N+1,N+1), dtype= np.complex64)
    I= 1j
    dx= (xmax-xmin)/N
    dy= (ymax-ymin)/N
    for i in range(N+1):
        for j in range(N+1):
            A[i,j]= (xmin + j * dx) + I*(ymax - i*dy)
    return A
In [19]:
tab_complexe(0,4,0,4,4)
Out[19]:
array([[0.+4.j, 1.+4.j, 2.+4.j, 3.+4.j, 4.+4.j],
       [0.+3.j, 1.+3.j, 2.+3.j, 3.+3.j, 4.+3.j],
       [0.+2.j, 1.+2.j, 2.+2.j, 3.+2.j, 4.+2.j],
       [0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j, 4.+1.j],
       [0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j]], dtype=complex64)

Voici la correction de la version rapide (cf exercice suivant):

In [20]:
def tab_rapide(xmin, xmax, ymin, ymax, N):
    L= np.linspace(xmin, xmax, N+1, dtype= np.complex64)
    LL= [L for i in range(N+1)]
    X= np.stack(LL)
    
    L= np.linspace(ymax, ymin, N+1, dtype= np.complex64)
    LL= [L for i in range(N+1)]
    Y= np.stack(LL, axis= 1)
    
    return X + 1j*Y
In [21]:
tab_rapide(0,4,0,4,4)
Out[21]:
array([[0.+4.j, 1.+4.j, 2.+4.j, 3.+4.j, 4.+4.j],
       [0.+3.j, 1.+3.j, 2.+3.j, 3.+3.j, 4.+3.j],
       [0.+2.j, 1.+2.j, 2.+2.j, 3.+2.j, 4.+2.j],
       [0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j, 4.+1.j],
       [0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j]], dtype=complex64)
In [22]:
%time T= tab_complexe(0,1,0,1,10000)

CPU times: user 26.6 s, sys: 279 ms, total: 26.9 s
Wall time: 26.9 s
In [23]:
%time T= tab_rapide(0,1,0,1,10000)

CPU times: user 800 ms, sys: 803 ms, total: 1.6 s
Wall time: 1.71 s

Opérations¶

On peut demander à Numpy d'appliquer une fonction à tous les éléments d'un tableau. Dans la plupart des cas, il s'agit de fonctions particulières, qui font partie du module Numpy. Par exemple, pour calculer le sinus de tous les nombres dans un tableau, on utilisera np.sin.

In [24]:
A= np.array([[1,2], [3,4]])
In [25]:
A
Out[25]:
array([[1, 2],
       [3, 4]])
In [26]:
np.sin(A)
Out[26]:
array([[ 0.84147098,  0.90929743],
       [ 0.14112001, -0.7568025 ]])

Une chose très agréable, c'est que les objets Numpy ont redéfini les opérations arithmétiques telles que + et * (nous verrons comment faire ça dans le dernier TP!). On peut donc écrire :

In [27]:
A+A
Out[27]:
array([[2, 4],
       [6, 8]])
In [28]:
A*A
Out[28]:
array([[ 1,  4],
       [ 9, 16]])

Attention, il ne s'agit pas du produit matriciel, mais du produit terme à terme. Il faut utiliser np.matmul(A,B) pour le produit des deux matrices $A$ et $B$ (celui qui n'est défini que quand $A$ possède autant de colonnes que $B$ a de lignes) ; nous n'en parlerons plus dans ce TP.

Numpy comprend aussi :

In [29]:
A-1
Out[29]:
array([[0, 1],
       [2, 3]])

De sorte que si vous avez une fonction telle que :

In [30]:
def f(x):
    return x**2 - 1

Alors on peut faire :

In [31]:
f(A)
Out[31]:
array([[ 0,  3],
       [ 8, 15]])

Au total, la fonction f a été appliquée aux éléments de A, un par un. Mais ça ne marche que parce que la définition de f n'utilise que des opérations que Numpy sait interpréter!

Remarque. Il y a une syntaxe pour appliquer une fonction quelconque à tous les éléments d'un tableau, à savoir np.vectorize(f)(A), mais ça n'a pas un grand intérêt car c'est très lent, et il reviendrait au même de faire une boucle.

Voyons un peu si on y gagne en vitesse. Prenons une matrice au hasard:

In [32]:
np.random.randint(10, size=(3,9))  # pour comprendre la syntaxe
Out[32]:
array([[7, 0, 4, 6, 0, 6, 1, 2, 5],
       [6, 2, 2, 3, 2, 5, 0, 0, 0],
       [7, 8, 7, 4, 4, 1, 8, 5, 8]])

Mais prenons-la beaucoup plus grosse :

In [33]:
N= 1000
A= np.random.randint(300, size=(N,N))
In [34]:
# faisons une copie de A, mais comme une liste de listes:
B= [ [A[i,j] for j in range(N)] for i in range(N)]
In [35]:
# une fonction qui met tout au carré dans B :
def test():
    for i in range(N):
        for j in range(N):
            B[i][j]= B[i][j]**2
In [36]:
%time test()  # essayons

CPU times: user 294 ms, sys: 4 ms, total: 298 ms
Wall time: 296 ms
In [37]:
%time A= A**2 # maintenant comparons avec la fonction "puissance" de Numpy

CPU times: user 1.39 ms, sys: 1.91 ms, total: 3.31 ms
Wall time: 2.05 ms

Avec la machine sur laquelle ce test a été fait la première fois, la deuxième méthode est 100 fois plus rapide. Et ça augmente avec la taille des tableaux.

Cette amélioration impressionnante explique pourquoi, dans certains cours de calcul numérique avec Python, on va jusqu'à dire que "les bons programmeurs ne font jamais de boucle for". Ils utilisent Numpy. C'est un peu exagéré.

Exercice (comme le précédent en plus rapide) Essayez les commandes suivantes, sur des exemples ou bien en regardant l'aide: np.stack(LL) où LL est une liste de listes, np.stack(LL, axis= 1), et puis

np.linspace(xmin, xmax, N, dtype= np.complex64).

Puis, en combinant ces commandes, écrire une fonction tab_rapide(xmin, xmax, ymin, ymax, N) qui renvoie rigoureusement la même chose que tab_complexe(xmin, xmax, ymin, ymax, N), mais sans faire de boucle for.

Comparez la vitesse (c'est normalement 10 fois plus rapide, au bas mot).

(Attention, la correction est ci-dessus à côté de l'exercice précédent.)

Numpy et les booléens¶

Voyons maintenant comment Numpy gère les "booléens" -- les conditions.

In [38]:
A= np.array([[1,2], [3,4]])
print(A == 2)   # renvoie un tableau de booléens

[[False  True]
 [False False]]
In [39]:
np.logical_or(A == 2, A == 3)  # voilà comment on fait un "ou"
Out[39]:
array([[False,  True],
       [ True, False]])
In [40]:
np.logical_and(1<A, A<4)    # et voilà le "et"
Out[40]:
array([[False,  True],
       [ True, False]])

On peut récupérer les éléments qui satisfont une condition. Par exemple si:

In [41]:
condition= A > 1
print(condition)

[[False  True]
 [ True  True]]

Alors on peut faire:

In [42]:
A[condition]
Out[42]:
array([2, 3, 4])

La réponse est de dimension 1 (on ne peut pas prévoir à l'avance combien d'éléments vont satisfaire la condition, et comment ils seront disposés!). Mais les deux syntaxes suivantes fonctionnent aussi, et c'est très pratique :

In [43]:
A[condition]= 0  # on donne une seule valeur
print(A)

[[1 0]
 [0 0]]

La réponse est bien encore un tableau de dimension 2! Prenez le temps de méditer sur le fait que les auteurs de Numpy ont redéfini les crochets et le signe =, entre autres.

La deuxième façon est :

In [44]:
A[condition]= np.array([11, 12, 13])  # on donne une liste
print(A)

[[ 1 11]
 [12 13]]

Exercice (booléens Numpy 1) Ecrire une fonction f(A) qui prend en paramètre un tableau Numpy A (n'importe quelle taille, n'importe quelle nombre de dimensions) et qui renvoie un tableau B qui est une copie de A dans laquelle les nombres $<0$ ont été remplacés par $0$ et les nombres $>1$ ont été remplacés par $1$. Vous pouvez faire A.copy() pour obtenir une copie de A.

Exercice (booléens Numpy 2) Si A est un tableau Numpy, expliquez ce que donne la commande

A[ A<0 ]= - A[ A<0 ]

et pourquoi.

In [45]:
def f(A):
    B= A.copy()
    B[ A < 0 ]= 0
    B[ A > 1 ]= 1
    return B
In [46]:
A= np.array([[-1, 2], [-3, 4]])
In [47]:
A
Out[47]:
array([[-1,  2],
       [-3,  4]])
In [48]:
f(A)
Out[48]:
array([[0, 1],
       [0, 1]])
In [49]:
A[ A < 0 ]= - A[ A < 0 ]
In [50]:
A
Out[50]:
array([[1, 2],
       [3, 4]])

Images à partir de matrices¶

Nous allons produire quelques images. Pour cela, nous allons faire appel au module matplotlib, que nous étudierons plus en détails dans le TP suivant.

Pour commencer, il faut importer l'objet suivant -- il est standard de l'appeler plt:

In [52]:
from matplotlib import pyplot as plt

La commande suivante permet d'agrandir un peu les figures, par défaut:

In [53]:
plt.rcParams["figure.figsize"] = (10,10)

Supposons que vous ayez une matrice $A$. Matplotlib possède une syntaxe toute simple pour faire la chose suivante, qui revient étonnamment souvent : on interprète les nombres dans $A$ comme des couleurs, et on affiche directement ça comme une image, chaque entrée de la matrice contrôlant un pixel. Pour interpréter les nombres comme des couleurs, matplotlib utilise des "colormaps" ; ce sont des listes de couleurs, et la première sera attribuée à la valeur minimale dans la matrice, la dernière à la valeur maximale, le reste étant interpolé.

Un exemple va rendre ça beaucoup plus clair!

In [54]:
N= 15
A= np.array([ [i]*N for i in range(N)   ])
print(A)

[[ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]
 [ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1]
 [ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2]
 [ 3  3  3  3  3  3  3  3  3  3  3  3  3  3  3]
 [ 4  4  4  4  4  4  4  4  4  4  4  4  4  4  4]
 [ 5  5  5  5  5  5  5  5  5  5  5  5  5  5  5]
 [ 6  6  6  6  6  6  6  6  6  6  6  6  6  6  6]
 [ 7  7  7  7  7  7  7  7  7  7  7  7  7  7  7]
 [ 8  8  8  8  8  8  8  8  8  8  8  8  8  8  8]
 [ 9  9  9  9  9  9  9  9  9  9  9  9  9  9  9]
 [10 10 10 10 10 10 10 10 10 10 10 10 10 10 10]
 [11 11 11 11 11 11 11 11 11 11 11 11 11 11 11]
 [12 12 12 12 12 12 12 12 12 12 12 12 12 12 12]
 [13 13 13 13 13 13 13 13 13 13 13 13 13 13 13]
 [14 14 14 14 14 14 14 14 14 14 14 14 14 14 14]]

La commande en question s'appelle imshow :

In [55]:
plt.imshow(A)
Out[55]:
<matplotlib.image.AxesImage at 0x7fd260ff8070>

Ici la "colormap" est celle par défaut (elle va du bleu au jaune). Vous pouvez augmenter $N$ ci-dessus, c'est plus joli. Essayons autre chose :

In [56]:
plt.imshow(A, cmap= "hot")
Out[56]:
<matplotlib.image.AxesImage at 0x7fd291e55b20>

La "colormap" qui s'appelle "hot" va du noir au blanc en passant par les couleurs des flammes, de la plus froide à la plus chaude.

Pour sauvegarder ça comme une image, c'est un peu compliqué, il faut importer un sous-module:

In [57]:
import matplotlib.image as mpl

Et maintenant :

In [58]:
mpl.imsave("mon-image.png", A, cmap= "hot")

Pour une liste de colormaps, allez voir :

https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html

Vous pouvez ignorer le texte!

Une chose bien utile, que l'on peut mettre en valeur ici, c'est qu'en ajoutant _r au nom d'une colormap, on l'obtient inversée (reverse). Exemple :

In [59]:
plt.imshow(A, cmap= "hot_r")
Out[59]:
<matplotlib.image.AxesImage at 0x7fd2724478b0>

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 d'un TP précédent.) Attention, pour $k > n$ on veut que la fonction renvoie $0$.

Puis, faire un tableau Numpy A carré, de taille $2^a\times 2^a$ où $a$ est un entier que vous choisissez, tel que A[n,k] = binomial(n,k) % 2. Puis, affichez A comme une image. Augmentez $a$.

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).

In [60]:
def factorielle(n):
    if n == 0:
        return 1
    else:
        return n*factorielle(n-1)
In [62]:
def binomial(n,k):  # on pourrait faire bien plus efficace que ça...
    if k > n:
        return 0
    else:
        return factorielle(n) // (factorielle(k)*factorielle(n-k))
In [115]:
N= 2**7
pascal= [ [binomial(n,k)%2 for k in range(N)] for n in range(N) ]
image= np.array(pascal)
In [116]:
plt.imshow(image, cmap= "hot")
Out[116]:
<matplotlib.image.AxesImage at 0x7fd292d6c7c0>

Note: essayez d'augmenter $N$. Pour des valeurs plus grandes, on voit plusieurs couleurs apparaître, c'est un bug. Je crois que quand l'image à afficher est plus grande que la place qu'on lui donne à l'écran, matplotlib prend tout seul des décisions étranges. Preuve que c'est un bug: si vous sauvegardez l'image dans un fichier (voir instructions ci-dessus) et que vous l'ouvrez séparément, c'est bel et bien en noir et blanc, ouf.

L'ensemble de Mandelbrot¶

Voici un grand classique des cours de programmation. On va dessiner (une approximation de) l'ensemble de Mandelbrot. Il s'agit de

$$ M= \{ c \in \mathbb{C} ~~|~~ \text{la suite définie par}~ z_0 = 0, z_{n+1} = z_n^2 + c ~\text{reste bornée} \}.$$

Pour un $c$ donné, il est exclu de vérifier vraiment si la suite en question reste bornée, et en pratique, une façon de faire qui produit des résultats corrects est la suivante (sachant que l'on veut utiliser Numpy pour accélérer les calculs). On prend un entier $N$ "assez grand", disons $N=100$, et pour chaque $c$ que l'on veut tester, on calcule les termes de la suite ci-dessus jusqu'à $z_N$. Si $|z_N] > 10$, la suite va tendre vers $+\infty$ et $c \not\in M$; dans le cas contraire, on suppose que $N$ est assez grand pour en conclure que la suite reste bornée et que $c\in M$.

On peut placer les résultats dans une matrice, avec un 1 si $c\in M$ et un 0 sinon, puis afficher la matrice comme dans les exemples ci-dessus.

Ecrivons déjà une fonction pour calculer les termes de la suite -- notez bien que le code ci-dessous marche aussi bien pour des nombres complexes que pour des tableaux Numpy!

In [117]:
def iterations(nmax,c):
    """On regarde la suite z0= 0, z_{n+1} = z_n^2 + c, et on calcule son terme z_nmax."""
    z= 0*c  # astuce pour que ça soit un tableau si c en est un!
    for i in range(nmax):
        z= z**2 + c
    return z

Si vous essayez cette fonction avec par exemple $c= 0.5$ et $N=100$, vous allez voir que Python se plaint que les nombres deviennent trop grands. Heureusement, dans un tableau Numpy, ça ne produit qu'un "warning" et non pas une erreur, et on peut même désactiver les "warnings" comme ceci:

In [118]:
import warnings
warnings.filterwarnings("ignore")

(Ce n'est pas très élégant. Ne répétez pas que vous avez appris ça en cours.)

In [119]:
c= np.array([0.5])
In [120]:
iterations(100, c)
Out[120]:
array([inf])

Voilà, le résultat fait $+\infty$.

Tout ce qu'il nous reste à faire, c'est prendre un tableau de nombres exactement comme tab_rapide en produit, appliquer la fonction dessus, et voir quelles nombres dans le tableau sont $>10$ en valeur absolue.

On commence par le carré $[- 1.5, 0.5] \times [-1, 1]$:

In [121]:
C= tab_rapide(-1.5,0.5, -1,1, 1000)  # image 1000 x 1000 pixels
In [122]:
resultat= iterations(100, C)
In [123]:
dessin= np.abs(resultat) < 10 # tableau avec des False et des True
In [124]:
plt.imshow(dessin, cmap= "Blues")  # matplotlib va interpréter False comme 0, True comme 1
Out[124]:
<matplotlib.image.AxesImage at 0x7fd1e80523d0>

Cet ensemble est très riche (vous irez chez vous sur youtube pour voir des images incroyables). Voici la "vallée des hippocampes":

In [125]:
a= .05
C= tab_rapide(-.75-a,-0.75+a, 0.1-a,0.1+a, 1000)
resultat= iterations(100, C)
dessin= np.abs(resultat) < 10 
plt.imshow(dessin, cmap= "Blues")
Out[125]:
<matplotlib.image.AxesImage at 0x7fd262cb81f0>

Ces dessins s'affichent en quelques secondes seulement, et on peut remercier Numpy pour ça.

Exercice (Mandelbrot en couleurs, plus difficile) Cette fois ci on va remplacer les appels à Numpy par des boucles for, ce qui va ralentir énormément les choses, mais va permettre de produire des images en couleurs très jolies.

(1) Ecrire une fonction temps(nmax,c) qui commence à itérer la suite comme iterations mais s'arrête lorsque $|z_n|>10$, si toutefois ça arrive pour un $n < n_{max}$ ; et la fonction renvoie l'entier $n$ (ou $n_{max}$ en cas d'échec après $n_{max}$ tentatives), ce qui est une mesure du "temps qu'il faut pour que la suite donne un terme de module $>10$".

(2) Ecrire une fonction mandelbrot(xmin, xmax, ymin, ymax, N, nmax) qui crée un tableau de taille N x N qui contient les diverses valeurs de temps(nmax, c) pour $c$ dans le rectangle $[x_{min}, x_{max}] \times [y_{min}, y_{max}]$.

In [126]:
def temps(cx, cy, R, nmax):
    """nombre d'itérations nécessaires pour que la suite, pour la valeur
    c= cx + cy*i, donne un terme de module > R. On fait nmax tentatives, 
    et on renvoie nmax si ça n'a pas marché.    
    """
    n= 0
    x= 0
    y= 0
    xx= 0
    yy= 0
    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
In [127]:
def mandelbrot(x, y, a, res= 100, R=10, nmax= 100):
    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

    mandel= np.zeros(shape= (res, res), dtype= np.int64)
    for i in range(res):
        for j in range(res):
            mandel[i,j]= temps(gauche+j*dx, haut-i*dy, R, nmax)

    return mandel
In [128]:
a= 0.05
x= -0.75
y= 0.1
In [129]:
%time M= mandelbrot(x, y, a, res= 1000)

CPU times: user 17.9 s, sys: 28.4 ms, total: 17.9 s
Wall time: 17.9 s
In [130]:
plt.imshow(M, cmap= 'Spectral')
Out[130]:
<matplotlib.image.AxesImage at 0x7fd26292ee80>

Exercice (ensemble de Julia, plus difficile) 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.

In [131]:
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
In [132]:
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
In [133]:
J= julia(-0.8, 0.156, 0, 0, 1,res= 600, nmax= 1000)
In [134]:
plt.imshow(J, cmap= 'gnuplot2')
Out[134]:
<matplotlib.image.AxesImage at 0x7fd262cf4b80>