Tableaux et matrices¶
L'un des points communs fondamentaux de Numpy et de Scipy est leur aptitude à gérer les tableaux et les matrices. Vous pouvez créer des tableaux à partir de listes Python (list ou []) grâce à la fonction objet array() :
array([1,2,3,4,5,6])
Vous pouvez passer en second argument de la fonction array() le type numérique souhaité. Vous trouverez une liste exhaustive de ces types à cet endroit. Certains types ont un équivalent qui se résume en une seule lettre. Les plus courants sont 'd' (nombre à virgule flottante à double précision), 'D' (nombre complexe à double précision) et 'i' (int32, entier signé codé sur 32 bits, valeurs allant de -2 147 483 648 à +2 147 483 647).
Exemples :
array([1,2,3,4,5,6],'d')
array([1,2,3,4,5,6],'D')
array([1,2,3,4,5,6],'i')
Pour créer une matrice, vous pouvez utiliser array() avec des listes de listes Python :
array([[0,1],[1,0]],'d')
Vous pouvez aussi générer des matrices vides (remplies de zéros) de tailles arbitraires, y compris des tableaux unidimensionnels, que Numpy traitera comme des tableaux à une seule ligne grâce à la fonction objet zeros() :
zeros((3,3),'d')
Le premier argument de la fonction est un tuple (m_lignes, n_colonnes) et le second argument définit le type de données à utiliser, comme pour la fonction objet array(). Ainsi, vous pouvez initialiser un tableau simple (une ligne) de cette façon :
zeros(3,'d')
zeros((1,3),'d')
ou des tableaux en colonne :
zeros((3,1),'d')
Il y a aussi la fonction objet identity() qui produit, comme vous pourriez vous y attendre :
identity(4,'d')
de même que la fonction objet ones(), qui fait la même chose que zeros(), sauf qu'elle initialise la matrice avec des uns (1).
Espace linéaire, fonctions matricielles et tracés¶
La fonction objet linspace(start, end) crée un tableau espace linéaire de points de valeurs comprises entre start (inclus) et end (inclus).
linspace(0,1)
Si vous fournissez un troisième argument à linspace(start, end, nb_points), celui-ci définira le nombre de points à calculer dans cet espace linéaire. Si vous omettez ce troisième argument, il prend par défaut la valeur nb_points=50.
linspace(0,1,11)
La fonction objet linspace() est un moyen bien pratique de générer des valeurs de coordonnées pour un graphique (NdT : pour un intervalle d'abscisses, par exemple). Les fonctions de la librairie Numpy – qui sont toutes intégrées d'office dans l'interface IPython Notebook – peuvent traiter un tableau entier (ou même une matrice) de points en une seule fois.
Exemple :
x = linspace(0,2*pi)
sin(x)
En alliant Numpy à la puissance de Matplotlib, nous obtenons de grandes capacités pour tracer des graphiques facilement :
plot(x,sin(x))
Opérations sur les matrices¶
Les objets matrice font ce qu'il faut lorsqu'ils sont multipliés par des scalaires :
0.125*identity(3,'d')
de même lorsque vous additionnez deux matrices entre elles (sous réserve toutefois qu'elles soient de mêmes dimensions).
identity(2,'d') + array([[1,1],[1,2]])
Une chose qui pourrait dérouter les utilisateurs de Matlab™ : l'opérateur arithmétique dénoté par l'astérisque (*) applique une multiplication élément par élément plutôt qu'une vraie multiplication entre matrices :
identity(2)*ones((2,2))
Pour obtenir une vraie multiplication entre matrices, utilisez la fonction dot() :
dot(identity(2),ones((2,2)))
la fonction dot() peut aussi faire des produits scalaires :
v = array([3,4],'d')
sqrt(dot(v,v))
tout comme la multiplication matrice-tableau.
Les fonctions determinant(), inverse() et transpose() produisent les résultats attendus. Une transposition peut être abrégée en plaçant .T à la fin d'un nom d'objet matrice :
m = array([[1,2],[3,4]])
m.T
Vous avez aussi la fonction diag() qui place une liste ou un tableau le long de la diagonale d'une matrice carrée.
diag([1,2,3,4,5])
Nous verrons son utilité plus tard.
Solveurs matriciels¶
Vous pouvez résoudre des systèmes d'équations linéaires avec la fonction solve() :
A = array([[1,1,1],[0,2,5],[2,5,-1]])
b = array([6,-4,27])
solve(A,b)
Plusieurs fonctions pour calculer des valeurs propres ainsi que des vecteurs propres :
eigvals()retourne les valeurs propres d'une matrice ;eigvalsh()retourne les valeurs propres d'une matrice hermitienne ;eig()retourne les valeurs propres et les vecteurs propres d'une matrice ;eigh()retourne les valeurs propres et les vecteurs propres d'une matrice hermitienne.
A = array([[13,-4],[-4,7]],'d')
eigvalsh(A)
eigh(A)
Exemple : différences finies¶
À présent que nous avons tous ces outils dans notre boîte à outils, nous pouvons commencer à faire des choses intéressantes. La plupart des équations que nous cherchons à résoudre en Physique implique les équations différentielles. Nous aimerions pouvoir calculer les dérivées de fonctions \(y'(x) =\lim\limits_{\substack{h \to 0 \\ h\neq 0}} \frac{y(x+h)-y(x)}{h}\) en discrétisant la fonction \(y(x)\) sur un espace de points \(x_0, x_1, \dots, x_n\) répartis uniformément qui produiront \(y_0, y_1, \dots, y_n\). Grâce à la discrétisation des valeurs, nous pouvons approcher la dérivée \(y'(x)\) par \(y_i' \approx \frac{y_{i+1}-y_{i-1}}{x_{i+1}-x_{i-1}}\)
En Python, nous pourrions écrire une fonction dérivée comme ceci :
def nderiv(y,x):
"Différence finie, dérivée de la fonction f"
n = len(y)
d = zeros(n,'d') # virgule flottante à double précision (double)
# différences de part et d'autre
# centrées sur les points intérieurs
for i in range(1,n-1):
d[i] = (y[i+1]-y[i-1])/(x[i+1]-x[i-1])
# différences sur un seul côté pour les extrémités
d[0] = (y[1]-y[0])/(x[1]-x[0])
d[n-1] = (y[n-1]-y[n-2])/(x[n-1]-x[n-2])
return d
Voyons si cela fonctionne avec l'exemple sin(x) précédent :
x = linspace(0,2*pi,15)
dsin = nderiv(sin(x),x)
plot(x,dsin,label='numerical')
plot(x,cos(x),label='analytical')
title("Comparison of numerical and analytical derivatives of sin(x)")
legend()
Très proches ! (NdT : avec seulement 15 points pour discrétiser)
Oscillateur harmonique unidimensionnel utilisant les différences finies¶
Voyons si nous pouvons utiliser ce schéma aux différences finies pour calculer l'oscillateur harmonique unidimensionnel avec une bonne approximation.
Nous aimerions résoudre l'équation indépendante du temps de Schrödinger :
\(-\frac{\hbar^2}{2m}\frac{\partial^2\psi(x)}{\partial x^2} + V(x)\psi(x) = E\psi(x)\)
avec \(\psi(x)\) quand \(V(x)=\frac{1}{2}m\omega^2x^2\) représente le potentiel de l'oscillateur harmonique. Nous utiliserons le tour de passe-passe habituel, qui consiste à transformer l'équation différentielle en équation matricielle, en multipliant chaque terme de l'équation par l'expression conjuguée \(\psi^*(x)\) puis en intégrant sur \(x\). Ce qui nous donne :
\(-\frac{\hbar}{2m}\int\psi^*(x)\frac{\partial^2}{\partial x^2}\psi(x)\mathrm{d}x + \int\psi^*(x)V(x)\psi(x)\mathrm{d}x = E\)
NdlR : E étant une constante on peut la sortir de l'intégrale, reste donc pour l'expression de droite \(E\int \psi^*(x)\psi(x)\mathrm{d}x\) or, \(\psi^*(x)\psi(x)=\parallel\psi(x)^2\parallel\) est la densité de probabilité de présence. Celle-ci, étant intégrée sur l'ensemble du domaine, vaut 1. Donc \(\int \psi^*(x)\psi(x)\mathrm{d}x= \int\parallel\psi(x)^2\parallel=1\), d'où la simplification de l'expression de droite à \(E\).
Nous allons de nouveau recourir à l'approche par différences finies. La formule de différences finies pour la dérivée seconde est :
\(y''_{i} = \frac{y_{i+1}-2y_i+y_{i-1}}{x_{i+1}-x_{i-1}}\)
Nous pourrions penser au premier terme de l'équation de Schrödinger comme étant le chevauchement de la fonction d'onde \(\psi(x)\) avec la dérivée seconde de la fonction d'onde \(\frac{\partial^2}{\partial x^2}\psi(x)\). Avec la formule de la dérivée seconde que nous venons d'écrire, nous pouvons voir – si nous prenons le chevauchement des états \(y_1,\dots,y_n\) avec la dérivée seconde – que nous aurons seulement trois points où le chevauchement est non nul, à \(y_{i-1}\), \(y_{i}\) et \(y_{i+1}\). Sous la forme matricielle, cela nous mène à la matrice tridiagonale de Laplace, qui a pour valeurs des -2 dans la diagonale principale et des 1 dans les diagonales situées au-dessus et en-dessous de cette diagonale principale.
Le second terme de l'équation de Schrödinger mène à une matrice diagonale avec pour éléments diagonaux \(V(x_i)\). En assemblant toutes les pièces de ce puzzle, nous obtenons :
def Laplacian(x):
h = x[1]-x[0] # on suppose les points répartis uniformément
n = len(x)
M = -2*identity(n,'d')
for i in range(1,n):
M[i,i-1] = M[i-1,i] = 1
return M/h**2
x = linspace(-3,3)
m = 1.0
ohm = 1.0
T = (-0.5/m)*Laplacian(x)
V = 0.5*(ohm**2)*(x**2)
H = T + diag(V)
E,U = eigh(H)
h = x[1]-x[0]
# on trace le potentiel harmonique
plot(x,V,color='k')
for i in range(4):
# Pour chacune des premières solutions, on trace le niveau d'énergie
axhline(y=E[i],color='k',ls=":")
# de même que les valeurs propres, décalées par le niveau d'énergie
# de sorte qu'elles ne s'empilent pas les unes sur les autres :
plot(x,-U[:,i]/sqrt(h)+E[i])
title("Eigenfunctions of the Quantum Harmonic Oscillator")
xlabel("Displacement (bohr)")
ylabel("Energy (hartree)")
Nous avons bricolé un peu pour obtenir les orbitales telles qu'attendues. Tout d'abord, nous avons inséré un facteur -1 devant les fonctions d'onde pour régler la phase à son état le plus bas. La phase (signe) d'une fonction d'onde quantique ne fournit aucune information, seul le carré de la fonction d'onde le fait, donc tout ceci est sans grande conséquence.
En revanche, les fonctions propres telles que nous les avons générées ne sont pas correctement normalisées. La raison en est que la différence finie n'est pas un véritable principe de base au sens de la mécanique quantique. C'est le principe de base des fonctions delta δ de Dirac à chaque point ; nous interprétons les espaces entre les points comme étant « remplis » par la fonction d'onde, mais le principe de la différence finie n'a de solution que pour les points eux-mêmes. Nous pouvons résoudre ce problème en divisant les fonctions propres de notre différence finie hamiltonienne par la racine carrée de l'espacement entre les points et ainsi nous obtenons des fonctions correctement normalisées.
Fonctions spéciales¶
Les solutions de l'oscillateur harmonique sont supposées être des polynômes hermitiens. La page Wikipedia a les états HO donnés par :
\(\psi_n(x) = \frac{1}{\sqrt{2^n n!}} \left(\frac{m\omega}{\pi\hbar}\right)^{1/4} \exp\left(-\frac{m\omega x^2}{2\hbar}\right) H_n\left(\sqrt{\frac{m\omega}{\hbar}}x\right)\)
Voyons s'ils ressemblent à cela. Il existe des fonctions spéciales dans la librairie Numpy et quelques autres encore plus spéciales dans Scipy. Les polynômes hermitiens se trouvent dans Numpy :
from numpy.polynomial.hermite import Hermite
from math import factorial
def ho_evec(x,n,m,ohm):
vec = [0]*9
vec[n] = 1
Hn = Hermite(vec)
return (1/sqrt(2**n*factorial(n)))*pow(m*ohm/pi,0.25)*exp(-0.5*m*ohm*x**2)*Hn(x*sqrt(m*ohm))
Comparons notre première fonction à notre solution.
plot(x,ho_evec(x,0,1,1),label="Analytic")
plot(x,-U[:,0]/sqrt(h),label="Numeric")
xlabel('x (bohr)')
ylabel(r'$\psi(x)$')
title("Comparison of numeric and analytic solutions to the Harmonic Oscillator")
legend()
L'accord est presque parfait.
Nous pouvons utiliser la fonction matplotlib subplot(rows, cols, index) pour afficher en grille plusieurs tracés dans un même graphique.
phase_correction = [-1,1,1,-1,-1,1]
for i in range(6):
subplot(2,3,i+1)
plot(x,ho_evec(x,i,1,1),label="Analytic")
plot(x,phase_correction[i]*U[:,i]/sqrt(h),label="Numeric")
Hormis les erreurs de phase – que j'ai corrigées par un petit bidouillage, vous voyez lequel ? – la concordance entre les courbes est plutôt pas mal, bien que cela empire au fur et à mesure que le niveau d'énergie augmente, en partie parce que nous avons utilisé une résolution de 50 points seulement.
Le module Scipy a encore plus de fonctions spéciales :
from scipy.special import airy,jn,eval_chebyt,eval_legendre
subplot(2,2,1)
x = linspace(-1,1)
Ai,Aip,Bi,Bip = airy(x)
plot(x,Ai)
plot(x,Aip)
plot(x,Bi)
plot(x,Bip)
title("Airy functions")
subplot(2,2,2)
x = linspace(0,10)
for i in range(4):
plot(x,jn(i,x))
title("Bessel functions")
subplot(2,2,3)
x = linspace(-1,1)
for i in range(6):
plot(x,eval_chebyt(i,x))
title("Chebyshev polynomials of the first kind")
subplot(2,2,4)
x = linspace(-1,1)
for i in range(6):
plot(x,eval_legendre(i,x))
title("Legendre polynomials")
Par exemple les polynômes de Jacobi, Laguerre, Hermite, les fonctions hypergéométriques et plein, plein d'autres. Vous trouverez une liste complète sur la page fonctions spéciales Scipy.
Résolution par les moindres carrés¶
Très souvent, nous traitons des données que nous voulons faire cadrer avec un comportement attendu. Supposons que nous ayons ceci :
raw_data = """\
3.1905781584582433,0.028208609537968457
4.346895074946466,0.007160804747670053
5.374732334047101,0.0046962988461934805
8.201284796573875,0.0004614473299618756
10.899357601713055,0.00005038370219939726
16.295503211991434,4.377451812785309e-7
21.82012847965739,3.0799922117601088e-9
32.48394004282656,1.524776208284536e-13
43.53319057815846,5.5012073588707224e-18"""
Vous trouverez plus tard une rubrique qui traite de l'analyse de données au format CSV (Comma-Separated Values). Nous récupérerons l'analyseur qui s'y trouve. Pour une explication détaillée, passez la rubrique ici présente et lisez plus loin. Sinon, considérez simplement qu'il s'agit d'un moyen d'extraire un tableau Numpy de la chaîne de caractères raw_data, tableau que nous pourrons tracer en graphique et soumettre à diverses analyses par la suite.
data = []
for line in raw_data.splitlines():
words = line.split(',')
data.append(map(float,words))
data = array(data)
title("Raw Data")
xlabel("Distance")
plot(data[:,0],data[:,1],'bo')
Comme nous nous attendons à ce que nos données décroissent exponentiellement, nous pouvons utiliser une échelle semi-logarithmique.
title("Raw Data")
xlabel("Distance")
semilogy(data[:,0],data[:,1],'bo')
Pour une décroissance exponentielle pure comme celle-ci, nous pouvons assimiler les points représentés à une droite. Le tracé ci-dessus suggère que nous avons là une bonne approximation, étant donné que la fonction \(y = Ae^{-ax} \Leftrightarrow \log(y) = \log(Ae^{-ax})\), \(\log(y)=\log(A)-ax\) est du type affine. Ainsi, en ajustant l'échelle par rapport aux abscisses \(x\), nous devrions obtenir une droite de coefficient directeur \(-a\) et d'ordonnée à l'origine \(\log(A)\).
Il existe une fonction numpy nommée polyfit() qui ajuste les données pour cadrer avec une forme polynomiale. Utilisons-la pour obtenir une droite (polynôme d'ordre 1) :
params = polyfit(data[:,0],log(data[:,1]),1)
a = params[0]
A = exp(params[1])
Voyons à présent si cette courbe s'ajuste aux données que nous avons :
x = linspace(1,45)
title("Raw Data")
xlabel("Distance")
semilogy(data[:,0],data[:,1],'bo')
semilogy(x,A*exp(a*x),'b-')
Avec des fonctions plus compliquées, il est possible que l'on n'arrive pas à s'ajuster sur un simple polynôme. Par exemple, avec les données suivantes :
gauss_data = """\
-0.9902286902286903,1.4065274110372852e-19
-0.7566104566104566,2.2504438576596563e-18
-0.5117810117810118,1.9459459459459454
-0.31887271887271884,10.621621621621626
-0.250997150997151,15.891891891891893
-0.1463309463309464,23.756756756756754
-0.07267267267267263,28.135135135135133
-0.04426734426734419,29.02702702702703
-0.0015939015939017698,29.675675675675677
0.04689304689304685,29.10810810810811
0.0840994840994842,27.324324324324326
0.1700546700546699,22.216216216216214
0.370878570878571,7.540540540540545
0.5338338338338338,1.621621621621618
0.722014322014322,0.08108108108108068
0.9926849926849926,-0.08108108108108646"""
data = []
for line in gauss_data.splitlines():
words = line.split(',')
data.append(map(float,words))
data = array(data)
plot(data[:,0],data[:,1],'bo')
Ces données ressemblent plus à une courbe de Gauss qu'à une courbe exponentielle. Si nous le voulions, nous pourrions utiliser la fonction polyfit() pour ce cas aussi, mais utilisons plutôt la fonction curve_fit() du module Scipy, qui peut s'ajuster à toute sorte de fonctions arbitraires. Pour en savoir plus, tapez help(curve_fit) dans une cellule de code (ou en console interactive IPython).
Tout d'abord, définissons une fonction gaussienne générique à ajuster :
def gauss(x,A,a): return A*exp(a*x**2)
À présent, ajustons en utilisant curve_fit() :
from scipy.optimize import curve_fit
params,conv = curve_fit(gauss,data[:,0],data[:,1])
x = linspace(-1,1)
plot(data[:,0],data[:,1],'bo')
A,a = params
plot(x,gauss(x,A,a),'b-')
La fonction curve_fit() que nous venons d'utiliser est basée sur les excellents algorithmes d'optimisation (minimisation) de Scipy. Pour en savoir plus, consultez la documentation Scipy.
Monte Carlo, nombres aléatoires et calcul de \(\pi\)¶
Plusieurs techniques de calcul scientifique reposent sur la méthode de Monte-Carlo, où une séquence de nombre (pseudo) aléatoires est utilisée pour approcher l'intégrale d'une fonction. Python a un bon générateur de nombres aléatoires dans sa librairie standard. La fonction random() génère des nombres pseudo-aléatoires uniformément répartis entre 0 et 1.
from random import random
rands = []
for i in range(100):
rands.append(random())
plot(rands)
La fonction random() se base sur l'algorithme de Mersenne Twister (page officielle), qui est un générateur de nombres pseudo-aléatoires particulièrement apprécié. Le module Python standard random contient aussi des fonctions pour générer des entiers aléatoires (random.randint(a, b)) des listes d'entiers aléatoires (random.randrange(n)), pour mélanger une liste (random.shuffle(mutable)) ainsi que des fonctions pour extraire des nombres aléatoires d'une distribution particulière, telle que la distribution normale, par exemple :
from random import gauss
grands = []
for i in range(100):
grands.append(gauss(0,1))
plot(grands)
Il est généralement plus efficace de générer une liste complète de nombres aléatoires en une seule opération, surtout si vous extrayez des valeurs d'une loi non uniforme. La librairie Numpy possède des fonctions pour générer des tableaux et des matrices de plusieurs lois de probabilité.
plot(rand(100))
L'un des tous premiers programmes que j'ai écrit devait servir à calculer le nombre \(\pi\) en prenant des nombres aléatoires pour les coordonnées de points (x, y) et en comptant les occurrences situées à l'intérieur du cercle unitaire. Par exemple :
npts = 5000
xs = 2*rand(npts)-1
ys = 2*rand(npts)-1
r = xs**2+ys**2
ninside = (r<1).sum()
figsize(6,6) # donne une forme carrée au graphique
title("Approximation to pi = %f" % (4*ninside/float(npts)))
plot(xs[r<1],ys[r<1],'b.')
plot(xs[r>1],ys[r>1],'r.')
figsize(8,6) # change les dimensions pour le reste du document
L'idée qui sous-tend derrière ce programme est que le ratio de l'aire du cercle unitaire divisée par l'aire du carré circonscrit est de l'ordre de \(\dfrac{\pi}{4}\). Alors, en comptant la portion de points aléatoires qui se trouvent à l'intérieur du cercle unitaire par rapport à l'ensemble du carré, nous obtenons une assez bonne estimation de la valeur de \(\pi\).
Le bout de code précédent utilise des notations Numpy de haut niveau pour calculer le rayon de chaque point de coordonnées (x, y) en une seule ligne de code, puis pour compter combien de rayons sont inférieurs à 1 – toujours en une seule ligne de code – et enfin pour filtrer les points (x, y) par rapport à leurs rayons respectifs. Pour être honnête, j'écris rarement du code sous cette forme : je trouve que ces techniques Numpy sont par trop pointues pour que l'on puisse s'en souvenir et je préfère personnellement utiliser les techniques de liste en compréhension (voir plus loin) pour extraire les points que je veux, ce qui est plus facile à mémoriser.
À comparer des méthodes actuelles de calcul de \(\pi\), celle-ci est parmi les pires. Une bien meilleure méthode consiste à utiliser la série alternée de Leibniz (développement de arctan(1)) :
\(\dfrac{\pi}{4} = \sum_k \dfrac{(-1)^k}{2\times k+1}\)
n = 100
total = 0
for k in range(n):
total += pow(-1,k)/(2*k+1.0)
print 4*total
Si vous recherchez une méthode remarquable, consultez la méthode de Srinivasa Ramanujan. Elle converge si vite que vous aurez réellement besoin d'une précision mathématique arbitraire pour afficher les décimales. Vous pouvez faire ceci avec le module Python standard decimal, si cela vous intéresse.
Intégrales¶
Le calcul numérique d'une intégrale peut s'avérer épineux et parfois, il est plus facile d'évaluer le résultat par le biais d'une approximation. Par exemple, supposons que nous voulions voir ce que donne l'intégrale :
\(\int_0^\infty\exp(-x)\mathrm{d}x=1\)
def f(x): return exp(-x)
x = linspace(0,10)
plot(x,exp(-x))
La librairie Scipy possède une fonction de calcul numérique d'intégrale nommée quad() (on appelle parfois quadrature ce genre de calculs), que nous pouvons utiliser ici :
from scipy.integrate import quad
quad(f,0,inf)
Vous avez aussi des intégrateurs 2D et 3D dans Scipy. Pour en savoir plus, consultez la documentation.
Transformée de Fourier rapide et traitement du signal¶
Nous avons très souvent recours aux techniques FFT (Fast Fourier Transform – transformée de Fourier rapide) pour nous aider à extraire un signal à partir de données « bruitées ».
from scipy.fftpack import fft,fftfreq
npts = 4000
nplot = npts/10
t = linspace(0,120,npts)
def acc(t): return 10*sin(2*pi*2.0*t) + 5*sin(2*pi*8.0*t) + 2*rand(npts)
signal = acc(t)
FFT = abs(fft(signal))
freqs = fftfreq(npts, t[1]-t[0])
subplot(211)
plot(t[:nplot], signal[:nplot])
subplot(212)
plot(freqs,20*log10(FFT),',')
show()
Dans Scipy, il existe des fonctions additionnelles qui gèrent le traitement du signal.
