Livre_silo 30 août 2013 16:32 Page 186
¨
©
¨
©
¨
©
¨
©
C o p y r i g h t E y r o l l e s
186
Informatique pour tous
En voici une utilisation élémentaire, qui reprend la résolution de l’exemple 1 :
In [7]: numpy.linalg.solve([[2,2,-3],[-2,-1,-3],[6,4,4]],[[2],[-5],[16]])
Out[7]: array([[-14.], [ 21.], [ 4.]])
On note que le résultat renvoyé est un array, mais que numpy accepte en entrée des listes de
listes ⁸. La précision peut sembler meilleure qu’avec le programme précédent, mais c’est en
fait un leurre d’affichage :
In [8]: numpy.linalg.solve([[2,2,-3],[-2,-1,-3],[6,4,4]],[[2],[-5],[16]]) [0][0]
Out[8]: -14.000000000000023
D’une manière générale, on peut faire le pari que les résolutions effectuées par numpy seront
plus performantes que celles que l’on va coder soi-même sans finesse, tant pour la qualité
du résultat que pour le temps de calcul.
Pour autant, on évitera de s’y fier comme à une vérité indiscutable. Un crash-test classique ⁹
en calcul numérique matriciel consiste à inverser les matrices de Hilbert. On verra que ces
matrices sont très mal conditionnées et que tout calcul numérique portant sur de telles
matrices peut induire de grosses erreurs d’approximation.
La matrice de Hilbert d’ordre n ∈ N
∗ est la matrice H n ∈ M n (R) de terme général
h i,j =
1
i+j−1 pour 1 ⩽ i, j ⩽ n. On peut la définir par exemple ainsi (en pensant au
décalage d’indice) :
def hilbert(n):
return [[1./(i+j+1) for j in range(n)] for i in range(n)]
En notant C la matrice de M n,1 (R) avec des zéros partout, sauf un « 1 » en dernière
position, la résolution de H n X = C va renvoyer théoriquement la dernière colonne de
H
−1
n . Un résultat classique est le caractère entier de H
−1
n . Par exemple, la « dernière »
composante de H
−1
10 (calculée avec un logiciel de calcul formel) est 44 914 183 600 :
In [9]: Y = [[0]]*9 + [[1]]
In [10]: resolution(hilbert(10),Y)[9]
Out[10]: 44909661051.85882
In [11]: numpy.linalg.solve(hilbert(10),Y)[9][0]
Out[11]: 44908633925.999527
Pour n = 20, c’est bien pire. Le résultat attendu est 48722219250572027160000
et ceux renvoyés sont respectivement 368362187621445.56 (resolution) et
−7669635576970050.0 (solve) !
8. D’une manière générale, Python est plutôt de bonne composition avec les types !
9. [Higham], chapitre 28, explique que c’est plus subtil.
¨
©
¨
©
¨
©
¨
©
C o p y r i g h t E y r o l l e s
186
Informatique pour tous
En voici une utilisation élémentaire, qui reprend la résolution de l’exemple 1 :
In [7]: numpy.linalg.solve([[2,2,-3],[-2,-1,-3],[6,4,4]],[[2],[-5],[16]])
Out[7]: array([[-14.], [ 21.], [ 4.]])
On note que le résultat renvoyé est un array, mais que numpy accepte en entrée des listes de
listes ⁸. La précision peut sembler meilleure qu’avec le programme précédent, mais c’est en
fait un leurre d’affichage :
In [8]: numpy.linalg.solve([[2,2,-3],[-2,-1,-3],[6,4,4]],[[2],[-5],[16]]) [0][0]
Out[8]: -14.000000000000023
D’une manière générale, on peut faire le pari que les résolutions effectuées par numpy seront
plus performantes que celles que l’on va coder soi-même sans finesse, tant pour la qualité
du résultat que pour le temps de calcul.
Pour autant, on évitera de s’y fier comme à une vérité indiscutable. Un crash-test classique ⁹
en calcul numérique matriciel consiste à inverser les matrices de Hilbert. On verra que ces
matrices sont très mal conditionnées et que tout calcul numérique portant sur de telles
matrices peut induire de grosses erreurs d’approximation.
La matrice de Hilbert d’ordre n ∈ N
∗ est la matrice H n ∈ M n (R) de terme général
h i,j =
1
i+j−1 pour 1 ⩽ i, j ⩽ n. On peut la définir par exemple ainsi (en pensant au
décalage d’indice) :
def hilbert(n):
return [[1./(i+j+1) for j in range(n)] for i in range(n)]
En notant C la matrice de M n,1 (R) avec des zéros partout, sauf un « 1 » en dernière
position, la résolution de H n X = C va renvoyer théoriquement la dernière colonne de
H
−1
n . Un résultat classique est le caractère entier de H
−1
n . Par exemple, la « dernière »
composante de H
−1
10 (calculée avec un logiciel de calcul formel) est 44 914 183 600 :
In [9]: Y = [[0]]*9 + [[1]]
In [10]: resolution(hilbert(10),Y)[9]
Out[10]: 44909661051.85882
In [11]: numpy.linalg.solve(hilbert(10),Y)[9][0]
Out[11]: 44908633925.999527
Pour n = 20, c’est bien pire. Le résultat attendu est 48722219250572027160000
et ceux renvoyés sont respectivement 368362187621445.56 (resolution) et
−7669635576970050.0 (solve) !
8. D’une manière générale, Python est plutôt de bonne composition avec les types !
9. [Higham], chapitre 28, explique que c’est plus subtil.