-- Introduction au Python scientifique --
Les bases du Python et de numpy
Introduction au python scientifique
-
simon.marie@lecnam.net
Ce notebook présente les bases du language Python qui seront utiles pour les TP du cours. Vous devez déjà avoir quelques notions en python. Ci ce n'est pas le cas vous pouvez par exemple regarder cette page.
Plan du Notebook:
Exécution d'un script Python
Le langage Python est un langage précompilé ne nécessitant donc aucune compilation. Il y a 3 façons principales d'utiliser Python:
- Ecrire un script dans un fichier .py et le lancer en ligne de commande:
\$> python script.py
Pour éditer des scripts, on utilise un éditeur spécialisé qui reconnait la syntaxe par exemple Visual Studio Code, Atom, SublimeText, vim, emacs, gedit, notepad++... - Utiliser directement Python en ligne de commande avec Ipython ou Spyder
- Utiliser un Notebook Jupyter (C'est ce moyen qui est utilisé ici):
\$> jupyter notebook filename.ipynb
Pour un projet large et complexe, on préfèrera la première solution, pour faire des petits tests ou pour une utilisation locale et restreinte on préfèrera la seconde solution. Pour les TP, les rapports et tous les support nécessitant une diffusion large et une utilisation de texte, le Notebook constitue la solution idéale.
Importation des modules
Les premères lignes d'un script Python sont consacrées à l'appel des modules nécessaires aux calculs. Les modules sont des librairies de fonctions qui permettent une utilisation plus simple de certaines notions. On utilise pour cela la commande import. Pour les besoins des TP de simulation numérique, les modules numpy, matplotlib et math seront nécessaires. Le début du script ressemblera donc à:
import numpy as np # contitent des fonctions utiles pour le traitement des matrices
import matplotlib.pyplot as plt # contitent des fonctions utiles pour la representation graphique
import time # Pour évaluer les temps de calcul
%matplotlib inline
# permet une manipulation des figures dans le notebook directement.
Les lignes ci-dessus importent le module numpy en le nommant np. Ainsi tous les appels à des fonctions propres aux modules numpy seront précédés de np.:
y=np.sin(np.pi) # calcul sin(pi)
print(y)
print(np.sqrt(np.log(np.exp(4))))
Pour obtenir de l'aide sur n'importe quelle fonction, il suffit de taper help(nom_fonction):
help(np.sin)
Manipulation des vecteurs / matrices
Par défaut, les listes créées nativement en python à l'aide des crochets ne sont pas très utiles:
A=[0,1,2,3,4,5,6,7,8]
B=[[0,1,2],[3,4,5],[6,7,8]]
print(A)
print(B)
print(B[1][1])
Les objets A et B sont des listes on ne peut pas faire d'opérations avec:
print(type(A),type(B))
print(2*A)
print(B*B)
Le module $\textit{numpy}$ permet la manipulation des matrices dans python. Avec numpy, on peut facilement transformer ces objets en matrices:
D=np.array(A)
E=np.array(B)
print(type(E))
print(E[1,1])
print()
print(2*D)
print()
print(E*E) # produit terme à terme ! ce n'est pas un produit matricielle
print()
print(np.dot(E,E)) # produit matriciel
Pour créer des vecteurs et des matrices on pourra donc utiliser des fonction natives du module numpy. Ainsi pour initialiser et créer une matrice pleine de 0 on utilisera la fonction $\textit{zeros}$ de $\textit{numpy}$.
nx=10
ny=4
A=np.zeros((nx,ny)) # Matrice de 10 lignes et 4 colonnes ne contenant que des 0
print(A)
La fonction arange(n) de $\textit{numpy}$ sert a créer une liste d'entier allant de $0$ à $n-1$.
B=np.arange(nx) # Vecteur de nx elements allant de 0 à 9
print(B)
Attention, en python les indices commencent à 0 ! Pour sélectionner un seul élément, il suffit de passer l'incie entre crochet:
$$ B[indice] $$print(B[0]) # Le premier élément
print(B[1]) # Le deuxième élément
print(B[nx-1]) # Le dernier
print(B[-1]) # Le dernier
print(B[nx-2]) # L'avant dernier
print(B[-2]) # L'avant dernier
On peut aussi sélectionner ou modifier une partie de matrice en sélectionnant une tranche. Pour cela on utilise la syntaxe général:
$$ Matrice[ début : fin : pas ] $$Si on ne précise pas de valeurs pour ces bornes, les valeurs extrêmes sont choisies par défaut comme présenté ci-dessous:
print(B[:]) # Tous les elements: équivalent de B[::]
print(B[0:]) # Equivalent à B[:]
print(B[0:nx]) # Equivalent à B[:]
print()
print(B[:5]) # 5 premiers Elements
print(B[2:5]) # tranche du 3ème au 5ème
print(B[-5:]) # 5 derniers
# B[n1:n2] # selection de n1 a n2-1 donc ligne n1+1 a n2
On peut enfin faire varier la valeur du pas qui peut même prendre des valeurs négatives:
print(B[::2]) # 1 élément sur 2 => équivalent de B[0:nx:2]
print(B[1::2]) # 1 élément sur 2 en partant du deuxième => équivalent de B[1:nx:2]
print()
print(B[::-1]) # tous les élements dans l'ordre inverse : pas négatif
Pour créer des vecteurs contenant des suites de nombres régulièrement espacés (par exemple pour des maillages) on peut utiliser:
- Soit la fonction $\textit{arange(start,end,step)}$ de $\textit{numpy}$ qui impose le début la fin et le pas:
A=np.arange(-2,2,0.5) # Vecteur contenant tous les nombres entre -2 et 1.5 par pas de 0.5
print(A)
- Soit la fonction linespace(start,end,nombre) de numpy qui impose le début, la fin et le nombre d'élément:
A=np.linspace(-2,2,10) # Vecteur contenant 10 élément régulièrement espacés entre -2 et 2
print(A)
Pour créer une matrice diagonale on utilise la fonction $\textit{eye(nx,ny)}$ de $\textit{numpy}$:
B=np.eye(3,3) # Matrice identite de 3*3 elements
print(B)
Enfin, pour créer une matrice complexe à partir d'un vecteur ou d'une autre matrice de taille différente, on peut utiliser une fonction très puissante, la fonction reshape. La seule condition est que le nombre d'élements total reste le même:
C=np.linspace(2.1,1.2,16)
D=C.reshape(4,4)
print(C)
print()
print(D)
Le premier objet est une matrice 1x16, le second une matrice 4x4.
Pour inverser une matrice on peut utiliser la fonction inv de linalg:
Dinv=np.linalg.inv(D)
print(np.dot(D,Dinv)) # produit matriciel avec l'inverse
Les boucles python
L’ensemble de la syntaxe python est basée sur l'indentation du texte. Ainsi toutes les lignes devront commencer au même endroit. Le début et la fin d'une boucle (for, if, while...) sera repérée par une indentation différente. C'est une syntaxe propre au Python qui permet de s'affranchir des syntaxe de type end for mais qui impose d'être rigoureux dans l'écriture du code ! Par convention on utilise une indentation de 4 espaces. Ainsi, une boucle for s'écrira de la façon suivante:
A=np.arange(10)
j=1
for i in A:
j+=2 # equivalent de j=j+2
print(j,i,i**2,i**3)
Attention a ne pas oublier le signe $:$ à la fin de la première ligne de boucle.
Notes sur les Boucles et le slicing en Python
Si l'on doit effectuer des boucles sur un très grand nombre de points, le temps de restitution risque de devenir très long. Dans ce cas on préferrera faire du "slicing" c'est à dire crire une seule ligne en décalant les indices. On écrit des tranches d'indices (slice en anglais).
L'exemple suivant illustre le même calcul réalisé en boucle puis en slicing. Pour le slicing, on remplace i par la tranche 1:-1 (du deuxième indice à l'avant dernier). Donc:
i => 1:-1
i+1 => 2:
i-1 => 0:-2
nx=1000000 # On prend un nombre important de points pour que la démonstation soit convainquante:
f1=np.zeros(nx)
f2=np.zeros(nx)
g=np.ones(nx)
# boucle
t0=time.time()
for i in np.arange(1,nx-1):
f1[i]=g[i+1]+g[i]+g[i-1]
t1=time.time()-t0
# slicing
t0=time.time()
f2[1:-1]=g[2:]+g[1:-1]+g[0:-2]
t2=time.time()-t0
print(t1,t2,t1/t2)
Ici le slicing est donc trois-cent fois plus rapide !!
On peut vérifier que les fonctions $f_1$ et $f_2$ sont bien identique
np.max(abs(f1-f2))
Utilisation des fonctions
Dans tous les langages de programmation, il est possible d'effectuer une tache répétitive à l'aide d'une fonction. En python on utilise pour cela le mot clé def suivit du nom de la fonction que l'on veut créer. Essayons par exemple de créer une fonction qui calcul la valeur d'un polynôme en x:
def ma_fonction(x):
return 2.*x**2-3.*x+1
print(ma_fonction(-1.))
print(ma_fonction(0.))
print(ma_fonction(1.))
On peut également créer des fonctions un peu plus évoluées comme par exemple la fonction factorielle qui utilise la récusrion (le fait d'appeller la fonction dans la définition de la fonction):
def fact(n):
if n==0:
return 1
else:
return n*fact(n-1)
fact(4)
Représentation graphique
Le module matplotlib permet l'utilisation des fonctions graphiques de type plot(), pcolormesh() ou contour().
plot¶
On peut par exemple tracer la réprésentation graphique de la fonction créée précédemment:
x=np.arange(-2.,2.,0.1)
plt.plot(x,ma_fonction(x))
On peut également créer des axes plus riches, ajouter des légendes et utiliser les notations scientifiques:
x=np.arange(-10,10,0.5)
z=np.exp(-x**2/16)
z2=np.exp(-x**2/9)
fig=plt.figure(1,figsize=(12,6)) # Creation de l'objet figure 1 on force la taille de la figure
plt.plot(x,z,label="$z=exp(-x^2/16)$") # Courbe de y en fonction de x
plt.plot(x,z2,label="$z=exp(-x^2/9)$") # Courbe de y en fonction de x
plt.xlabel('On peut etiqueter les axes $\epsilon$',fontsize=20.)
plt.ylabel('Utiliser des symboles: $\sum_\infty \epsilon$',fontsize=20.)
plt.title('Titre du graphique',fontsize=20.)
plt.grid(True);plt.legend(fontsize=20.)
Pour enregistrer une figure, on peut utiliser la fonction savefig:
fig.savefig('Ma_belle_figure.png', bbox_inches='tight',dpi=120) # Le format peut etre choisi librement (png,jpg,eps...)
pcolormesh¶
Pour les grandeurs 2D on utilise souvent des cartographies de couleur pour représenter l'évolution d'une grandeur dans le plan. Par exemple, en reprennant l'exemple précédent, on ajoute une dimension y et crée une donnée 2D z3:
y=np.arange(-10,10,0.5)
z3=np.einsum('i,j->ij',np.exp(-x**2/16),np.exp(-y**2/16))
print(z3.shape)
fig=plt.figure(2,figsize=(16,6)) # Creation de l'objet figure 2 en specifiant la taille
fig.add_subplot(121) # On crée une subfigure dans une division de 1 ligne et 2 colonne
plt.pcolormesh(x,y,z3.T,shading='nearest')
# Representation 2D de la matrice A avec la colormap "spectral"
# A.T signifie que l'on affiche la transposee de A
# vmin et vmax permettent de fixer l'echelle de couleur
plt.xlim((-10,10)) # Restriction de l'affichage a un domaine particulier
plt.ylim((-10,10))
# !! Attention en dehors des notebook il faut utiliser plt.show() pour afficher une figure
plt.xlabel(r'$x$',fontsize=18.)
plt.ylabel(r'$y$',fontsize=18.)
plt.colorbar()
fig.add_subplot(122) # Sur la deuxieme subfigure on changhe les options d'affichage
plt.pcolormesh(x,y,z3.T,shading='gouraud',cmap="jet") # on change la colormap et on lisse l'affichage
plt.xlim((-10,10)) # Restriction de l'affichage a un domaine particulier
plt.ylim((-10,10))
plt.xlabel('$x$',fontsize=18.)
plt.colorbar()
contour()¶
Lorsque l'on souhaite représenter les isovaleurs d'une nappe de données 2D, on peut utiliser la fonction contour ou contourf (contour remplis).
fig=plt.figure(2,figsize=(16,6)) # Creation de l'objet figure 2 en specifiant la taille
fig.add_subplot(121) # On crée une subfigure dans une division de 1 ligne et 2 colonne
plt.contour(x,y,z3,levels=[0.2,0.5,0.8],colors=['k','b','r'])
plt.axis('equal') # Restriction de l'affichage a un domaine particulier
# !! Attention en dehors des notebook il faut utiliser plt.show() pour afficher une figure
plt.xlabel(r'$x$',fontsize=28.)
plt.ylabel(r'$y$',fontsize=28.)
fig.add_subplot(122) # Sur la deuxieme subfigure on changhe les options d'affichage
plt.contourf(x,y,z3.T,cmap="Greys") # on change la colormap et on lisse l'affichage
plt.axis('equal') # Restriction de l'affichage a un domaine particulier
# !! Attention en dehors des notebook il faut utiliser plt.show() pour afficher une figure
plt.xlabel(r'$x$',fontsize=28.)
plt.colorbar()
Conclusion
Le python est un language simple très utilisé en science et dans plein d'autres domaines. Il permet de se familiariser facilement avec les problématiques liées à la programmation scientifique, au calcul mathématique et à la simulation numérique.
Note: Vous pouvez télécharger ce notebook (fichier ipynb) et l'executer (modifier) à votre guise sur une machine disposant d'une version de Python (3.x) ou bien l'éxecuter sur le JupyterHUB du CNAM
Pour aller plus loin voici quelques liens utiles:
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)