En physique, pour faire des simulations numériques, du post-processing de données expérimentales, même pour les calculs les plus simples, l'outil le plus répandu est matlab. Cependant, en plus d'être propriétaire, matlab est une grosse usine gaz, lourde, lente, et sous linux encore un peu buggé. Et puis le rendu des graphiques en 3D est franchement pas terrible. En plus de ça, le langage de programmation de matlab est très haut niveau et pas forcément bien optimisé pour les simulations.

Après être passé par Octave qui résoud les problèmes de licence et de lourdeur du logiciel, je me suis tourné vers python.

Pourquoi python pour des applications scientifiques ?

  • Je connais déjà un peu python et j'apprécie sa syntaxe simple et très lisible.
  • C'est un véritable langage de programmation, donc plus rapide que le langage matlab, avec des librairies de haut niveau et des librairies de bas niveau, ce qui permet d'intégrer ses programmes dans des applications d'intérêt plus large que le seul but scientifique (interface graphique, liaison avec le matériel, ...).
  • Python est moderne mais possède déjà une très grosse communauté, une documentation bien fournie et à jour, et surtout un très grand nombre de bibliothèques dans de très nombreux domaines, et notamment, numpy, scipy pour les maths et la simulation, matplotlib et gnuplot-py pour dessiner des graphiques.

J'ai donc fais un petit essai de simulation numérique d'un système dynamique simple mais qui donne des résultats intéressants : le modèle de Lorenz. C'est ce modèle qui permit la découverte et l'étude des comportements chaotiques, à l'origine développé par Lorenz donc, pour étudier les comportements météorologiques.

Pour les physiciens, il s'agit d'un système d'équations aux dérivées ordinaires du 3eme ordre : g12809.png

Pour résoudre ce système, il suffit de définir le système d'équation et d'utiliser la fonction d'intégration de scipy : integ.odeint(function, conditions initiales, temps). L'algorithme utilisé est très performant et rapide, largement plus qu'un RungeKutta 4 à pas variable dans matlab. Après ça, je mets quelques fonctions pour dessiner le portrait de phase et sauver les points dans un tableau.

Pour faire tourner ce code, il vous faudra python bien sur, numpy, scipy et gnuplot-py.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy
from scipy import arange, pi
import scipy.integrate as integ
import Gnuplot, Gnuplot.funcutils

sigma=10  #
b=8.0/3.0  #  Parameters for the system
r=28         #

y_0=[-10,-10,10] #  Initial Condition

def func(y,t):
    y0=sigma*(y[1]-y[0])
    y1=r*y[0] - y[1] - y[0]*y[2]
    y2=y[0]*y[1] - b*y[2]
    return [y0, y1, y2]

def resol():
    t=arange(0, 100.0, 0.01)
    y=integ.odeint(func, y_0, t)
    plot(y)
    tableau(y)
    
def plot(y):
    g = Gnuplot.Gnuplot(debug=1)
    g('set parametric')
    g('set data style lines') 
    g.xlabel('x')
    g.ylabel('y')
    g.zlabel('z')
    g('set border 1+4+16') #Choose which axis will be displayed
    g('set xtics nomirror; set ytics nomirror') #Ticks only on the front axis
    g('set ticslevel 0') #No offset for the ztics
    #g('set hidden')
    #g('set contour base')
    g.splot(Gnuplot.Data(y, binary=0))
    nom=raw_input("""Enter the name for the postcript image ("C" to cancel):\n""")
    if not (nom=='C' or nom=='c'):
        g.hardcopy(nom + '.ps', enhanced=1, color=1)
          
def tableau(y):
    table=raw_input("""Enter the name for the data file ("C" to cancel) :\n""")
    if not (table=='C' or table=='c'):
        numpy.savetxt(table,y)

if __name__ == '__main__':
    resol()

Les paramètres sont ceux utilisés par Lorenz lorsqu'il fit ce calcul. Une fois lancé, la simulation dure quelques secondes (pour calculer 10000 points !) et voici le résultat obtenu pour le portrait de phase : gp_lorenz.png

Joli, hein ! Un mouvement chaotique avec 2 ailes où les trajectoires tournent autour de 2 points fixes instables, et des passages d'une aile à l'autre complètement imprévisibles.

Conclusion : Python et ses librairies pour le calcul scientifique, c'est simple, efficace, libre et ça rend bien. Tout bénef donc !