TP 4 : Modules mathématiques, partie 1 : Numpy et Matplotlib

1 - Utiliser des modules

Nous avons déjà vu, dans le premier TP, comment utiliser une fonction contenue dans un module, c'est-à-dire une collection de programmes écrits par quelqu'un d'autre. Mais nous allons explorer ça un peu plus avec le module math, qui donne accès à des fonctionnalités type calculatrice (sinus et cosinus, etc).

Rappelons comment on importe une fonction isolée :

In [1]:
from math import sqrt # sqrt pour 'square root', racine carrée
In [2]:
sqrt(2)
Out[2]:
1.4142135623730951

On peut aussi changer le nom si on préfère :

In [3]:
from math import sqrt as racine
In [4]:
racine(2)
Out[4]:
1.4142135623730951

Il est possible d'utiliser from math import * pour tout importer, mais il n'est vraiment pas recommandé de faire ça. (Il y a bien des raisons, mais imaginez déjà que le module math contienne une fonction qui s'appelle exactement comme l'une des vôtres, sans que vous le sachiez...)

On préfèrera importer le module lui-même en faisant :

In [5]:
import math
In [6]:
math.sqrt(2)
Out[6]:
1.4142135623730951

Pour voir toutes les fonctions du module math, maintenant qu'il est importé, tapez math. avec le point et utilisez la touche TAB.

In [7]:
math.pi
Out[7]:
3.141592653589793

Pour les plus paresseux, on peut également abréger le nom du module :

In [8]:
import math as m
In [9]:
m.sin( m.pi )
Out[9]:
1.2246467991473532e-16

L'écriture ci-dessus indique $1.22 \cdot 10^{-16}$, qui est vraiment très proche de 0. Le nombre pi tel que le définit le module math est un "float", et la fonction sinus de ce module renvoie des floats également, qui sont des nombres approchés.

2 - Numpy

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 [10]:
import numpy as np   # il est standard de l'appeler 'np'

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

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

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

In [13]:
A.dtype  # A est un objet, évidemment !
Out[13]:
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 nombres pour modéliser les éléments de $\mathbb{Z}$.

Pour choisir le type :

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

[[1 2]
 [3 4]]
In [15]:
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 [16]:
np.zeros(shape= (4,9), dtype= np.int8)
Out[16]:
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 [17]:
np.ones(shape= (1,10), dtype= np.float)   # float sur le nombre de bits par défaut
Out[17]:
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 [18]:
np.array([1,2,3,4])
Out[18]:
array([1, 2, 3, 4])
In [19]:
np.ones(shape= (10,)) # notez le "10," ; un tuple avec un seul nombre, 10
Out[19]:
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
In [20]:
troisD= np.zeros(shape= (3,4,2))  # plus fort qu'une matrice : il y a trois dimensions
In [21]:
troisD  # affiche une liste de matrices
Out[21]:
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 [22]:
A= np.array([[1,2], [3,4]])
In [23]:
A[0,0]
Out[23]:
1
In [24]:
A[1,1]
Out[24]:
4
In [25]:
troisD[2,1,1]   # trois dimensions donc trois indices
Out[25]:
0.0

Et pour faire des modifications, sans surprise:

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

[[100   2]
 [  3   4]]

Exercice (Numpy et complexes) Ecrire une fonction tab_complexe(N) qui renvoie un tableau Numpy de taille NxN, dont les éléments sont de type np.complex, et tel que si A= tab_complexe(N) alors A[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= 1j pour que I soit égal au complexe $i$, et ensuite faire x + 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).

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 [27]:
A= np.array([[1,2], [3,4]])
In [28]:
A
Out[28]:
array([[1, 2],
       [3, 4]])
In [29]:
np.sin(A)
Out[29]:
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 * (comme nous dans le dernier TP). On peut donc écrire :

In [30]:
A+A
Out[30]:
array([[2, 4],
       [6, 8]])
In [31]:
A*A
Out[31]:
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 [32]:
A-1
Out[32]:
array([[0, 1],
       [2, 3]])

De sorte que si vous avez une fonction telle que :

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

Alors on peut faire :

In [34]:
f(A)
Out[34]:
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 [35]:
np.random.randint(10, size=(3,9))  # pour comprendre la syntaxe
Out[35]:
array([[6, 4, 7, 5, 5, 2, 3, 4, 5],
       [0, 2, 3, 6, 3, 2, 9, 4, 2],
       [0, 6, 0, 5, 8, 2, 9, 7, 2]])

Mais prenons-la beaucoup plus grosse :

In [36]:
N= 1000
A= np.random.randint(300, size=(N,N))
In [37]:
# 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 [38]:
# 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 [39]:
%time test()  # essayons

CPU times: user 326 ms, sys: 2.82 ms, total: 329 ms
Wall time: 332 ms
In [40]:
%time A= A**2 # maintenant comparons avec la fonction "puissance" de Numpy

CPU times: user 3.48 ms, sys: 5.13 ms, total: 8.62 ms
Wall time: 9.29 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é.

Numpy et les booléens

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

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

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

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

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

[[False  True]
 [ True  True]]

Alors on peut faire:

In [45]:
A[condition]
Out[45]:
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 [46]:
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! Ici on voit que les auteurs de Numpy ont été habiles quand ils ont redéfini les crochets et le signe = (ils ont même dû utiliser des choses que l'on n'a pas décrites dans le TP précédent).

La deuxième façon est :

In [47]:
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.

3 - Matplotlib

Il s'agit d'un module pour faire des graphiques, au sens large. Le nombre de possibilités est énorme, mais nous ne verrons que les choses les plus simples.

On notera que matplotlib s'appuie énormément sur numpy.

Syntaxe de base

Commençons par importer l'objet pyplot, en le renommant plt, ce qui est standard

In [48]:
from matplotlib import pyplot as plt

Voici une syntaxe de base : la méthode plt.plot(x,y), où x et y sont des listes (ou tableaux, ou tout ce qu'on veut) de nombres, de même longueur ; ceci trace une ligne brisée dans le plan, passant par les points correspondants. Un exemple vaudra mieux :

In [49]:
x = [5, 2, 9, 4, 7]  
y = [10, 5, 8, 4, 2]  
# ligne brisée passant par (5,10) puis (2,5) etc :
plt.plot(x, y)  
#plt.show()
Out[49]:
[<matplotlib.lines.Line2D at 0x109c0d4a8>]

Il y a une remarque importante à faire tout de suite. Tout d'abord, la ligne plt.show() n'est pas obligatoire (essayez de l'enlever!) : c'est un peu comme de finir une cellule par un print, c'est souvent inutile parce qu'on sait bien que le notebook jupyter essaie d'afficher le résultat de la dernière ligne ; ici de toute façon, peu importe la dernière ligne, si vous avez touché à plt dans la cellule, il va y avoir un appel à show.

Pourquoi est-ce si important? Parce qu'à chaque fois qu'on appelle show, la figure est effacée, et l'objet plt est prêt à partir sur une nouvelle figure. D'ailleurs, essayons :

In [50]:
plt.show()

Cette commande ne produit aucun effet : il n'y a plus de figure "en cours" -- ou plutôt si, il y en a une, qui est vide.

Moralité : quand on fait une figure, on la fait dans une seule cellule, et après elle est "perdue". (Matplotlib offre la possibilité de faire autrement, mais on ne va pas s'embêter avec les détails de ça.)

Et ajoutons tout de suite qu'on peut sauvegarder la figure, par exemple en PDF :

In [51]:
# il faut la refaire, du coup!
x = [5, 2, 9, 4, 7]  
y = [10, 5, 8, 4, 2]  
plt.plot(x, y)  

plt.savefig("ma-figure.pdf")

Vous voyez bien que la figure est affichée à l'écran, alors qu'on n'a rien demandé. (Au fait, en cliquant à gauche de la figure sur la bande bleue verticale qui apparaît, vous pouvez la cacher.) Mais surtout, si vous êtes dans jupyterlab, vous verrez à gauche, dans quelques secondes, le fichier ma-figure.pdf apparaître. Double-cliquez dessus, vous verrez le résultat. Essayez de zoomer.

Ajoutons ici la syntaxe pour que les figures soient plus grosses dans le notebook. Ça n'a pas d'incidence sur les figures sauvegardées.

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

Maintenant, exécutez à nouveau les cellules précédentes : les figures sont plus grosses, jusqu'à nouvel ordre.

Graphes de fonctions

Voyons maintenant comment tracer le graphe d'une fonction simple $f$ définie sur un intervalle de $\mathbb{R}$. On le fait simplement avec une ligne brisée qui relie les points $(x, f(x))$ pour un grand nombre de valeurs de $x$. Pour ça, utilisons Numpy :

In [53]:
x = np.linspace(-1, 2, 100)
print(x)

[-1.         -0.96969697 -0.93939394 -0.90909091 -0.87878788 -0.84848485
 -0.81818182 -0.78787879 -0.75757576 -0.72727273 -0.6969697  -0.66666667
 -0.63636364 -0.60606061 -0.57575758 -0.54545455 -0.51515152 -0.48484848
 -0.45454545 -0.42424242 -0.39393939 -0.36363636 -0.33333333 -0.3030303
 -0.27272727 -0.24242424 -0.21212121 -0.18181818 -0.15151515 -0.12121212
 -0.09090909 -0.06060606 -0.03030303  0.          0.03030303  0.06060606
  0.09090909  0.12121212  0.15151515  0.18181818  0.21212121  0.24242424
  0.27272727  0.3030303   0.33333333  0.36363636  0.39393939  0.42424242
  0.45454545  0.48484848  0.51515152  0.54545455  0.57575758  0.60606061
  0.63636364  0.66666667  0.6969697   0.72727273  0.75757576  0.78787879
  0.81818182  0.84848485  0.87878788  0.90909091  0.93939394  0.96969697
  1.          1.03030303  1.06060606  1.09090909  1.12121212  1.15151515
  1.18181818  1.21212121  1.24242424  1.27272727  1.3030303   1.33333333
  1.36363636  1.39393939  1.42424242  1.45454545  1.48484848  1.51515152
  1.54545455  1.57575758  1.60606061  1.63636364  1.66666667  1.6969697
  1.72727273  1.75757576  1.78787879  1.81818182  1.84848485  1.87878788
  1.90909091  1.93939394  1.96969697  2.        ]

La commande ci-dessus crée une liste (qui est tableau Numpy, en fait) comprenant 100 valeurs régulièrement espacées dans l'intervalle $[-1, 2]$. En combinant Numpy et Matplotlib, on peut maintenant faire :

In [54]:
plt.plot(x, x**2 + 1)
Out[54]:
[<matplotlib.lines.Line2D at 0x109ed3668>]

Et voilà le graphe de $x\mapsto x^2+1$, tout simplement. Pour faire celui, disons, de la fonction sinus sur $[0, 2\pi]$, ne pas oublier que c'est la fonction sinus de Numpy :

In [55]:
x= np.linspace(0, 2*np.pi, 100)
plt.plot(x, np.sin(x))
Out[55]:
[<matplotlib.lines.Line2D at 0x10a133828>]

Voici un exemple plus complet, sans commentaires, qui illustre quelques unes des options disponibles.

In [56]:
x = np.linspace(0, 2, 100)

plt.plot(x, x, label='linéaire')  
plt.plot(x, x**2, label='quadratique')  
plt.plot(x, x**3, label='cubique')  
plt.xlabel('abscisses')  
plt.ylabel('ordonnées')  
plt.title("Un exemple")  
plt.legend()  

plt.show()

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.plot pour 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$.

Images à partir de matrices

Terminons ce petit (mini) tour d'horizon par le traitement des tableaux comme des images. 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 [57]:
N= 256
A= np.array([ [i]*N for i in range(N)   ])
print(A)

[[  0   0   0 ...   0   0   0]
 [  1   1   1 ...   1   1   1]
 [  2   2   2 ...   2   2   2]
 ...
 [253 253 253 ... 253 253 253]
 [254 254 254 ... 254 254 254]
 [255 255 255 ... 255 255 255]]

La commande en question s'appelle imshow :

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

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 [59]:
plt.imshow(A, cmap= "hot")
Out[59]:
<matplotlib.image.AxesImage at 0x10aa9aa20>

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 [60]:
import matplotlib.image as mpl

Et maintenant :

In [61]:
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 [62]:
plt.imshow(A, cmap= "hot_r")
Out[62]:
<matplotlib.image.AxesImage at 0x10aafabe0>

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

4 - 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 de bons résultats est la suivante. On se donne un nombre $R>0$ et un entier $N$. Ensuite pour un $c$ donné, on regarde si on trouve un $n$ tel que $|z_n| > R$, où $(z_n)$ est la suite comme ci-dessus ; si on n'en trouve aucun après avoir essayé tous les $n$ avec $n < N$, on laisse tomber. On note $n_c$ le plus petit entier $n$ que l'on a trouvé, avec $n_c=N$ si on n'en a pas trouvé. On se retrouve avec une collection d'entiers $n_c$ indexés par les nombres complexes $c$ que l'on a essayés, ce qui donne une matrice (lignes et colonnes donnant les parties imaginaires et réelles). Finalement on affiche cette matrice avec une colormap bien sentie.

Il est assez raisonnable de penser que, si on n'a pas dépassé le rayon $R$ après $N$ itérations, l'entier $N$ étant assez (!) grand, alors la suite est probablement bornée. Un peu plus mathématiquement (mais sans être très rigoureux non plus), on peut commencer par montrer que $M$ est borné ; ensuite, on en déduit qu'en choisissant $R$ assez grand, alors la condition $|z_n| > R$ entraîne que la suite en question va tendre vers $\infty$ en module, et donc le $c$ qu'on regarde n'est effectivement pas dans $M$. Donc avec l'algorithme ci-dessus, tous les $c$ pour lesquels on trouve un $n_c$ qui n'est pas $N$ ne sont pas dans $M$.

In [63]:
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

Quelques commentaires sur cette première fonction. On n'utilise pas de nombres complexes ici, on travaille avec $x$ et $y$, les parties réelle et imaginaire de $z= z_n = x +iy$. Ensuite, on a utilisé la formule $$(x+iy)^2 = x^2 - y^2 + 2xyi.$$

Notez que l'on a bien fait attention à ne pas calculer $x^2$ et $y^2$ deux fois, c'est pourquoi on les stocke dans xx et yy.

Enfin, on a utilisé la syntaxe Python que vous ne connaissez pas encore :

a, b = val1, val2

Ça permet de mettre val1 dans x et val2 dans y d'un seul coup ; mais ça n'est pas équivalent à :

a= val1

b= val2

parce que val2 peut être une expression qui dépend de a ! Dans la fonction ci-dessus, vous avez un bel exemple. Prenez le temps de bien comprendre.

Traçons maintenant notre ensemble de Mandelbrot approximatif. On va se donner un $a>0$ et deux nombres $x, y$ qui donnent le centre du carré qui va de $(x-a, y-a)$ à $(x+a, y+a)$. (Comme ça, en gardant $x$ et $y$ et en diminuant $a$, on pourra zoomer.)

On parcourt le carré pixel par pixel, on appelle la fonction temps à chaque fois, qui nous donne le "temps nécessaire pour dépasser $R$", et on stocke tout ça dans un tableau. L'entier res (comme résolution) donne la taille du tableau.

Notez bien que, pour la première fois, on définit une fonction avec des valeurs par défaut pour les paramètres : par exemple, dans l'exemple qui suit on ne donnera ni R ni nmax, et alors Python prendra 10 et 100.

In [64]:
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 [65]:
x= 0
y= 0
a= 1.5
In [66]:
%time M= mandelbrot(x, y, a, res= 400)

CPU times: user 1.01 s, sys: 7.75 ms, total: 1.01 s
Wall time: 1.04 s
In [67]:
plt.imshow(M, cmap= 'hot')
Out[67]:
<matplotlib.image.AxesImage at 0x10acc7f98>

Ici, tous les pixels où la couleur n'est pas le blanc correspondent à des points qui ne sont pas dans $M$. Pour les points en blanc, on peut seulement supposer qu'ils sont "probablement" dans $M$.

Amusez vous à zoomer et à changer de colormap. Envoyez vos plus belles images à la classe.

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.

Exercice (défi Mandelbrot) Utilisez Numpy habilement pour accélérer le calcul de l'ensemble de Mandelbrot (en évitant toutes les boucles for). On peut accélérer d'un facteur 10 et même plus, mais c'est délicat.