Pourquoi en venir à Python, alors que Maple, Scilab, ou encore C ou FORTRAN sont des outils bien adaptés à nos besoins de programmation en physique numérique ?
Pour plusieurs raisons:
Je ne vous ferai pas un cours de Python. Vous trouverez tout ce qu'il faut sur le net pour installer Python sur votre Mac ou sur votre PC, sous Windows, Mac OSX ou Linux.
Vous trouverez également les tutorials nécessaires pour débuter. Je vous encourage à abuser du site officiel de Python.
Installer Python n'est pas très compliqué. Il suffit de se rendre sur le site officiel de Python, de sélectionner le produit à télécharger et de le télécharger. Vous aurez à choisir entre windows, linux ou Mac OS, et la version de Python. Il faut savoir qu'il y a un vrai gap entre la version 2.7 et la version 3.0. Dans un premier temps, je vous suggère de rester en 2.7.
Il vous faudra aussi charger les packages scientifiques dont on aura besoin : NumPy, SciPy et MatPlotLib. Ils sont disponible sur leur site, que vous trouverez immédiatement en faisant un coup de Google. Là aussi, il suffira de les télécharger et de déclencher l'installation (cela diffère en fonction de l'OS). Mais bon, je suppose que mes fidèles lecteurs ne sont pas des tanches en informatique...
Je vous recommande vivement de lire la suite de cette page après avoir ouvert une fenêtre commande Python (lancez l'IDLE puis faites Run/Python Shell dans le menu). Vous obtiendrez une fenêtre avec le prompt de commande de Python (>>>). N'hésitez surtout pas essayer les commandes mentionnées, profitez du fait que Python soit un interpréteur...
Un des grands attraits de Python, c'est que beaucoup de scientifiques l'utilisent et donc le nombre de librairies (je devrais écrire bibliothèques mais bon...) est important et leur nature très diversifiée.
Vous trouverez ci-dessous une liste des principales bibliothèques à usage scientifique. Elle est bien incomplète !
On trouve des librairies techniques:
et aussi des librairies métier:
Les types de variables Python qui nous intéressent en calcul sont : les "integer", les "long integer", les "float" et les "complex". On rencontera aussi des variables de type "list", "array" et "matrix". Il existe d'autres types (string, dictionary,tuple,list, etc.) dont vous pourrez prendre connaissance dans un cours Python.
Il est prudent de n'utiliser les opérateurs mathématiques pour faire des calculs que sur les quatre types mentionnés. Avec les autres types, les résultats peuvent être inattendus! Essayez par exemple de faire la somme de deux variables de type "list" !
les opérateurs de Python sont standard:
Attention, le comportement de la division est un peu curieux. Elle fonctionne normalement lorsque vous divisez deux floats, par exemple 5./2. donne comme résultat 2.5. Mais si vous divisez 5/2 (des integer donc) vous obtiendrez 2 ! Il restitue la partie entière du résultat sous forme d'un integer. Source d'erreur garantie !
Les règles de prévalence des opérateurs sont classiques.
Les opérateurs +=, -= , *= et /= sont disponibles. Par exemple, x = x + 1 pourra s'écrire x += 1. Rien de bien original...
Le package math de Python contient les principales fonctions mathématiques. Toutefois, je vous conseille d'utiliser plutôt le package NumPy beaucoup plus complet, car nous aurons à utiliser des fonctions absentes du package math.
Dans le package NumPy, on retrouve:
Vous trouverez dans la doc en ligne de NumPy les autres fonctions et la syntaxe précise des fonctions.
Il est tentant d'utiliser les possibilités "primitives" de Python pour définir une matrice, c'est à dire d'écrire simplement m1 = [[0,1], [1,2]] pour définir une matrice carrée d'ordre 2, c'est à dire une liste de listes.
Cela semble simple, et c'est la méthode que présente certains bouquins de Python. Et ça peut fonctionner... Avec des aléas ! Par exemple, si vous voulez calculer la somme de deux matrices, m1 = [[0,1], [1,2]] et m2 = [[1,0], [4,5]] par exemple, vous seriez tenté d'écrire m1 + m2. Quoi de plus naturel que de surcharger l'opérateur d'addition?
Sauf que ça ne fonctionne pas! Essayez donc, vous obtiendrez [[0, 1], [1, 2], [1, 0], [4, 5]] avec notre exemple, ce qui n'est pas vraiment le résultat attendu ! Pourquoi: parce que Python ne gère pas ses listes comme des matrices! Idem pour la soustraction et la multiplication par un scalaire. Je ne vous parle même pas de la multiplication de deux matrices.
Comment faire ? Il y a deux solutions :
Nous allons donc utiliser le package (la librairie c'est pareil) numpy.matlib, qui offre beaucoup de possibilités en calcul matriciel.
Dans votre programme, il faudra tout d'abord importer le package numpy.matlib, en le renommant np pour faciliter l'écriture, avec la commande
import numpy.matlib as np
Puis pour déclarer une matrice, vous écrirez simplement:
m1 = np.matrix([[0,1], [1,2]])
et
m2 = np.matrix([[1,0], [4,5]]).
Maintenant, si vous tapez:
print m1+m2
non seulement vous obtiendrez le résultat attendu, mais en plus écrit sous forme matricielle. Pas mal, non?
Pour accéder à l'élement (i,j) de la matrice, rien de plus simple, il suffit de le désigner par m1[i][j]. Par exemple, tapez print m1[1][1], vous obtiendrez 2. Si vous tapez m1[0][0], vous obtiendrez 0.
En utilisant le type matrix, rien de plus simple: vous tapez :
print m1*m2
et vous obtiendrez le résultat attendu [[4, 5] [9,10]].
La fonction np.transpose permet d'obtenir facilement la transposée d'une matrice. Essayez:
print np.transpose(m2)
Vous obtenez la matrice transposée de m2. On peut aussi écrire m2.transpose(). Essayez donc avec m1. Cette matrice a la curiosité d'être égale à sa transposée...
Pour inverser une matrice, supposée inversible, on utilise la fonction np.linalg.inv, qui est appelle à un algorithme de la librairie bien connue LinPack. Essayez:
print np.linalg.inv(m1)
vous obtiendrez la matrice inverse [[-2.,1.][ 1.,0.]]. Notez que les termes de m1 sont des entiers, alors que les termes de sa matrices inverse sont des flottants...
Le déterminant d'une matrice inversible se calcule avec la fonction np.linalg.det. Essayez:
print np.linalg.det(m1)
vous obtiendrez -1, le déterminant de la matrice m1.
La commande:
print np.linalg.eig(m1)
vous fournit les valeurs propres et les vecteurs propres de la matrice m1 : [-0.41421356, 2.41421356] et [[-0.92387953, -0.38268343],[ 0.38268343, -0.92387953]].
Nous avons déjà étudié sur une autre page de TangenteX la résolution d'un système linéaire de type Ax = B. NumPy offre un moyen très simple de les résoudre, après avoir saisie les matrices A et B.
Essayez les commandes suivantes pour décrire le système avec ses deux matrices et afficher le vecteur solution x:
A = np.matrix([[1,1,3],[2,2,5],[1,0,9]])
B = np.array([1,1,1])
print np.linalg.solve(A,B)
On obtient le vecteur solution [-8. 6. 1.]. Plutôt simple, non?
Pour info, numpy.linalg.solve utilise la méthode de la décomposition LU que nous avons étudié dans la page mentionnée ci-dessus.
Du point de vue de Python, un vecteur est représenté par une matrice unidimensionnelle. Je peux donc définir tranquillement un vecteur v1 par:
v1 = np.array([[0,1,1])
Pour la suite, je définirai de la même façon un vecteur v2 de coordonnées [1,0,1].
NumPy nous fournit bien sur de quoi calculer le produit scalaire de deux vecteurs de dimension quelconque. Par exemple, le produit scalaire de v1 et v2 se calcule et s'affiche par:
print np.inner(v1,v2)
Plutôt simple non? Le résultat est un réel(un float plus exactement), ici 0.
Le calcul d'un produit vectoriel est tout aussi simple. Essayez la commande:
print np.cross(v1,v2)
Le résultat est bien sur un vecteur (un pseudo-vecteur), dans notre exemple [ 1 1 -1].
Nous avons abordé sur la page les méthodes d'intégration numérique les différentes méthodes d'intégration d'une fonction, dont la méthode de Simpson et celle de Romberg.
Python, plus précisement le package scipy.integrate, nous fournit quelques fonctions pour intégrer très simplement, des intégrales simples ou multiples. Je ne vais pas toutes les reprendre ici, je vous invite à consulter la documentation scipy en ligne.
Signalons, et c'est important en physique numérique, que Python offre des outils pour intégrer des fonctions mais aussi des tableaux d'échantillons. Dans la pratique, il est en effet rare que nous connaissions la fonction à intégrer. Nous avons plutôt des séries de mesures ! Le package offre les mêmes méthodes pour les fonctions et les séries de mesures, j'en resterais donc aux fonctions.
Pour utiliser les fonctions suivantes, vous devez d'abord importer le package scipy.integrate avec la commande
import scipy.integrate as sp
Dans les exemples ci-dessous, j'ai choisi d'intégrer la fonction sin(x) entre 0 et pi, la simplicité quoi...
La fonction standard de sp pour calculer les intégrales simples est la fonction quad. Essayez la commande:
print sp.quad(np.sin,0,np.pi)
La fonction retourne deux données dans un tuple : la valeur estimée de l'intégrale, ici 2.0, et la valeur haute de l'erreur, ici 2.220446049250313e-14
Notez les préfixes np poir la fonction sinus et la valeur de pi, deux objets qui sont inclus dans le package NumPy.
Comment faire si vous vouliez intégrer une fonction de votre composition, nommons la "funcToto", que vous auriez bien sur défini auparavant ? Il suffit de passer son nom à la fonction quad, en lieu et place de np.sin !
Prenons la même fonction sin et intégrons la avec la méthode de Romberg sur les mêmes bornes. J'utilise la commande:
print sp.romberg(np.sin,0,np.pi)
Vous constaterez que le résultat est le même (ouf !) mais que l'erreur n'est pas indiquée. Par défaut, elle est bornée à 1.48e-8.
Les équations différentielles ordinaires (EDO ou ODE en anglais) sont très classiques en physique. Elles décrivent le comportement d'un système physique, très généralement simplifié, linéaire et ne dépendant que d'une variable. Dans la réalité, c'est plutôt rare !
Nous les avons déjà rencontré dans plusieurs expériences sur TangenteX, par exemple lors de l'étude de la désintégration radioactive . Dans cette même page, nous avions étudié une méthode numérique pour les résoudre, la méthode d'Euler, méthode que nous avions implémenté en C.
Vous vous souvenez que dans le programme scilab qui résolvait l'EDO, nous avions utilisé la fonction ode(), après avoir défini la fonction à intégrer. On retrouve dans beaucoup de langages une fonction similaire (maple, matlab, etc.).
Python est riche ! Il offre deux possibilités pour résoudre une EDO: une fonction nommée odeint() et une classe ODE. Nous ne parlerons pas ici de cette dernière, mais je vous encourage vivement à consulter la documentation à son sujet.
Bien sur, vous pouvez également coder votre méthode favorite en Python ou même appeler une fonction C déjà écrite (il y en a quelques unes qui trainent sur tangenteX...).
La fonction odeint() de Python est très similaire à la fonction ode() de scilab ou de maple. Nous avons déjà utilisé cette fonction pour simuler l'expérience de Millikan. Revoyons cela.
Il convient d'abord de définir l'équation différentielle que l'on veut intégrer. Nous devons la définir dans une fonction, ici Millikan, qui retournera pour chaque valeur de v, sa valeur à l'instant suivant:
# definition de l'équation différentielle
def Millikan(v,t):
return -K*v + C
K et C sont des constantes. v est la variable d'intégration. Puis on définit l'intervalle de temps sur lequel on intégre:
time = arange(t0, tmax, pastemps) # créé un vecteur de pastemps points entre t0 et tmax
La variable time est un vecteur contenant donc toutes les valeurs de temps utilisées lors de l'intégration de l'EDO. Enfin on intègre en utilisant la fonction odeint:
ve = odeint(Millikan,v0,time)
La fonction odeint() retourne un vecteur, nommé ici ve, de même dimension que time, qui contient les solutions. Il ne reste plus qu'à afficher ou à "plotter" ce vecteur, par exemple en utilisant l'instruction Python:
plot (time,ve,'g-')
La fonction odeint() est plus riche que ce que je viens de décrire. Vous vous en convaincrez en lisant la doc relative à celle-ci. En particulier, l'algorithme utilisé sait traiter les EDO "raides", c'est à dire celles qui décrivent un phénomène variant très brusquement dans le temps.
Nous avons déjà rencontré plusieurs fois sur ce site des EDO de second ordre, par exemple dans l'étude du pendule simple . Dans le code Scilab utilisé pour résoudre l'équation, nous avions utilisé une technique classique : transformer une EDO de second ordre en un système de deux EDO de premier ordre, ce qui nous avait permis d'utiliser la fonction ode() de Scilab.
Comment cela fonctionne-t-il? Reprenons l'équation différentielle du pendule:
\(\frac{d^2\theta}{dt^2} + \omega^2 \sin(\theta) = 0\)
Posons \(\frac{d\theta}{dt} = \theta'\) et \(\frac{d\theta'}{dt} = \omega^2 \sin(\theta)\). Nous obtenons un système de deux équations différentielles de premier ordre, qui est strictement équivalent à notre équation différentielle du second ordre.
Avec Python, nous allons traiter ce système selon le même principe qu'avec Scilab. La fonction \(\theta(t) \) et sa dérivée première \(\theta'(t) \) seront stockées dans un tableau de deux éléments que j'appellerai T\(\theta \). Posons que T\(\theta \)[0] = \(\theta(t) \) et T\(\theta \)[1] = \(\theta'(t) \). Notre système différentiel s'écrira donc : dT\(\theta \)[0]/dt = T\(\theta \)[1] et dT\(\theta \)[2]/dt = \( -\omega^2 \sin(T\theta[0])\).
Je vous propose un petit programme Python PendulePython, qui sera très similaire au programme scilab mentionné ci-dessus. Vous constaterez que ces deux programmes se ressemblent beaucoup...
Vous noterez en particulier la similitude dans la définition du système différentiel et l'usage des fonctions ode() et odeint(). Même les instructions de visualisation sont semblables.
Les équations aux dérivées partielles (PDE en anglais) sont d'un abord numérique un peu plus compliqué! Je survole le sujet dans cette page .
Il existe beaucoup de librairies Python spécialisées dans la résolution des PDEs selon leur type. Personnellement, je vous recommande la librairie FiPy, mais il y en a d'autres (Dolfin, MayaVi, GiNaC, etc.).
Ou 1D, cela dépend des auteurs ! Disons la courbe bête dans un repère Oxy, ce qui fait bien 2 dimensions si je ne m'abuse... Dans les petits programmes que nous avons écrit précédement, nous avons déjà eu l'occasion de tracer ces courbes.
Comme vous avez pu le noter, dans PendulePython.rar par exemple, j'utilise le module pylab de la librairie MatPlotLib. Pas question ici d'aborder toute la richesse de cette librairie, il existe sur le net un tas de tutorials et de manuels de référence pour se faire. Si vous ne trouvez pas, je vous en indiquerais.
Je vais plutôt vous commenter dans le source un exemple classique de tracé de deux courbes dans un repère, avec légende et labels des axes. Faisons dans le classique : je vais tracer les courbes de sinus et cosinus entre -pi et +pi. Voici ce que cela donne:
On peut enjoliver la chose, tracer les courbes dans deux fenêtres séparées (avec subplot) et faire plein d'autres choses. C'est décrit dans la doc...
Il existe d'autres librairies graphiques encore plus riches, mais MatPlotLib convient très bien à nos besoins...
ATTENTION : pour des raisons techniques, les scripts Python sont disponibles sous format rar compressé. Pensez à les décompresser...
Dans le même esprit, voyons comment tracer une courbe 3D (ou 2D selon les auteurs)... J'utilise le module mplot3d de MatPlotLib. Voici ce que le programme exemple nous donne pour le tracé de la courbe z = sin(sqrt(x2 + y2)):
Pourquoi ne pas revenir sur un grand classique de la mécanique : le pendule simple! Nous avons déja simulé avec Scilab et Maple le pendule simple. Nous pourrions recommencer avec Python, mais cette fois, je vais corser un peu (un tout petit peu) le problème: nous allons simuler un pendule simple amorti..
Soit donc un pendule de masse m, de longueur l, qui oscille dans l'air. On décide de ne pas négliger le frottement fluide de l'air sur la masse oscillante. Et l'on veut étudier l'évolution de l'angle \( \theta \) que fait le pendule avec la verticale.
Dans le cas d'amortissement nul, l'équation différentielle est très simple, comme nous l'avons vu ci-dessus :
\(\frac{d^2\theta}{dt^2} = -\omega^2 \sin(\theta)\)
Dans le cas d'un amortissement fluide, un terme supplémentaire s'ajoute à droite de l'équation pour donner:
\(\frac{d^2\theta}{dt^2} = -\omega^2 \sin(\theta) - K\frac{d\theta}{dt} \)
K est une constante, qui possède la dimension de l'inverse d'un temps, et qui vaut K = k/ml2, avec k le coefficient d'amortissement fluide.
Je me propose d'écrire un petit programme Python qui va simuler le comportement de notre pendule et aussi d'en étudier les caractéristiques. Une simulation numérique nous permettra de nous affranchir de l'approximation des petits angles, et donc de fournir des résultats plus réalistes.
Je suis parti du programme PendulePython dont j'ai pour l'instant modifié seulement le système différentiel, en ajoutant le terme d'amortissement:
def Pendule(theta,t):
return array([theta[1], -Omega2*theta[0] - K*theta[1]])
Et j'obtiens la courbe d'oscillation suivante, avec \(\theta_0 = 1 \) rd.s-1, une vitesse initiale nulle et un coefficient d'amortissement égal à 0.5:
Vous notez que, comme attendu, l'amplitude des oscillations décroit rapidement. Vous pouvez d'ailleurs faire plusieurs essais avec des valeurs différentes de K pour évaluer son influence sur la décroissance de l'amplitude. Dans notre cas, le pendule oscille en régime pseudo-périodique. Cherchez donc les conditions pour qu'il oscille en régime apériodique et critique...
Je vais maintenant m'intéresser au portrait de phase de mon pendule simple amorti. Pour rappel, un portrait de phase est la courbe \(\frac{d\theta}{dt} = f(\theta) \). Vous vous souvenez que dans le cas d'un pendule simple, cette courbe est une courbe fermée (voir ici). Voyons ce que cela donne dans le cas d'un pendule simple amorti, avec les mêmes valeurs des paramètres que ci-dessus:
La courbe est ouverte, en forme de spirale, qui plonge vers le point (0,0), le point attracteur, qui est la position d'équlibre stable. Vous pouvez vérifier l'influence qu'a le régime d'oscillation sur la forme du portait de phase...
Voilà un petit exemple simple de ce qu'on peut faire avec Python. Le code est simple, les outils puissants et donc, on peut se consacrer à ce que nous aimons vraiment : la physique. Je prépare quelques autres programmes pour essayer d'autres outils, la FFT par exemple...
L'objectif de cette page est de vous donner quelques outils pour faire de la physique numérique avec Python. Et aussi de vous faire découvrir, peut-être, un langage puissant et assez plaisant à utiliser.
J'espère avoir atteint cet objectif. Je vous encourage à vous plonger dans les manuels de référence des différents packages cités. C'est le seul moyen, avec la programmation évidemment, de vous familiariser avec eux. Bon courage !
Contenu et design par Dominique Lefebvre - tangenteX.com mars 2014 -- Vous pouvez me joindre par ou sur PhysiqueX
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France .