Le langage Python se répand à grande vitesse dans la communauté scientifique et c’est une bonne chose. C’est un langage simple, universel et d’appropriation rapide. Les étudiants en physique et les élèves de prépa sont largement encouragés à s’y investir pour acquérir les bons réflexes d’algorithmique et de programmation.
Bien, mais notre préoccupation majeure est quand même de faire de la physique. Et donc, au-delà de la phase indispensable d’apprentissage des algorithmes de base en calcul numérique, il ne me semble pas indispensable de coder et surtout de debugger et d’optimiser toutes nos routines de calcul. D’autant que les mêmes opérations reviennent très souvent : dérivation, intégration, résolutions d’équations, etc.
Il fut un temps, assez lointain quand même, où le numéricien passait un certain temps à écrire et à mettre au point sa propre librairie de routines, avec plus ou moins de bonheur et l’appui précieux des Numerical Recipes. Ce temps est fini. Vive les packages de calcul de Python, bien écrits, optimisés et bien documentés.
Voici donc le package scipy, l’un des plus classiques du calcul numérique, suffisamment complet pour les besoins des numériciens débutants, et très répandu. Nous allons y piocher touts les outils dont nous avons besoin. Pour information, le package scipy utilise un package « père », numpy, qui définit essentiellement les tableaux multidimensionnels et les fonctions de base de calcul.
Nous utiliserons également dans notre boîte à outils le package matplotlib.pyplot qui, comme son nom l’indique, est un package graphique dont les fonctions sont très comparables à ce que vous trouverez dans Matlab ou SciLab.
Il faut bien sur avoir d’abord installé Python sur sa machine, car l’installation de scipy se fait dans l’arborescence de Python. Pour ma part j’utilise Python 2.7.8 et donc j’ai installé le package correspondant à cette version. Adaptez-le à votre version de Python.
Pour télécharger et installer scipy, rendez-vous sur le site SciPy. Bien sur, scipy est un logiciel open source, gratuit, sous licence BSD.
J’utilise IDLE ou Spyder 2 dans Anaconda (sur mon Mac) pour travailler en Python. Pour contrôler la bonne installation de scipy, je vais dans la fenêtre 'Python shell' et je tape les commandes :
import scipy
puis
help (scipy.integrate)
par exemple. Vous devez voir apparaître un long texte d'aide. Sinon, c'est que le module scipy n'est pas installé ou mal installé.
Evidemment, vous devrez adapter cette méthode à votre environnement de développement sous Python…
Comme d’habitude en Python, pour utiliser les fonctions d’un package, il faut d’abord indiquer à Python le package et éventuellement les fonctions utilisées. Vous pouvez par exemple écrire en tête de votre script :
from scipy import *
ça marche mais ce n’est pas conseillé : vous venez d’importer les très nombreuses fonctions de scipy dans votre espace de travail. Cela consomme beaucoup de ressources système, la mémoire en particulier, sans aucune raison ! D’autre part, certaines fonctions existent avec le même nom dans plusieurs modules (numpy et math, par exemple), ce qui peut provoquer des choses bizarres ! Il est de beaucoup préférable de préciser les fonctions que vous voulez importer, par exemple :
from scipy import odeint
pour importer la fonction de résolution d’une équation différentielle de scipy.
Une boîte à outils, c’est une boîte dans laquelle on pioche l’outil nécessaire pour réaliser un travail donné. Vous aurez compris que cette boîte sera le package scipy. Mais cette boîte est vaste, et il faut trouver l’outil adéquat pour faire notre boulot.
Dans la suite, je vais donc partir d’un besoin : calculer une dérivée, intégrer une fonction, chercher la racine d’une autre, pour présenter l’outil, la fonction scipy, qui rend le service attendu. Et pour que cela soit plus parlant, j’illustrerai l’usage de cet outil dans un exemple, de physique bien sur !
Si vous souhaitez davantage de renseignements sur chacune des fonctions de scipy, je vous invite à consulter la documentation en ligne sur SciPy/Docs
Tous les scripts Python étudiés ci-dessous sont disponibles au téléchargement, voir à la fin de cette page.
Vous pouvez parcourir cette page, un peu longue, en séquence. Vous pouvez aussi accéder directement au paragraphe qui vous intéresse :
Bonne lecture !
Je vais commencer par l’outil de tracé de courbes en 1D car nous nous en servirons dans tous les exemples ! Il s’agit en fait de tracer la courbe représentative d’une fonction y = f(x) dans un certain domaine choisi, en précisant les caractéristiques des axes (type, graduations, légende, couleur), et mettant un titre, une légende, etc.
Nous devons d’abord déterminer le domaine de variation de x, en définissant les bornes inférieure et supérieure du domaine, que je noterai debut et fin. Nous devons aussi définir l’intervalle pas entre chaque point dans le domaine, car nous sommes dans un domaine discret : notre ordinateur ne sait pas calculer en réels ! En pratique, nous désirons donc obtenir un tableau x = [debut,fin] de (fin – debut)/pas éléments. Nous utiliserons pour cela la fonction arange(debut,fin,pas), qui nous retournera le vecteur x.
Pour tracer la courbe 1D, j’utiliserai différentes fonctions de matplotlib.pyplot, que je vous laisse découvrir dans l’exemple ci-dessous.
Imaginons donc que nous voulions tracer la courbe représentative de la fonction sinus entre 0 et 2*pi. Avant de tracer cette courbe, il faut importer les fonctions nécessaires puis la définir, ce qui est l’objet des instructions suivantes :
from scipy import arange, pi, sin
from matplotlib.pyplot import *
Je définis maintenant le vecteur x à l'aide de la fonction arange :
debut = 0.0
fin = 2*pi
pas = 0.01
x = arange(debut, fin, pas)
et je calcule le vecteur y = sin(x) :
y = sin(x)
Puis nous allons définir la fenêtre du tracé et les quelques fioritures : titre de la courbe, nom des axes, graduation, à l’aide des instructions :
figure()
title("Fonction sinus(x)")
xlabel('x')
ylabel('y')
grid(True)
axis([x.min(),x.max(), y.min(), y.max()])
Notez dans cette dernière instruction, l'utilisation des fonctions min() et max() sur un vecteur, qui permettent de définir les axes Ox et Oy précisément.
Enfin, nous allons tracer la courbe avec les instructions :
plot(x,y, 'r')
show()
Dans la fonction plot(), je définis la couleur du tracé, ici en rouge (red). Les couleurs sont désignées par leur initiale, en anglais!
Et voilà ce que cela donne :
Il s’agit de tracer une courbe z = f(x,y), une surface équipotentielle par exemple.
Pour tracer des surfaces en 2D et plus encore, nous utiliserons un package particulier de matplotlib, qui s'appelle mpl_toolkits. C'est un package très complet, assez pointu. Vous trouverez un tutoriel et la documentation sur le site qui lui est consacré.
Pour ce que je veux en faire ici, le principe d'utilisation n'est pas très différent de ce que nous avons vu plus haut. On définit une figure muni d'un système de 3 axes, puis on trace les trois vecteurs x, y et z qui définissent la surface.
Pour illustrer le tracé d'une surface avec Python, j'ai choisi de tracer la surface représentative de la fonction \( z = \sin(x^2 + y^2)\), pour x et y variant entre -pi/2 et +pi/2. Rien de bien passionnant, mais c'est un exemple...
Je commence par importer les fonctions nécessaires :
from scipy import arange, meshgrid,sin, pi
from matplotlib.pyplot import *
from mpl_toolkits.mplot3d import Axes3D
Puis je définis les vecteurs x et y :
xmin = -pi/2.0
xmax = pi/2.0
ymin = xmin
ymax = xmax
pas = 0.1
x = arange(xmin, xmax, pas)
y = arange(ymin, ymax, pas)
et je calcule le vecteur \( z = \sin(x^2 + y^2) \) :
x,y = meshgrid(x,y)
z = sin(x**2 + y**2)
A noter la fonction meshgrid() qui construit une matrice de coordonnées (dimension 2) à partir de deux tableaux de dimension 1, opération indispensable pour que la fonction plot_surface fonctionne.
Je prépare maintenant la figure et les axes à tracer. C'est un peu plus compliqué qu'en 1D, mais cela est quand même ressemblant :
fig = figure()
Fig3D = Axes3D(fig)
title('Fonction z = sin(x^2 + y^2)')
Fig3D.set_xlabel('x')
Fig3D.set_ylabel('y')
Fig3D.set_zlabel('sin(x^2 + y^2)')
Je peux enfin tracer la surface par les instructions:
Fig3D.plot_surface(x,y,z,rstride=1,cstride=1, color = 'r
show()
Et voilà ce que cela donne :
Chercher le ou les zéros d’une fonction f revient à résoudre l’équation f(x) = 0. Cela peut être trivial, si on sait le faire analytiquement, par exemple dans le cas d’un polynôme de degré 1 ou 2, plus difficilement pour les degrés supérieurs. Et cela peut être très complexe et impossible analytiquement dans le cas de fonctions non linéaires ou carrément exotiques.
Or, nous retrouvons ce besoin assez fréquemment en physique numérique. Il existe plusieurs méthodes numériques pour rechercher les racines d’une équation, nous en avons vu quelques exemples sur cette page. scipy nous évite de chercher plus loin…
Il faut faire attention dans la recherche des racines d'une fonction ! Certaines peuvent nous réserver des surprises désagréables. Un conseil donc : tracez d'abord la fonction pour vérifier visuellement sa continuité, la présence de racines dans le domaine, mais pas une infinité , et autres comportements exotiques. Bon, en physique, c'est assez rare, mais on ne sait jamais !
Le package scipy.optimize fournit une fonction nommée bisect() qui permet de calculer la ou les racines d'une fonction sur un intervalle donné. Vous trouverez la documentation de cette fonction sur la page qui lui est consacrée. Cette fonction est performante mais présente un inconvénient certain : il faut avoir une idée relativement précise de l'intervalle dans lequel se trouve la racine recherchée, ce qui signifie qu'il faut d'abord se faire une idée de la fonction, en la traçant par exemple...
Il existe une autre fonction plus puissante, la fonction fsolve(), qui ne necéssite pas de connaitre l'intervalle dans lequel la fonction change de signe. Elle utilise l'algorithme de Newton-Raphson (si on lui fournit la dérivée de la fonction) ou l'algorithme de la sécante.
Dans les exemples, nous utiliserons les deux fonctions pour les comparer.
Pour illustrer l'utilisation de ces fonctions, je vais calculer le temps de vol d’un projectile avant sa chute.
Imaginons un projectile lancé avec une vitesse initiale v0 et un angle initial \( \alpha_0 \). Le frottement de l’air est négligé. Il est lancé d’une altitude z0 = 0. Le référentiel d’étude est le référentiel terrestre supposé galiléen. Je cherche à établir le temps de vol du projectile avant sa chute sur le sol.
Je suppose que vous savez établir l’équation de la trajectoire de ce projectile, qui est \( z = -\frac{1}{2}gt^2 + v_0\sin(\alpha_0)t + z_0 \). Notre problème consiste donc à chercher pour quelle valeur de t, z = 0, en éliminant la solution triviale t = 0 (avec z0 = 0).
Pour traiter ce problème, je vais d’abord utiliser la fonction bisect(). Comme je l’ai dit plus haut, il faut indiquer à la fonction bisect() un intervalle de recherche dans lequel le signe de la fonction s’annule. Ce qui implique que l’on ait une petite idée de la « gueule » de la fonction. Le meilleur moyen est d’abord la tracer…
Le script BObisect.py met en oeuvre la fonction bisect() pour traiter le problème. Voyons ses caractéristiques principales. Pour commencer, j'importe les librairies dont j'ai besoin, rien de nouveau:
from scipy import sin, pi, linspace
from scipy.optimize import bisect
from matplotlib.pyplot import *
puis je définis les constantes et paramètres de l'expérience:
g = 9.81
v0 = 10.0
z0 = 0.0
alpha0 = pi/4.0
Maintenant, je définis la fonction dont je cherche la racine, en fait l'équation de la trajectoire z = f(t):
def f(t):
z = -0.5*g*t**2 + v0*sin(alpha0)*t + z0
return z
Je construit le vecteur temps entre 0 et 2 secondes, avec 100 points:
t0 = 0.
tmax = 2.0
nbpoint = 100
t = linspace(t0, tmax, nbpoint)
puis je calcule z = f(t) pour l'affichage :
z = f(t)
et enfin je calcule ma racine par bisect(). Remarquez que je lui fournis l'intervalle dans laquelle je suppute que la racine se trouve. J'ai estimé cet intervalle en traçant la fonction avec les instructions qui suivent...
tvol = bisect(f,1.0,2.0)
print tvol
Les instructions suivantes sont bien connues maintenant: elles tracent sans fioriture la courbe représentative de la fonction f(t).
figure()
grid(True)
plot(t,z)
show()
Le temps de vol obtenu est de 1.44160403912 s.
La courbe obtenue nous indique qu'en effet, la racine est bien dans ce coin...
Le script BOfsolve.py ne différe de BObisect.py qu'en deux lignes. D'abord dans l'importation des librairies qui devient:
from scipy import sin, pi, linspace
from scipy.optimize import fsolve
from matplotlib.pyplot import *
puis dans le calcul de la racine, où j'appelle fsolve() plutôt que bisect() :
tvol = fsolve(f,1.0)
Je dois fournir à fsolve() un point proche de la racine, j'ai ici choisi arbitrairement 1.
Après exécution, le temps de vol obtenu est de 1.44160404 s, identique à celui obtenu par bisect().
C'est un usage détourné de la fonction fsolve(). Imaginons que vous vouliez calculer la valeur de \( \sqrt(2) \) et que vous n'aviez pas de calculette sous la main, mais seulement Python. Il suffit de remarquer que \( \sqrt(2) \) est la solution de l'équation \( x^2 - 2 = 0 \) ! Vous la traitez par fsolve et le résultat est garanti. Voyez le script BORacine2.py, qui est simplissime. Les parties significatives du code sont :
def f(x):
return x**2 - 2
et :
racine2 = fsolve(f,1.0)
print racine2
Plutôt simple non, et vous noterez que le résultat obtenu est très bon !
La modélisation d'un phénomène physique aboutit souvent à l'écriture d'un système d'équations linéaires. Souvenez-vous de ces problèmes d'électricité remplis de noeuds et de mailles, dans lesquels il fallait calculer des courants ou des tensions... Lorsque le système ne comprend que deux équations à deux inconnues, ça va encore, mais au delà, le problème se complique!
On trouve aussi le besoin de résoudre un système lorsqu'on veut calculer le déterminant d'une matrice, inverser une matrice, résoudre un simplexe et que sais-je encore! Bref, il est bon de savoir comment résoudre un système d'équations. Python va nous y aider.
Le package scipy.optimize nous propose la fonction fsolve(), qui a de multiples usages comme nous l'avons vu ci-dessus, dont la résolution des systèmes d'équations. Il suffit de lui passer en paramètres la définition de votre système d'équations sous la forme \( a_ix_i = 0 \) et le vecteur d'initialisation xi0 et le tour est joué: elle vous retourne la solution sous la forme d'un vecteur contenant les xi.
Vous trouverez la documentation de la fonction fsolve() sur cette page .
Imaginons que vous étudiez un circuit électrique alimenté en 5 volt et composé d'un réseau de résistances dont vous connaissez les valeurs (ce sont les constantes de mon système). Les lois d'Ohm et de Kirchhhoff vous permettent d'écrire le système linéaire suivant:
\( \left\{ \begin{array}{r c l} 5 - 3i_1 - 6(i_1 - i_2) &=& 0 \\ 6(i_1 - i_2) - 10i_2 -5(i_2 +i_3) &=& 0 \\ -5(i_2 +i_3) - 3i_3 &=& 0 \end{array} \right. \)
A vous de trouver le circuit si cela vous amuse... En attendant, voyons comment résoudre ce système, c'est à dire trouver les valeurs i1, i2 et i3.
Le script BOsyslin.py va nous aider à résoudre ce système en un tour de main! Il est plutôt simple. Je passe sur l'importation des librairies, il suffit d'importer dans le package scipy.optimize, la fonction fsolve().
Vient ensuite la définition du système à résoudre, dans la fonction systeme():
def systeme(var):
i1,i2,i3 = var[0], var[1], var[2]
eq1 = 5 - 3*i1 - 6*(i1 - i2)
eq2 = 6*(i1 - i2) - 10*i2 - 5*(i2 + i3)
eq3 = -5*(i2 + i3) - 3*i3
r = [eq1, eq2, eq3]
return r
La fonction fsolve() a besoin de valeurs initiales pour démarrer son calcul. Je lui fixe arbitrairement 0 pour les 3 courants:
init = [0,0,0]
Puis je lance le calcul et j'affiche la solution :
S = fsolve(systeme, init)
print S
Pour information, le résultat obtenu pour les 3 courants i1, i2 et i3 est : [ 0.71571572 0.24024024 -0.15015015].
Le calcul d'une intégrale est un acte quotidien du physicien numéricien. Il peut être confronté à deux situations:
L'intégration numérique est une opération délicate dès que l'on s'éloigne des fonctions triviales. Les mathématiciens ont développé un nombre considérable de méthodes pour intégrer dans les cas les plus exotiques, des fonctions rencontrées parfois en physique. Dans cette page, nous en resterons aux cas simples, avec des fonctions qui sont intégrables au sens de Riemann, sans discontinuité ou divergence...
scipy met à notre disposition plusieurs fonctions pour calculer des intégrales simples, double ou triples. Vous en trouverez la liste en tapant :
help(scipy.integrate)
dans votre shell Python, après avoir chargé scipy bien sur !
Pour notre besoin, j'ai choisi d'utiliser la fonction quad(), basée sur une librairie FORTRAN très réputée dans le milieu. Vous trouverez la documentation de la fonction quad() sur cette page .
Son utilisation est des plus simples : vous lui passez le nom de la fonction à intégrer, que vous aurez défini au préalable et les bornes d'intégration, qui peuvent être +infini ou -infini. Eventuellement, les paramètres de la fonction à intégrer, et voilà ! quad() vous retourne un tuple (un vecteur) contenant, entre autres, le résultat de l'intégrale et une estimation de l'erreur commise.
Pour commencer, un exemple très simple. Supposons que je veuille calculer l'intégrale de la fonction sin(x) entre 0 et pi/2. Voyons comment procéder avec la fonction quad().
Le script BOintegrationSinus.py va nous permettre de faire simplement ce calcul. Il commence comme d'habitude par l'importation des librairies. Dans notre cas, la librairie particulière est :
from scipy.integrate import quad
Puis je définis la fonction à intégrer :
def f(x):
return sin(x)
et j'appelle la fonction quad() en lui passant en paramètres le nom de la fonction à intégrer, soit f dans notre cas, et les bornes d'intégration :
t = quad(f,0,pi/2)[0]
Une précision, la fonction quad() retourne un tuple qui contient le résultat de l'intégrale et la précision. Comme je désire ne récupérer que le résulat, je le dis au programme en indiquant [0], qui désigne le premier terme du tuple qui est le résultat de l'intégration.
Enfin, j'affiche le résultat :
print 'integrale entre 0 et PI/2 de sin(x) : ',t
Vous pouvez modifier ce script pour intégrer d'autres fonctions. Il vous suffit de modifier la définition de f(x) et aussi de ne pas oublier d'importer les fonctions nécessaires de scipy.
Faites attention tout de même à ne choisir que des fonctions qui soient intégrables sur le domaine que vous définirez, et aussi à éviter les fonctions "bizarres". Si c'est le cas, scipy vous retournera un message d'alerte en vous expliquant ce qui ne va pas. Essayez, par exemple, d'intégrer tang(x) entre 0 et PI/2...
Nous avons déjà étudié le pendule simple. Nous savons que l’équation différentielle qui décrit son mouvement est de la forme \( \frac{d^2\theta}{dt^2} + \omega_0^2 \sin(\theta) = 0 \). Il s’agit d’une équation non linéaire que l’on a l’habitude de simplifier en \( \frac{d^2\theta}{dt^2} + \omega_0^2 \theta = 0 \) en faisant l’approximation dite « des petits angles », qui pose \( \sin(\theta) \approx \theta \). Sous cette approximation, la période du pendule s’exprime simplement \( T_0 = 2 \pi \sqrt(\frac{l}{g}) \).
Mais qu’en est-il de la période T du pendule dans le cas où l’amplitude ne respecte plus l’approximation des petits angles ? Pour se faire une idée, nous allons étudier le rapport T/T0 pour un angle initial \(\alpha_0 \) variant entre 0 et pi/2. Il faut donc calculer T par chaque valeur de l’angle \(\alpha_0 \) que nous choisirons.
Je ne vais pas vous infliger la démonstration, qui n’est pas très difficile en partant de considérations énergétiques, en vous donnant l’équation qui permet de calculer T. Il s’agit, vous l’aurez deviné, d’une intégrale : \( T = 4 \sqrt(\frac{l}{2g}) \int_0^{\theta_0}\frac{d\theta}{\sqrt(\cos(\theta) - \cos(\theta_0))}\)
Cette intégrale se calcule très bien avec la fonction quad(). Voici donc le script BOintegration.py qui, pour différentes valeurs d’angle initial comprises entre 0.01 et pi/2 radians, calcule T et trace les deux courbes T et T0 en fonction de l’angle initial. Nous verrons ainsi l’erreur que l’on commet en adoptant l’approximation des petits angles, en fonction de l’angle initial.
Analysons les particularités de ce script. Je commence par définir la fonction qui calcule la période en fonction de l'angle \( \theta \), après l'appel des librairies bien sur :
def f(theta,theta0):
x = 1.0/sqrt(cos(theta) - cos(theta0))
return x
Je définis ensuite le domaine de variation de l'angle initial, fixé entre 0.01 et pi/2 et je calcule la période T0 avec l'approximation des petits angles. Rien de particulier. Vous noterez que j'ai posé k = l/g = 1 pour simplifier le calcul.
domaine = arange(0.01, pi/2.0, 0.01)
T0 = [2*pi*sqrt(k)]*len(domaine)
Cette dernière instruction me permet de construire un vecteur de même dimension que domaine, afin de tracer la courbe T0 constante.
Les intructions suivantes sont le coeur du script : il s'agit de la boucle de calcul de la période T pour chaque valeur de l'angle initial. Je définis d'abord le tableau T de stockage des résultats, puis la boucle de calcul de T avec l'appel à quad().
T = [0.0]*len(domaine)
i=0
for angle in domaine:
T[i] = (4.0/sqrt(2))*quad(f,0,angle,args=(angle))[0]
i = i+1
Je ne reviens pas sur les instructions suivantes qui permettent le tracé des courbes T et T0, qui sont classiques.
Voici l'affichage obtenu avec ce script:
Vous constatez sur cette courbe que l'erreur commise sur la période avec l'approximation des petits angles est très vite significative ! A ne pas oublier lorsque dans un problème vous faites cette approximation. On dit classiquement qu'elle est valable jusqu'à 10°, soit environ 0,2 radian. On observe sur la courbe que c'est en effet le domaine acceptable pour cette approximation.
A propos de l’approximation des petits angles, j’aimerais vous soumettre un calcul très simple qui permet d’estimer l’erreur relative que nous commettons en la faisant.
L’idée de base part du développement limité de la fonction \( \sin(\theta) \). Les étudiants de prépa et de licence le savent, on peut démontrer l’égalité suivante : \( \sin(\theta) = \theta -\frac{\theta^3}{6} + O(\theta^5)\). On appelle cette égalité un développement limité (DL), dans lequel \( O(\theta^5)\) représente un reste d'ordre 5 négligeable.
Partant de ce DL, je peux écrire : \( \sin(\theta) - \theta = -\frac{\theta^3}{6} + O(\theta^5)\), ce qui me permet d’affirmer que \(|\sin(\theta) - \theta| <= -\frac{|\theta|^3}{6} \).
Imaginons maintenant que je fixe l’erreur acceptable pour l’approximation des petits angles à 10^-2. L’équation précédente me permet facilement de calculer l’angle (en radian bien sur !) auquel correspond cette erreur, en posant que \( 10^-2 <= -\frac{\theta^3}{6} \), soit \( \theta^3 = 6*10^-2 \) ou encore \( \theta = 0,391 \) rad.
Donc, pour un angle inférieur ou égal à 0,391 rad (22° environ), l’approximation des petits angles est acceptable à 10^-2 près. Ce qui n’est pas si mal ! Vous commettez une erreur relative de 10^-3 si vous vous limitez à un angle maximum de 0.182 rad soit 10° environ. C’est l’erreur communément acceptée en calcul de physique.
Comme l'intégration, la dérivation d'une fonction est une opération de routine en physique numérique. Elle sert à presque tout! Et comme pour l'intégration, on peut se retrouver à devoir dériver une fonction donnée sous sa forme analytique ou bien une série de nombres issus de mesures. Dans ce dernier cas, comme pour intégrer, il faut d'abord interpoler les mesures pour obtenir une fonction, que l'on dérivera ensuite.
Notre besoin se limitera donc ici à dériver, sur un intervalle donné ou bien en un point, une fonction, que l'on supposera continue et dérivable.
Le package scipy propose une routine de calcul de dérivée d'ordre n en un point, construite sur l'algorithme des différences centrées. Il s'agit de la routine derivative() que vous trouverez dans scipy.misc.
Les principaux arguments sont : le nom de la fonction à dériver, le point auquel vous voulez la dériver, le pas de calcul et l'ordre de la dérivation. Sur ce dernier point, assurez-vous que votre fonction est bien de classe Cn, sinon vous aurez droit à un beau message d'erreur. Attention au pas de calcul dx, qui doit être assez petit, le "assez" dépendant de la fonction. Faites plusieurs essais...
Vous trouverez la documentation de la fonction derivative() sur cette page .
Commençons par un exemple très simple : imaginons que je veuille calculer la dérivée de la fonction f(x) = sin(x) sur l'intervalle de 0 à 2*PI. Cette fonction est continue et dérivable sur le domaine, donc pas de problème.
Comme nous l'avons vu ci-dessus, la fonction derivative() calcule la dérivée de f(x) en un point. Il faut donc découper le domaine en un grand nombre de points, calculer la dérivée en chaque point et construire la courbe dérivée à l'aide de ces calculs. Pas très pratique...
Heureusement, Python qui sait travailler avec des vecteurs nous permet de réaliser ce calcul très simplement. Voyons comment dans le script BODerivationSinus.py.
Pour commencer les habituels appels de librairie :
from scipy import sin, pi, arange
from scipy.misc import derivative
from matplotlib.pyplot import *
La fonction derivative() est dans le package scipy.misc.
Je passe sur la définition de la fonction à dériver, que vous devez maitriser maintenant.
Je construit maintenant le vecteur x, à l'aide de la fonction arange(). Puis je calcule le vecteur y = f(x):
x = arange(0,2*pi,0.01)
y = f(x)
Et maintenant l'objet de la manoeuvre : je vais dériver la fonction f(x) pour tous les points du domaine et construire le vecteur yprime qui contiendra ces dérivées. J'utilise la fonction derivative(), en lui indiquant que je souhaite une dérivée première (le n = 1) et que le pas de dérivation est égal à 0.01 (le dx = 0.01). C'est bien sur le même que pour arange(). Cela donne :
yprime = derivative(f,x,n=1,dx=0.01)
Ensuite viennent les instructions de tracé de courbe, que vous maitrisez également, rien de nouveau.
Et voilà ce que cela donne :
Vous reconnaissez bien sur la courbe représentative de cos(x) en la courbe de la fonction dérivée !
Le calcul des dérivées intervient en optimisation lorsqu’il s’agit de rechercher le minimum d’une fonction. Il est alors associé au calcul des racines d’une équation que nous avons vu plus haut. Voyons un exemple simple et très classique. Une usine fabrique des boites de conserve de 1 l en métal et voudrait économiser le plus possible sur ses coûts de fabrication. On va supposer que les coûts de procédé ont été optimisés et qu’il ne reste plus qu’à travailler sur les coûts de matières premières. L’objectif est donc de réduire au mieux la surface de métal nécessaire pour englober un volume d’un litre, en supposant l’épaisseur du métal constante.
Appelons h la hauteur de la boite et r son rayon. La surface du cylindre que constitue la boite s’exprime donc par \(S = 2 \pi r h + 2\pi r^2 \) et son volume \( \pi r^2h \), si j’exprime les grandeurs en décimètres. Rappelons que 1 litre = 1dm^3 et donc que \( \pi r^2 h = 1\). Je me sers de cette dernière égalité pour exprimer S en fonction de r uniquement, ce qui me donne \(S = \frac{2}{r} + 2\pi r^2 \).
Il me reste à calculer la dérivée de cette fonction, et à chercher sa racine. Le script BODeriveeOpt.py va le faire pour moi !
Après avoir défini la fonction de calcul de la surface du cylindre et le domaine de variation du rayon (attention à éliminer le rayon nul!), je calcule la dérivée de la fonction Surface et je cherche le rayon ropt pour lequel elle s'annule.
Le programme trace les deux courbes S et S', ce qui nous permet de lire le rayon optimum, qui est le même que celui qu'affiche le programme, soit 0,54 dm.
Voici les deux courbes :
Imaginez la situation suivante : vous disposez d’un tableau de nombres acquis au cours d’une expérience, qui décrit la variation d’une grandeur physique en fonction d’une autre, la tension en fonction du temps par exemple. Bien sur, vous aimeriez en déduire la loi physique de variation et donc la fonction mathématique qui relie ces nombres. Autre situation, rencontrée plus haut, vous disposez d’un tableau de nombres issus d’un calcul, et vous voudriez connaitre la fonction qui lie ces points.
Pour résoudre ce problème, nous disposons de deux méthodes : la régression, comme par exemple la méthode des moindres carrés et l’interpolation. Nous avons étudié la régression linéaire dans une autre page . Ici, nous allons voir ce qu’est l’interpolation.
Supposons donc que les couples de mesures dont nous disposons (xi,yi) soient des points représentatifs de la courbe d’une fonction f que nous cherchons à déterminer. Nous allons interpoler cette fonction f, c'est-à-dire l’approcher par un polynôme, en s’arrangeant pour que la courbe représentative de ce polynôme passe par tous les points (xi,yi) de notre tableau. Evidemment, on cherchera à ce que le polynôme d’interpolation soit le plus proche possible de la fonction f, dans le domaine de variation défini par notre tableau de mesures.
Il existe de nombreuses méthodes d’interpolation polynomiale. Vous étudierez sans doute en analyse numérique au moins les polynômes de Lagrange et ceux de l’Hermite. Il existe aussi d’autres méthodes d’interpolation, comme les splines, qui corrigent certains défauts de l’interpolation polynomiale classique. Pour tout ça, je vous renvoie à votre cours d’analyse numérique.
Le package scipy nous fournit un grand nombre de méthodes d’interpolation, plus ou moins élaborées. J’ai choisi de vous présenter ici l’une de plus courantes et néanmoins efficace pour les problèmes simples qui nous occupent. Il ‘agit de la routine interp1d().
La routine interp1d() permet de construire une fonction à partir de deux tableaux de nombres x[i] et y[i] dans le domaine de variation défini par ces tableaux. Selon le paramètre passé, elle travaille en différents types d’interpolations.
Vous trouverez la documentation de la routine interp1d() sur cette page .
Lors d’une expérience, j’ai mesuré la masse volumique de l’eau pure en fonction de la température, pour quelques températures données. Mes résultats de mesure sont :
Pour 0° C la masse volumique est de 999.87 kg/m3
Pour 5° : 999.99
Pour 10° : 999.73
Pour 15° : 999.13
Pour 20° : 998.23
A partir de ces mesures, je vais construire la fonction d’interpolation et pouvoir ainsi, par exemple, déterminer la masse volumique de l’eau à 18°.
Examinons le script qui me permet cette opération. Comme d’habitude, je commence par charger les packages utiles :
from scipy import arange, array
from scipy.interpolate import interp1d
from matplotlib.pyplot import *
Puis je stocke mes mesures dans deux tableaux qui seront les paramètres d’entrée de la routine interp1d() :
Temperature = array([0.0, 5.0, 10.0, 15.0, 20.0])
MasseVol = array([999.87, 999.99, 999.73, 999.13, 998.23])
Et j’appelle la routine d’interpolation, en lui demandant de faire une interpolation par spline cubique :
P = interp1d(Temperature,MasseVol, kind='cubic')
La routine interp1d() me retourne un polynôme dans la variable P. Je vais utiliser ce polynôme pour calculer la courbe d’interpolation pour x variant entre 0° et 20°, afin de pouvoir la tracer :
x = arange(0,20,0.1)
C = P(x)
Enfin, je trace mes points de mesures et la courbe d’interpolation :
figure()
grid(True)
plot(Temperature, MasseVol,'o', ms=8, label="Mesures")
plot(x,C,'r', label = "interpolation cubique")
legend()
show()
Voilà ce que cela donne :
Vous constatez que la courbe d’interpolation passe par tous les points de mesure, contrairement à ce que ferait une courbe de régression.
Vous pouvez aussi lire sur cette courbe, en utilisant le curseur mobile, la valeur de masse volumique de l’eau pure à 18° qui est de 998.62 kg/m3, par exemple.
Un des objets majeurs de la physique numérique est la simulation des phénomènes physiques. Dans beaucoup de cas, ces phénomènes présentent des comportements aléatoires ou du moins vus comme tels au niveau macroscopique, en physique statistique par exemple. La nécessité de créer de l'aléatoire est alors évidente.
Dans certains algorithmes de calcul, l'aléatoire est indispensable. C'est le cas de la méthode de Monte-Carlo, massivement utilisée pour le calcul d'intégrales et de physique statistique.
Le physicien numéricien se retrouve vite à devoir disposer d'un générateur de nombres (pseudo-)aléatoires. Nous avons déjà abordé ce sujet dans une page, à propos de la génération de nombres pseudo aléatoires. je vous y renvoie pour plus de détails.
Comme presque tous les langages de programmation, Python dispose d'un générateur de nombres float pseudo-aléatoires répartis selon une loi qui peut être choisie (uniforme, gaussienne ou autre) dans l'intervalle [0,1]. Il est disponible dans le package random, qui n'est pas un package de scipy, mais un package standard de Python. Vous trouverez la documentation de ce package dans cette page .
Dans ma page sur les générateurs de nombres pseudo-aléatoires , je propose un programme en C de calcul de Pi, à partir d'une méthode de tir statistique. Je vous propose ici le même programme en Python, afin que vous puissiez comparer la complexité du code et son efficacité.
Le script Python débute comme d'habitude par l'import des librairies :
from random import random
Puis je construis un tableau de N points (x,y) générés aléatoirement. Notez que je n'ai pas défini la "graine", le seed(), comme je l'ai fait en C. Comme vous le dira la doc du package, Python utilise par défaut l'horloge système comme "graine". c'est d'ailleurs ce que je fais dans le programme C. Comme en C, on peut choisir une "graine" particulière en utilisant random.seed().
N = 10**7
points = [(random(),random()) for _ in range(N)]
Puis vient une opération bien pratique et très puissante permise par Python : j'utilise la fonction filter() pour éliminer du tableau points tous les couples (x,y) qui ne répondent pas au critère de distance d <= (x^2 + y^2). Vous l'avez compris, il s'agit de tous les points qui ne sont pas à l'intérieur ou sur le bord du disque de rayon 1. Notez l'usage du mot clé lambda, que j'utilise rarement. Il permet de déclarer une fonction en ligne, ici le calcul de la distance.
disque = filter(lambda(x,y) : x**2 + y**2 <= 1, points)
Enfin, je calcule la valeur de PI, sachant que la probabilité qu'un point (x,y) soit compris dans le disque est de PI/4, et que celle-ci s'exprime le rapport entre le nombre de points de disque et le nombre de points du tableau points :
PI = 4.0*len(disque)/float(N)
print 'Valeur estimee de PI : ',PI
Attention à ne pas utiliser la variable pi mais PI ! En exécutant le script, vous constaterez que la valeur de PI affichée n'est pas trop mauvaise, à condition que N soit suffisamment grand. Vous noterez aussi que le délai d'exécution du script est assez long...
Si vous faites tourner les deux codes C et Python, en comparant les résultats et le temps de calcul, vous ne serez pas surpris de constater que le code C est plus rapide pour le même nombre de points. Les compilateurs C sont bien plus performants que l'interpréteur Python, mais ça on le savait déjà !
J'ai déjà abordé la simulation du mouvement brownien dans cette page. Vous pourrez vous y référer pour plus de détails. Dans cette page, je propose un programme Scilab pour simuler une marche aléatoire 2D bruitée, soit un mouvement brownien.
Je vous propose ici le script BOWalk2D.py qui simule une marche aléatoire 2D, qui mettra en oeuvre les fonctions disponibles de Python, en s'inspirant du même algorithme que le programme SciLab.
Après les traditionnels chargements de librairie et les déclarations des constantes de la simulation, je construis les tableaux x[] et y[] qui contiendront les coordonnées de la trajectoire:
x = zeros(N)
y = zeros(N)
Puis je construis les tableaux bruitX[] et bruitY[], qui contiennent les valeurs aléatoires. Notez que ces valeurs sont distribuées selon une loi gaussienne (normale) de moyenne 0 et d'écart-type sqrt(T).
bruitX = [gauss(0.,sqrt(T)) for i in range(N)]
bruitY = [gauss(0.,sqrt(T)) for i in range(N)]
Enfin, la boucle de calcul de la trajectoire :
for i in range(1,N):
x[i] = x[i-1] + sh*bruitX[i-1]
y[i] = y[i-1] + sh*bruitY[i-1]
L'affichage de la trajectoire n'appelle aucun commentaire. Ah si, peut-être l'instruction axis('equal') qui permet de dimensionner x et y de manière identique.
Avec les valeurs du code, nous obtenons :
N'hésitez pas à modifier le nombre N de pas et la variable T.
Nous avons abordé de nombreuses fois les équations différentielles ordinaires, les EDO, sur ce site. Elles sont un outil fondamental pour le physicien. Aussi, je vais vous renvoyer à ces multiples pages, et plus spécifiquement à cette page. Ici, je vous rappelerai brièvement comment résoudre une EDO de premier ordre, une EDO de second ordre et un système d'EDO avec Python.
Pour les besoins courants, Python propose la routine odeint() du package scipy.integrate. Vous en trouverez la documentation ici. Elle devrait convenir à tous nos besoins.
Si vous parcourez les pages de TangenteX.com, vous trouverez beaucoup d'EDO de premier ordre, chaque fois que l'on étudie la variation d'une grandeur physique en fonction d'une autre (souvent le temps). J'ai donc cherché une EDO qui change un peu de ce que l'on rencontre ici. Et je me suis souvenu d'une classe particulière d'EDO, les EDO périodiques. Elles sont de la forme x' = f(x,t) et on les rencontre dans les études sur les systèmes dynamiques, une grande spécialité d'un grand ancien, Henri Poincaré.
Leur mathématique n'est pas extrêmement compliquée mais sort du cadre de TangenteX.com. Si cela vous intéresse, je vous recommande l'excellent ouvrage "Equations différentielles et systèmes dynamiques" de John Hubbard et Beverly West, au §5.5.
Voici le script BOedo1.py, qui trace la courbe intégrale (la solution) d'une EDO périodique \( \frac{dx}{dt} = \cos(x^2 + t) \) avec x0 = 0. Ces grandes lignes :
L'importation des librairies:
from scipy import cos, arange,pi
from scipy.integrate import odeint
from matplotlib.pyplot import *
La définition de la fonction dérivée fprime, qui est le bout de code le plus important :
def fprime(x,t):
x_point = cos(x**2 + t)
return x_point
La définition du vecteur temps, rien de bien original :
t0 = 0.
tmax = 10*pi
pastemps = 0.1
time = arange(t0, tmax, pastemps)
Et l'appel à odeint(), après avoir défini la condition initiale:
x0 = 0
N = odeint(fprime,x0,time)
Je vous fais grâce du code d'affichage de la courbe intégrale, plus standard, tu meurs!
Voilà ce que cela donne :
Les systèmes dynamiques, linéaires ou non, constituent un chapitre essentiel de la physique mathématique. Leur étude est passionnante et très enrichissante. Je prévois d'ailleurs d'y consacrer une page particulière. Nous avons déjà examiné un de ces systèmes dans la page sur le modèle de Lotka-Volterra.
Ici, je vais juste vous donner un exemple de code pour traiter ces systèmes. Il vous suffira de l'adapter pour étudier d'autres systèmes.
J'ai choisi un système relativement simple, composé de deux équations différentielles \(\frac{dx}{dt} = -y + x(1 - x^2 - y^2) \) et \(\frac{dy}{dt} = x + y(1 - x^2 - y^2) \).
Le script BOedo2.py est très semblable au script BOedo1.py dont il reprend les grandes lignes. Bien sur, la définition de la fonction à intégrer est différente, la voici :
def EDO(sysEDO,t):
x = sysEDO[0]
y = sysEDO[1]
dx = -y + x*(1 - x**2 - y**2)
dy = x + y*(1 - x**2 - y**2)
return [dx, dy]
Vous retrouvez ici la définition de notre système différentiel.
Le code de déclaration des conditions initiales x0 et Y0 est aussi différent, car nous avons ici deux conditions initiales:
x0 = 0.01
y0 = 0.01
C0 = array([x0, y0])
Je vous invite vivement à tenter l'expérience avec des conditions initiales différentes. Vous noterez le comportement du système dynamique...
Pour les conditions initiales choisies, voici la courbe intégrale obtenue :
Vous notez que la trajectoire rejoint rapidement le cercle d'équation \( x^2 + y^2 = 1 \). C'est un attracteur du système dynamique. Nous aurons l'occasion d'aborder ce concept dans une page à venir.
La résolution d'une EDO de second ordre passe par la décomposition de l'EDO en deux EDO de premier ordre et donc par la résolution d'un système d'EDO de premier, ce qui nous ramène au problème précédent (comme on aime le faire à l'X !).
Pour illustrer cet exercice, j'ai choisi le classique modèle de l'oscillateur, un ressort avec une masse qui oscille. Mais pour agrémenter l'exercice, que vous avez fait 100 fois, j'ai décidé que mon ressort n'obéirait pas à la loi de Hook, c'est à dire que la force de rappel n'est égale à -kx mais à -x*sin(x/k). Pourquoi pas ! Voyons ce que cela donne...
Le script BOedo3.py commence par l'importation des librairies. je ne reviens pas sur ces lignes.
Puis je définis la loi de rappel de mon ressort:
def Rappel(x):
k = 4.0
y = -x*sin(x/k)
return y
et le système différentiel qui décrit mon EDO :
def EDO(X,t):
x = X[0]
v = X[1]
dx = v
dv = Rappel(x)
return array([dx, dv])
Rien que du classique !
Puis viennent les lignes de définition des conditions initiales. Attention, il y a en deux, puisque nous sommes dans le cas d'une EDO de second ordre. Je ne reviens pas dessus, c'est aussi du très classique.
Enfin j'intègre mon EDO avec odeint() :
X = odeint(EDO,C0,time)
Et je trace la courbe intégrale :
Notez la forme particulière de la courbe, qui est loin d'une sinusoïde. Faites donc varier la valeur de k pour obtenir d'autres courbes..
La Transformée de Fourier Discrète, la TFD, a déjà fait l'objet de plusieurs pages sur TangenteX, par exemple ici ou ici. Je ne vais donc pas m'étendre sur le sujet et je vous invite à lire ces pages.
Ici, je vais rappeler les outils Python nécessaires pour calculer la TFD d'un signal échantillonné et tracer son spectre d'amplitude.
Nous retrouverons deux fonctions que vous connaissez du package scipy.fftpack : fft() et fftfreq(). En fait seul l'usage de la première est indispensable, fft() permettant le calcul de la TFD par l'algorithme de FFT. La seconde, fftfreq(), n'est qu'une routine utilitaire qui produit le vecteur des fréquences de décomposition du spectre. Il est tout à fait possible de le produire à la main, et même assez facilement, mais bon cela nous évite quelques lignes de programmations.
Dans le code source du script, vous rencontrerez une nouveauté, tirée du package matplotlib.mlab, la fonction magnitude_spectrum(), qui en une seule ligne calcule le spectre d'amplitude d'un signal.
Je vais rester très classique en traçant le spectre d'amplitude d'un signal composé par une somme de deux sinusoïdes d'amplitude et de fréquence différentes.
Le script BOfft.py est très classique. Il commence par l'importation des librairies, parmi lesquelles :
from scipy.fftpack import fft,fftfreq
from matplotlib.mlab import magnitude_spectrum
Je définis ensuite la routine SpectreAmp() qui calcule la TFD du signal et trace son spectre:
# Fonction de tracé du spectre d'amplitude
# s : vecteur signal
# Te : période d'échantillonnage
def SpectreAmp(s,Te):
N = s.size
# vecteur fréquences
freq = fftfreq(N,Te)
freq = freq[range(N/2)]
# calcul de la FFT normalisée et extraction de la partie réelle positive
S = abs(fft(s))/N
S = S[range(N/2)]
# affichage du spectre d'amplitude
plot(freq,S,'ro')
xlabel('Frequence (Hz)')
ylabel('Amplitude')
return
Puis la routine MagSpectrum(), qui utilise la routine magnitude_spectrum() du package mlab :
def MagSpectrum(s,Fe):
S,F = magnitude_spectrum(s,Fe,sides='onesided',pad_to=512)
plot(F,S,'r')
xlabel('Frequence (Hz)')
ylabel('Amplitude')
Vient ensuite la routine qui permet de construire le vecteur signal :
def Signal(t):
A0 = 1.
A1 = 0.5
f0 = 5.
f1 = 20.
K = 2.*pi
s = A0*sin(K*f0*t) + A1*sin(K*f1*t)
return s
Vous pouvez bien sur modifier l'amplitude et la fréquence de chaque harmonique, mais faites attention à la fréquence d'échantillonnage, qui est définie ici :
Fe = 200
Te = 1./Fe
Comme d'habitude, je construis le vecteur temps :
t0 = 0.
tmax = 1.
t = arange(t0,tmax,Te)
Puis le vecteur s contenant le signal :
s = Signal(t)
J'affiche le signal dans la première partie de la fenêtre graphique (c'est l'objet du subplot()) :
subplot(2,1,1)
plot(t,s, 'b')
xlim(t0,tmax)
ylim(s.min(),s.max())
xlabel('Temps (s)')
ylabel('Amplitude')
Et le spectre d'amplitude dans la deuxième partie de la fenêtre. Vous pouvez l'afficher soit avec la routine SpectreAmp(), soit avec MagSpectrum():
subplot(2,1,2)
#SpectreAmp(s,Te)
MagSpectrum(s,Fe)
L'affichage obtenu est :
On retrouve bien sur notre spectre d'amplitude les deux pics à 5 et 20 Hz. Leur amplitude relative est bien également respectée.
Les problèmes aux valeurs propres (eigenvalues problem) sont très courants en physique. On les rencontre en mécanique, classique et quantique, en électromagnétisme et en physique de la matière condensée. Ils sont très généralement traités numériquement et constituent un chapitre important de la physique numérique.
Il n'est peut-être pas inutile de rappeler la notion de vecteur propre et valeur propre en algébre linéaire, et plus précisément ici en calcul matriciel.
Imaginons une équation matricielle de la forme Av = \( \lambda \)v, avec A une matrice carrée d'ordre N à coefficients réels ou complexes sans caractéristique particulière. Si l'équation est satisfaite, on appelle v vecteur propre de la matrice et \( \lambda \), un scalaire réel ou complexe, valeur propre. On en déduit presque immédiatement qu'il est équivalent à cette équation de poser det(A - \( \lambda \)I) = 0; I étant bien sur la matrice identité.
Le problème auquel on est confronté est de chercher les valeurs propres et vecteurs propres d'une matrice. Je n'ai pas la prétention de vous faire un cours de calcul matriciel ici, je vais donc passer directement aux moyens numériques offerts par Python pour résoudre ce problème, qui peut être très ... désagréable ! Si vous souhaitez plus de détails théoriques, je vous renvoie à votre cours d'algèbre.
Le package scipy.linalg fournit plusieurs fonctions de calcul des vecteurs propres et valeurs propres d'une matrice, selon la nature de la matrice. La page suivante vous éclairera sur ce choix.
Pour notre exemple, je choisiserai la fonction eig(), qui calcule les valeurs et vecteurs propres d'une matrice carrée. Vous trouverez sa documentation sur la page suivante.
Son emploi est simple : on lui passe en paramètre la matrice que l'on veut traiter (entre autres, mais pour nous cela suffira), et elle retourne deux tableaux (entre autres...) contenant les valeurs propres et les vecteurs propres de la matrice. Voyons cela sur un exemple.
Pour illustrer l'usage des valeurs propres et vecteurs propres, nous allons nous intéresser à un problème que nous retrouverons un autre jour : la vibration d'élogation des atomes dans une molécule diatomique. Imaginons une molécule de chlorure d'hydrogène HCl et essayons de calculer ses modes de vibration. Vous avez rencontré ,de très loin, ce problème en TS dans l'étude du spectre infrarouge d'une molécule...
Dans le modèle simplifié de mécanique classique (pas d'équation de Schrödinger ici !), nous allons considérer que nous avons affaire à deux atomes ponctuels, l'un de chlore de masse m1 et un d'hydrogène de masse m2, reliés par un ressort de raideur k, qui modélise la constante de liaison entre les deux atomes. Nous les faisons vibrer autour d'un point d'équilibre x0 sur un axe horizontal. Pour simplifier, on traitera le problème dans une seule dimension.
En appliquant le PFD dans le référentiel ainsi constitué, on obtient le système d'équations différentielles couplées suivant:
\( \left\{ \begin{array}{r c l} \frac{d^2x_1}{dt^2} &=& -\frac{k}{m_1} x_1 + \frac{k}{m_1} x_2 \\ \frac{d^2x_2}{dt^2} &=& \frac{k}{m_2} x_1 - \frac{k}{m_2} x_2 \end{array} \right. \)
On peut écrire ce système d'EDO sous forme matricielle : \( \frac{d^2 \vec{X}}{dt^2} = A\vec{X} \) avec :
\( A = \begin{pmatrix} -\frac{k}{m_1} & \frac{k}{m_1} \\ \frac{k}{m_2} & -\frac{k}{m_2} \end{pmatrix} \) et \( \vec{X} = \begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} \)
D'après la physique de l'oscillateur, nous recherchons des solutions de la forme \( \vec{X} = \vec{X_0}\exp(i \omega t) \). En dérivant deux fois cette équation, nous obtenons \( \omega^2 \vec{X_0} = A\vec{X_0} \), forme qui devrait vous rappeler quelque chose si je pose \( \lambda = \omega^2 \). Ce qui nous conduit à calculer les valeurs propres de la matrice A !
Il s'agit d'une matrice carrée d'ordre 2, donc nous devons obtenir deux valeurs propres et deux vecteurs propres. Le script BOEigen.py va faire ce calcul pour nous. Il est très simple. Nous débutons par l'appel des librairies, ici une seule:
from scipy.linalg import eig
Puis je définis les paramètres de mon système , je les ai choisi plus ou moins au hasard, surtout k, ce n'est qu'un exemple après tout :
k = 10. # raideur ressort
m1 = 1. # masse H
m2 = 35. # masse Cl
A partir de ces paramètres, je déclare la matrice A comme définie plus haut :
A = [[-k/m1,k/m1],[k/m2,-k/m2]]
Puis j'appelle la fonction de calcul eig()
la,v = eig(A)
Elle retourne deux tableaux: la contenant les deux valeurs propres et v contenant les vecteurs propres, tableaux que j'affiche :
l1,l2 = la
print 'Les valeurs propres sont :', l1, 'et', l2
print 'premier vecteur propre : ',v[:,0]
print 'second vecteur propre : ',v[:,1]
Le programme retourne les résultats suivants :
Les valeurs propres sont : (-10.2857142857+0j) et 0j
premier vecteur propre : [-0.99959209 0.02855977]
second vecteur propre : [-0.70710678 -0.70710678]
On retrouve la valeur propre triviale et nulle puis l'autre valeur propre égale à \( \lambda_2 = -\frac{k}{\mu} \) avec \( \mu = \frac{m_1 m_2}{m_1 + m_2} \), la masse réduite du système. Le système a une seule valeur propre non nulle et donc un seul mode propre de vibration, qui est proportionnel à une sinusoïde de pulsation propre \( \omega = \sqrt(\frac{k}{\mu}) \).
Les scripts Python étudiés dans cette page sont disponibles dans le package BOPython.rar :
J'espère que les différents outils Python que nous venons d'aborder, même très simples, vous rendront service dans vos expériences de physique numérique.
J'espère aussi qu'ils vous donneront envie de creuser plus loin les différents aspects du calcul numérique que nous venons d'effleurer. Le domaine est immense et plein de surprises. D'ailleurs, ce serait sympa de nous faire partager vos découvertes réalisées au fil de vos TIPE et autres TPE.
Contenu et design par Dominique Lefebvre - tangenteX.com mars 2015 -- 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 .