Sommaire

1 Résolution de AX = Y par l’algorithme du pivot de Gauss

Formation Python ILL 2 Formalisation

Calcul Numérique avec NumPy, SciPy et Matplotlib 3 Mise en œuvre

4 Le pivot de Gauss en NumPy

Matthieu Moy
5 Analyse, comparaison avec NumPy
Ensimag
6 Vers l’infini et au delà ...

octobre 2016 7 SciPy : algorithmes de calcul numérique

8 Tracés de courbes avec matplotlib

9 Conclusion

Matthieu Moy (Ensimag) NumPy octobre 2016 < 1 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 2 / 67 >

Pivot de Gauss Exemple : solution

 
Méthode pour résoundre un système linéaire. 1 −1 2 1 | 1
Exemple : résolution dans R4 du système 0 −2 −3 −1 | 4 
 Système final (S3 ) : 
0 0 −1 −2 | 1 

x − y + 2z + t =1
0 0 0 19 | −5


2x − 4y + z + t =6

(S0 ) 
20

 x − y + z − t =2
1  −22

3x + y + z + 2t =1
 
solution : X =  
19  −9 
1 Effectuer les opérations : −5
F L2 ← L2 − 2L1
F L3 ← L3 − L1 Remarques :
F L4 ← L4 − 3L1 I Dans la ligne L1 de (S0 ), le coefficient de x est le premier pivot.
F ; Résultat = (S1 ). I Dans la ligne L2 de (S1 ), le coefficient de y est le deuxième pivot.
2 Effectuer la transformation L4 ← L4 + 2L2 pour obtenir le système I Dans la ligne L3 de (S2 ), le coefficient de z est le troisième pivot.
(S2 ). I Existence et unicité de la solution si tous les pivots sont non nuls
3 Effectuer la transformation L4 ← L4 − 11L3 pour obtenir le système (système de Cramer).
triangulaire (S3 ).
4 Résoudre dans R4 le système (S3 ).

Matthieu Moy (Ensimag) NumPy octobre 2016 < 4 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 5 / 67 >

En résumé Stabilité numérique

Attention aux arrondis !

Deux phases :
mise sous forme triangulaire >>> 12*(1/3-1/4)-1
I Lj ← Lj − λLi = « transvection ». -2.220446049250313e-16
I En cas de pivot nul, échange avec une autre ligne (il y en a toujours
au moins une dans un système de Cramer). Notion de “pivot non nul” assez fragile !
résolution du système triangulaire (phase de remontée) Si le pivot est très petit (en valeur absolue), on va obtenir des
I Si on recontre l’équation 0 = β lors de la phase de remontée, le nombres très grands lors des transvections.
système n’a aucune solution. Solution : « pivot partiel ».
I Si on rencontre l’équation 0 = 0, le système a une infinité de I On sélectionne la ligne contenant le pivot le plus grand en valeur
solutions. absolue.
I Heuristique : marche bien en pratique, mais pas toujours !

Matthieu Moy (Ensimag) NumPy octobre 2016 < 6 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 7 / 67 >

Notations Structure de données « matrice » : tableaux NumPy

Entrée : A = (aij ) est une matrice n × n composée de lignes Li

Bibliothèque très performante (code C optimisé appelable depuis
pour i de 0 à n − 2 :
Python)
I trouver j, i ≤ j ≤ n − 1, tel que |aji | soit maximal
I échanger Li et Lj Structure de base = tableau (homogène) multidimensionnel :
I pour k de i + 1 à n − 1 : [Link]
aik
F Lk ← Lk − aii
Li [Link]
Sortie : le résultat est dans A

Matthieu Moy (Ensimag) NumPy octobre 2016 < 9 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 11 / 67 >
 Construire un tableau NumPy Opérations sur tableaux numpy
On peut créer un tableau (unidimensionnel : un vecteur) NumPy à
partir d’une liste Python ou même d’une liste de listes (il sera alors
bi-dimensionel : une matrice) :
>>> import numpy NumPy permet de faire des opérations directement sur ces tableaux,
>>> [Link]([1, 2, 3]) sans avoir à réécrire chaque fonction :
array([1, 2, 3]) >>> 2 * [Link]([1, 2, 3])
>>> [Link]([[1, 2], [3, 4], [5,6]]) + [Link]([8, 10, 12])
array([[1, 2], array([10, 14, 18])
[3, 4],
; ce à quoi on s’attend en algèbre linéaire.
[5, 6]])
>>> [Link][0]
3
>>> [Link][1]
2

Matthieu Moy (Ensimag) NumPy octobre 2016 < 12 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 13 / 67 >

Accéder aux éléments d’un tableau NumPy Accès aux lignes et aux colonnes
Accès aux éléments comme pour une liste : a[i]
Numérotation à partir de 0
>>> a = [Link]([[ 0.19, 0.14, 0.21],
Matrices : M[ligne][colonne] ou M[ligne, colonne]. ... [ 0.79, 0.43, 0.74],
... [ 0.12, 0.40, 0.30]])
>>> b = [Link]([1.414, 0.5, 3.14, 2.718])
>>> b[2]
>>> a[:, 1]
3.1400000000000001
array([ 0.14, 0.43, 0.4 ])
>>> a = [Link]([[ 0.19, 0.14, 0.21],
>>> a[0, :]
... [ 0.79, 0.43, 0.74],
array([ 0.19, 0.14, 0.21])
... [ 0.12, 0.40, 0.30]])
>>> a[2, 1]
0.40000000000000002

Matthieu Moy (Ensimag) NumPy octobre 2016 < 14 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 15 / 67 >

Accès aux lignes et aux colonnes Tranches de tableaux

NumPy ne retourne pas une copie des données dans ce cas mais
seulement un vue sur celles-ci : >>> v = [Link]([-1, -2, -3, -4, -5])
>>> a = [Link]([[ 0.19, 0.14, 0.21], # elements entre les indices 2 (inclus) et 4 (exclu)
... [ 0.79, 0.43, 0.74], >>> v[2:4]
... [ 0.12, 0.40, 0.30]]) array([-3, -4])
# elements entre les indices 0 (inclus) et 3 (exclu)
>>> v = a[:, 2] >>> v[:3]
>>> v array([-1, -2, -3])
array([ 0.21, 0.74, 0.3 ]) # elements entre les indices 2 (inclus)
>>> v[1] = 1.0
# et fin (inclus)
>>> v
array([ 0.21, 1. , 0.3 ]) >>> v[2:]
>>> a array([-3, -4, -5])
array([[ 0.19, 0.14, 0.21],
; Là aussi, on ne possède que des vues sur les données et pas des
[ 0.79, 0.43, 1. ],
[ 0.12, 0.4 , 0.3 ]])
copies de celles-ci.

Matthieu Moy (Ensimag) NumPy octobre 2016 < 16 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 17 / 67 >

Exercice 1 : échange de lignes Exercice 1 : échange de lignes

Solution naive Solution avec recopie 1
Étant donnés une matrice NumPy a et deux entiers i et j, échanger
les lignes i et j de a. >>> a = [Link]([[ 0.19, 0.14, 0.21],
1ere tentative : ... [ 0.79, 0.43, 1. ],
>>> i = 0
... [ 0.12, 0.4 , 0.3 ]])
>>> j = 2 >>> tmp = a[i, :].copy() # <-- là
>>> tmp = a[i, :] >>> a[i, :] = a[j, :]
>>> a[i, :] = a[j, :] >>> a[j, :] = tmp
>>> a[j, :] = tmp >>> a
>>> a array([[ 0.12, 0.4 , 0.3 ],
array([[ 0.12, 0.4 , 0.3 ], [ 0.79, 0.43, 1. ],
[ 0.79, 0.43, 1. ],
[ 0.19, 0.14, 0.21]])
[ 0.12, 0.4 , 0.3 ]])
Oups, ne fonctionne pas ! (copie au lieu d’échange) 1. En fait, la “bonne” manière de faire en NumPy n’est pas de faire une copie mais
I tmp = a[i, :] ne fait pas de copie, juste une vue plutôt d’utiliser l’indexation avancée : a[[i, j], :] = a[[j, i], :].
I a[i, :] = a[j, :] écrase la ligne
⇒ Il faut copier explicitement

Matthieu Moy (Ensimag) NumPy octobre 2016 < 18 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 19 / 67 >
 Exercice : transvection Opérations sur les tableaux NumPy (1/2)
Étant donnés une matrice NumPy a, deux entiers i et j et un flottant
x, appliquer la transvection Lj ← Lj + xLi sur la matrice a. Déjà vu : M1 + M2, scalaire * M
>>> a Opérations classiques, composante par composante : -, *, /, abs
array([[ 0.12, 0.4 , 0.3 ],
[ 0.79, 0.43, 1. ], >>> u = [Link]([1, 2, 3, 4])
[ 0.19, 0.14, 0.21]]) >>> v = [Link]([4, 3, 2, 1])
>>> i = 0 >>> u * v
>>> j = 1 array([4, 6, 6, 4])
>>> x = -a[j, i] / a[i, i] # pivot : a[i, i] >>> u / v
>>> a[j, :] = a[j, :] + x * a[i, :] # Lj ← Lj + xLi array([ 0.25, 0.66666667, 1.5, 4.])
>>> a >>> abs(u - v)
array([[ 0.12 , 0.4 , 0.3 ], array([3, 1, 1, 3])
[ 0. , -2.20333333, -0.975 ],
[ 0.19 , 0.14 , 0.21 ]])

Matthieu Moy (Ensimag) NumPy octobre 2016 < 20 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 21 / 67 >

Opérations sur les tableaux NumPy (2/2) Exercice : Valeur maximale dans une colonne
produits matriciels, produits scalaires : [Link]
>>> u = [Link]([1, 2, 3, 4]) Étant donnée une matrice A et un indice i, écrire la fonction
>>> v = [Link]([4, 3, 2, 1]) pivot_index(A,i) qui renvoie l’indice j ≥ i t.q. |aji | est maximal.
n−1
X
>>> [Link](u, v) # ui · vi def pivot_index(A, i):
i=0 n = [Link][0] # nombre de lignes
20 j = i
>>> w = u * v for k in range(i + 1, n):
>>> [Link]() if abs(A[k, i]) > abs(A[j, i]):
20 j = k
>>> a = [Link]([[1, 2], [3, 4]]) return j
>>> x = [Link]([1, -1])
>>> b = [Link](a, x)
>>> b
array([-1, -1])

Matthieu Moy (Ensimag) NumPy octobre 2016 < 22 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 23 / 67 >

Résolution de système triangulaire supérieur Algorithme du pivot de Gauss (1/3)

Étant donnés une matrice A triangulaire supérieure, et un vecteur b de
même taille, implémenter la résolution de Ax = b. import numpy Presque tout est fait !
def solve_triangular(A, b): def pivot_index(A, i):
n = len(b) n = [Link][0] # nombre de lignes
x = [Link]([n]) # vecteur de n zeros j = i
for i in range(n - 1, -1, -1): for k in range(i + 1, n):
 
Pn−1 if abs(A[k, i]) > abs(A[j, i]):
# xi = a1 bi − j=i+1 aij xj
ii j = k
sum = [Link](A[i, :], x) return j
x[i] = (b[i] - sum) / A[i, i]
return x def swap_lines(A, i, j):
tmp = A[i, :].copy()
>>> M = [Link]([[1, 2], [0, 4]]) A[i, :] = A[j, :]
>>> x = [Link]([1, -1]) A[j, :] = tmp
>>> v = [Link](M, x) # array([-1, -4])
>>> solve_triangular(M, v) def transvection_lines(A, i, k, factor):
array([ 1., -1.]) # on retrouve x ... A[k, :] = A[k, :] + A[i, :] * factor

Matthieu Moy (Ensimag) NumPy octobre 2016 < 24 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 26 / 67 >

Algorithme du pivot de Gauss (2/3) Algorithme du pivot de Gauss (3/3 : programme

def gauss(A0, b0): principal)
A = [Link]() # pour ne pas détruire A
def solve_triangular(A, b):
b = [Link]()
n = len(b)
n = len(b) # on pourrait aussi la taille de a
x = [Link]([n])
for i in range(n - 1):
for i in range(n - 1, -1, -1):
ipiv = pivot_index(A, i)
sum = [Link](A[i, :], x)
if ipiv != i: # echanges
x[i] = (b[i] - sum) / A[i, i]
swap_lines(A, i, ipiv)
return x
tmp = b[ipiv]
b[ipiv] = b[i]
# Fonction principale
b[i] = tmp
def solve(A0, b0):
for k in range(i + 1, n): # pivotage
A, b = gauss(A0, b0)
factor = -A[k, i] / A[i, i]
return solve_triangular(A, b)
transvection_lines(A, i, k, factor)
b[k] = b[k] + b[i] * factor
return A, b

Matthieu Moy (Ensimag) NumPy octobre 2016 < 27 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 28 / 67 >
 Algorithme du pivot de Gauss : test Algorithme du pivot de Gauss : expérimentons

>>> n = 5 gauss/[Link]
# les matrices inversibles sont denses
>>> M = [Link]([n, n]) Regarder le code
>>> x = [Link]([n]) Exécuter [Link]
>>> v = [Link](M, x)
Exécuter [Link] pour tester sur une matrice
>>> x_solved = solve(M, v)
problématique (erreur attendue ≈ 33)
>>> abs(x - x_solved).max()
7.2164496600635175e-16 Si le temps le permet : modifier [Link] pour remplacer le pivot
partiel (recherche du plus petit pivot) par un pivot naif (utilisation
du premier pivot non-nul) et re-tester (erreur attendu ≈ 1016 ).

Matthieu Moy (Ensimag) NumPy octobre 2016 < 29 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 30 / 67 >

Performances Performances : comparaison avec numpy

103 Lists
Numpy arrays
102 Numpy solve

Temps de calcul (s.)

Coût en O(n3 ) (démonstration dans Wack et al. pour les curieux). 101 O(n3 )
En général, il vaut mieux faire confiance aux algorithmes 100 O(n2 )
proposés par les bibliothèques de Python : plus rapides, plus 10−1
robustes aux cas limites
10−2
10−3
10−4
10−5
100 101 102 103 104
Taille du système

Matthieu Moy (Ensimag) NumPy octobre 2016 < 32 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 33 / 67 >

Et la stabilité numérique ? Et la stabilité numérique ?

       
1 1/4 0 0 0 1 1/4 0 0 0
 1 5/4 12 0   0   1 5/4 12 0   0 
 ·x =   ·x = 
 1 1/3 1 1   1e16   1 1/3 1 1   1e16 
1 5/4 13 1 0 1 5/4 13 1 0

       
1 1/4 0 0 0 1 1/4 0 0 0
 0 5/4 − 1/4 12 0   0   0 1 12 0   0 
 ·x =    ·x =  
 0 1/3 − 1/4 1 1   1e16   0 1/3 − 1/4 1 1   1e16 
0 5/4 − 1/4 13 1 0 0 1 13 1 0

Matthieu Moy (Ensimag) NumPy octobre 2016 < 34 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 34 / 67 >

Et la stabilité numérique ? Et la stabilité numérique ?

       
1 1/4 0 0 0 1 1/4 0 0 0
 1 5/4 12 0   0   1 5/4 12 0   0 
 ·x =   ·x = 
 1 1/3 1 1   1e16   1 1/3 1 1   1e16 
1 5/4 13 1 0 1 5/4 13 1 0

       
1 1/4 0 0 0 1 1/4 0 0 0
 0 1 12 0 
  0   0 1 12 0 
  0 
 ·x =    · x =  
 0 0 1 − 12 ∗ (1/3 − 1/4) = 2.2e-16 1   1e16   0 0 2.2e-16 1   1e16 
0 1 13 1 0 0 0 1 1 0

Matthieu Moy (Ensimag) NumPy octobre 2016 < 34 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 34 / 67 >
 Et la stabilité numérique ? Vous voulez des détails ?
   
1 1/4 0 0 0
 1 5/4 12 0   0 
 ·x = 
 1 1/3 1 1   1e16 
1 5/4 13 1 0 TP6 du cours d’informatique CPP Semestre 3 :
[Link]
    document/S2/[Link]
1 1/4 0 0 0
 0 1 12 0   0 
 ·x = 
 0 0 2.2e-16 1   1e16 
0 0 0 1 − 1/2.2e-16 = -4.5e15 4.7e31

Erreur = 1016 (33 seulement avec un pivot partiel)

Matthieu Moy (Ensimag) NumPy octobre 2016 < 34 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 35 / 67 >

« Infinity » et « Not a Number » Intégration

Rb
En maths : a f (x)dx
En Python :
Exercice : quel est précisément le type des éléments de
[Link]([0.0, 1.1]) ? (indice : utiliser type) import numpy as np
import [Link] as integrate
+∞ = [Link], −∞ = - [Link].
Exercice : combien valent
result = [Link](lambda x: [Link](x),
I ∞ − ∞?
1 0, [Link])
+∞ ?
I
1 print "Result =", result[0]
−∞ ?
I
I 1.0 / 0.0 ? print "Estimated error =", result[1]
I np.float64(1.0) / np.float64(0.0) ?
# [Link] for infinity
Exercice (dur) : comment distinguer 0.0 et -0.0 ?
result = [Link](lambda x: [Link](-x**2),
-[Link], [Link])
print "Result =", result[0]
print "Estimated error =", result[1]
print "sqrt(pi) =", [Link]([Link])

Matthieu Moy (Ensimag) NumPy octobre 2016 < 37 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 39 / 67 >

Exercice Préliminaire aux équations différentielles

R +∞ Opérations sur vecteurs
Calculer 0 e−t dt.
# import ...
# result = ...
# START CUT
import numpy as np [Link](a, b, n) : vecteur de n valeurs entre a et b
import [Link] as integrate (inclus).
Exemple :
result = [Link](lambda x: [Link](-x), [Link](0, 5, 3) == array([0., 2.5, 5.]).
0, [Link])
# END CUT
print "Result =", result[0]
print "Estimated error =", result[1]

Matthieu Moy (Ensimag) NumPy octobre 2016 < 40 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 41 / 67 >

Préliminaire aux équations différentielles Généralisation : équations différentielles

Rb 0
Opérations sur vecteurs
a f (t)dt = y (b) avec y (a) = 0 et y = f (t).
Les opérations NumPy marchent composante par composante Généralisation : y 0 = f (t) ; y 0 = f (y , t).
sur les vecteurs : En Python :
>>> [Link]([Link](0, 5, 3)) # y' = f(y, t)
array([ 0. , 0.59847214, -0.95892427]) def f(y, t):
Les autres opérations, pas toujours : return ...
>>> [Link]([Link](0, 5, 3)) # Intervalle et pas d'intégration
Traceback (most recent call last): t = [Link](a, b, n)
File "<stdin>", line 1, in <module> # Valeur initiale : y(a)
TypeError: only length-1 arrays can be converted y0 = 0
to Python scalars # Résolution numérique de l'équation
Solution : y = [Link](f, y0, t)
Rπ
>>> s = [Link]([Link]) Exercice : recalculer 0 sin(x)dx avec
>>> s([Link](0, 5, 3)) [Link]
array([ 0. , 0.59847214, -0.95892427]) (Squelette : scipy/integral_sin_odeint.py).

Matthieu Moy (Ensimag) NumPy octobre 2016 < 42 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 43 / 67 >
 Tracer une courbe avec un ensemble de points Tracé de Fonction
matplotlib/courbe_simple.py matplotlib/[Link]

import [Link] as plt

import numpy as np
# plt devient l'abrege de [Link]:
import [Link] as plt # Valeurs sur lesquelles on va calculer
x = [Link](0, 10 * [Link], 1000)
x = [0, 1, 3] # liste des abscisses # application d'une fonction à chaque point
y = [1, -1, 0] # liste des ordonnees y1 = [Link](x) # [Link] s'applique sur des tableaux
y2 = [Link](-x / 10.0) * y1
[Link](x, y) # trace y en fonction de x
[Link]() # affiche la fenetre du trace # y1 en fonction de x, y2 en fonction de x
[Link](x, y1)
[Link](x, y2)

Matthieu Moy (Ensimag) NumPy octobre 2016 < 45 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 46 / 67 >

Tracé de Fonction (suite) Matplotlib et Notebook

matplotlib/[Link]

# Titres sur les axes

[Link]('Oscillateurs harmoniques libre et amorti')
[Link]('x') %matplotlib inline pour afficher inline
[Link]('y = f(x)') %matplotlib pour revenir au comportement par défaut
Cf. [Link] pour un exemple
# Tracé de l'axe y=0 et grille
[Link](y=0, color='k')
[Link]()

[Link]()

Matthieu Moy (Ensimag) NumPy octobre 2016 < 47 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 48 / 67 >

Tracé de courbes paramétrées Tracé de courbes paramétrées

matplotlib/[Link]
Jusqu’ici :
x = [Link](...)
y1 = f1(x) import [Link] as plt
y2 = f2(x) import numpy as np
[Link](x, y1)
[Link](x, y2) t = [Link](0, 10 * [Link], 1000)
x = [Link](t) * t
Courbes paramétrées :
y = [Link](t) * t
x = [Link](...) # ou t = ...
y1 = f1(x) # ou x = ... [Link](x, y)
y2 = f2(x) # ou y = ...
[Link](y1, y2)

Matthieu Moy (Ensimag) NumPy octobre 2016 < 49 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 50 / 67 >

Animations
matplotlib/oscilateur_anim.py
import numpy as np
import [Link] as plt
import [Link] as animation
Joli, mais on peut faire bouger # figure = tout le contenu de la fenêtre
tout ça ? fig = [Link]() # A appeler avant de
# tracer autre chose
; Oui ! ! !
# Une courbe fixe
x = [Link](0, 10 * [Link], 500)
y1 = [Link](x)

parfait = [Link](x, y1)

y2 = y1 # Pas d'amortissement pour l'instant
amorti = [Link](x, y2)

Matthieu Moy (Ensimag) NumPy octobre 2016 < 51 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 52 / 67 >
 Animations Animations : à essayer
matplotlib/oscilateur_anim.py matplotlib/oscilateur_anim.py
[Link]('Oscillateurs libre et amorti')
[Link]('x')
[Link]('y = f(x)')
Ajouter une instruction print(i) dans la fonction animate
# Fonction appelée à chaque pas de l'animation Jouer avec les paramètres de FuncAnimation :
def animate(i): I interval (en miliseconde)
C = i + 1 # i vaut successivement 0, 1, 2, ... I frames (arguments passés à animate, par exemple
# On peut changer le titre [Link](...))
[Link]('Oscillateurs libre et amorti (C = ' I repeat=False (s’arrêter à la fin au lieu de répéter)
+ str(C) + ')')
[Link]
# Et changer les données
y2 = [Link](-x / float(C)) * y1 [Link]
# Attention, amorti est une liste a 1 élément.
amorti[0].set_data(x, y2)

ani = [Link](fig, animate, interval=100)

[Link]()

Matthieu Moy (Ensimag) NumPy octobre 2016 < 53 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 54 / 67 >

Animations : dessinons un peu Animations : dessinons un peu

matplotlib/pendule_graphique.py matplotlib/pendule_graphique.py
amplitude = 1 def animate(i):
dt = 0.05
t = i * dt # t = 0, 0.05, 0.10, ...
fig = [Link]() # On s'autorise quelques libertés par rapport
# 111 : découpage 1x1, sous-figure numéro 1 # à la physique ...
# (i.e. pas de découpage, mais on doit utiliser add_subplot angle = sin(t) * amplitude
# pour utiliser xlim et ylim)
x = (0, sin(angle))
ax = fig.add_subplot(111, autoscale_on=False,
xlim=(-1.5, 1.5), ylim=(-1.5, .5)) y = (0, -cos(angle))
[Link]() [Link]('Pendule parfait (t = ' +
str(round(t * 10) / float(10)) + ')')
x = [0, 0] # Attention, parfait est une liste a 1 élément.
y = [0, -1]
# y = f(x). 'o-' pour afficher des ronds sur les points. parfait[0].set_data(x, y)
parfait = [Link](x, y, 'o-')
ani = [Link](fig, animate,
[Link]('Pendule') interval=10)
[Link]('x')
[Link]('y')
[Link]()

Matthieu Moy (Ensimag) NumPy octobre 2016 < 55 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 56 / 67 >

Résolution d’équation différentielle Équations différentielles du second ordre

matplotlib/[Link]
import [Link]
import numpy as np Jusqu’ici : y 0 = f (y , t)
import [Link] as plt Comment faire y 00 = f (y 0 , y , t) ?

# y' = f(y, x) On pose Y = (y , y 0 )

def f(y, x):
return [Link](10 * y) * [Link](x)
On remarque que Y 0 = (y 0 , y 00 ) = (y 0 , f (y 0 , y , t))
On pose F (Y , t) = (Y [1], f (Y [1], Y [0], t))
x = [Link](0, 10, 100) et on a Y 0 = F (Y , t)
# Résolution numérique de l'équation (équation d’ordre 1 sur variables de dimensions 2)
y = [Link](f, 0, x)

# Tracé du résultat
[Link](x, y)
[Link]()

Matthieu Moy (Ensimag) NumPy octobre 2016 < 57 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 58 / 67 >

Application : pendule libre Animations : dessinons un peu

matplotlib/pendule_anim_simple.py

w02 = 10
y0 = [Link]([[Link] / 2.0, 0]) # angle, vitesse
θ00 + ω02 sin θ = 0
θ00 + ω02 sin θ = 0 (avec ω02 = constante) def f(y, t): On pose

Y = (y1 , y2 ) = (θ, θ0 )
# y[0] = angle, y[1] = vitesse angulaire

# On renvoie : (vitesse, accélération)
On pose Y = (y1 , y2 ) = (θ, θ0 )
return [Link]([y[1], -w02 * [Link](y[0])])
On remarque
Exprimer Y 0 en fonction de y1 et y2 (on n’a pas besoin de t) t = [Link](0, 100, 10000)
y = [Link](f, y0, t) Y 0 = (y2 , −ω02 sin y1 )
theta, thetadot = y[:, 0], y[:, 1]
Écrire la fonction F (Y , t) avec Y de dimension 2. On résout cette équation
fig = [Link]()
Résoudre le système avec [Link]. ax = fig.add_subplot(111, autoscale_on=False,
(θ en fonction de t) :

xlim=(-1.5, 1.5), ylim=(-1.5, .5))
Squelette : matplotlib/pendule_anim_simple.py [Link]() [Link]
[Link]('Pendule')
[Link]('x'); [Link]('y') On trace le segment
pendules = [Link]((0, sin(y0[0])), (0, -cos(y0[0])), 'o-')
(0, 0) → (sin(θ), − cos(θ)) à
def animate(i):

angle = theta[i] chaque instant :

x = (0, sin(angle))
[Link]

y = (0, -cos(angle))

pendules[0].set_data(x, y)

anim = [Link](fig, animate, interval=10)

Matthieu Moy (Ensimag) NumPy octobre 2016 < 59 / 67 > # Décommenter pour
Matthieu Moy enregistrer sous forme de vidéo:
(Ensimag) NumPy octobre 2016 < 60 / 67 >
#[Link]('pendule_anim_simple.mp4', writer='avconv', fps=30)
# En cas de problème avec writer, essayer writer='mencoder'
 Animations : amusons-nous un peu Diagramme de phase
matplotlib/pendule_anim.py matplotlib/diagramme_de_phase.py

# y' = f(y, t)
def f(y, t):
Même equa-diff

On a θ et θ0
# y[0] = angle, y[1] = vitesse angulaire

# On renvoie : (vitesse, accélération)

frottement = - .5 * y[1]

return [Link]([y[1], -w02 * [Link](y[0]) [Link](theta,
2 pendules
+ frottement])
thetadot)
Lois de la physique différentes pour les deux (frottement) t = [Link](0, 100, 10000)
y = [Link](f, y0, t)

Jouez avec les paramètres (position initiale, frottement, ...) theta, thetadot = y[:, 0], y[:, 1]

fig = [Link]()
ax = fig.add_subplot(111, autoscale_on=False,

xlim=(-5, 5), ylim=(-5, 5))
[Link]()

[Link]('Pendule : diagramme de phase')

# Le mode mathématique LaTeX est reconnu ($...$).
# Astuce Python : r'...' pour faire une chaîne de caractères où les \
# ne sont pas des caractères d'échappement : r'$\theta$ == '$\\theta$'
[Link](r'$\theta$')
[Link](r"$\theta'$")

# On n'utilise pas t, mais on trace theta et theta' l'un en fonction

# de l'autre.
[Link](theta, thetadot)
[Link]()

Matthieu Moy (Ensimag) NumPy octobre 2016 < 61 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 62 / 67 >

Diagramme de phase Autres outils intéressants

matplotlib/diagramme_de_phase_anim.py
SymPy : Calcul Symbolique
w02 = 10
y0 = [Link]([[Link], 0])
1 résolution d’équa-diff à [Link]
frict = 0.5
chaque itération SageMath (anciennement SAGE) : calcul symbolique et
# y' = f(y, t)
def f(y, t): numérique, accès unifié à ≈ 100 bibliothèques de calcul

# y[0] = angle, y[1] = vitesse angulaire
[Link]

# On renvoie : (vitesse, accélération)

frottement = - frict * y[1]

return [Link]([y[1], -w02 * [Link](y[0]) + frottement]) Aide au débug : pylint, pychecker, pyflakes, pep8, flake8
... [Link]
def animate(i):

global frict Interfaçage avec C ou C++

frict = i / float(200)

y = [Link](f, y0, t) [Link]

theta, thetadot = y[:, 0], y[:, 1]

[Link]('Pendule : diagramme de phase (frict = ' IDE plus évolués que Spyder : PyCharm, Eclipse (ou retour aux
+ str(frict) + ')')

curves[0].set_data(theta, thetadot) sources avec vi + terminal)
ani = [Link](fig, animate, interval=10)

[Link]()

Matthieu Moy (Ensimag) NumPy octobre 2016 < 63 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 65 / 67 >

Pour aller plus loin

Double pendule :
matplotlib/double_pendulum_animated.py,
Beaucoup d’autres exemples : Merci de votre attention !
[Link]
Une documentation très bien faite (NumPy, SciPy, Matplotlib, . . .) :
[Link]

Matthieu Moy (Ensimag) NumPy octobre 2016 < 66 / 67 > Matthieu Moy (Ensimag) NumPy octobre 2016 < 67 / 67 >