Le TP9 - Numericable

MPSI 831, PCSI 833
Lycée Masséna
TP 9 : Pivot de Gauss
Pour bien comprendre les références de tableaux en Python. Allez sur le site Python Tutor, section visualisation 1 . Tapez les lignes suivantes dans l’éditeur.
T=[[0]]*5
T[0][0]=1
print(T)
U=[[0] for i in range(5)]
U[0][0]=1
print(U)
Puis exécutez en pas à pas.
1
Le pivot de Gauss
L’algorithme du pivot de Gauss. On rappelle que les matrices sont représentées comme des tableaux de tableaux.
Par exemple dans le système suivant






x0
1 1 1
2
AX = Y
avec
A =  1 −1 2  , X =  x1  et Y =  9 
2 −1 1
x2
7
la matrice A et le vecteur Y sont donnés comme suit :
A=[[1, 1, 1], [1, -1, 2], [2, -1, 2]]
Y=[[2], [9], [7]]
On accède au coefficient en ligne i et en colonne j (numérotées à partir de 0) de A par A[i][j].
L’algorithme du pivot nécessite l’écriture des quatre fonctions de base suivantes :
— echange_lignes(A,i,j), qui réalise l’échange de lignes Li ↔ Lj sur la matrice A.
— transvection(A,i,j,mu), qui réalise l’opération de transvection Li ← Li + µLj sur les lignes de la matrice A.
— indice_pivot(A,i), qui cherche et renvoie l’indice de la ligne j tel que le coefficient aj,i soit maximal en valeur
absolue parmi les éléments ak,i avec i ≤ k ≤ n − 1 où n est le nombre de lignes de la matrice A.
— copie_matrice(A), qui renvoie une copie de la matrice A passée en entrée.
Il faut savoir recoder ces fonctions, et vite ! Téléchargez sur le site web 2 un squelette à compléter.
Question 1. Codez tout d’abord les trois fonctions incomplètes transvection(A,i,j,mu), indice_pivot(A,i) et
copie_matrice(A) (la fonction echange_lignes(A,i,j) est déja donnée, dans une version qui fonctionne aussi avec
les tableaux Numpy).
Question 2. On rappelle que l’algorithme du pivot de Gauss consiste à exécuter des transvections et des échanges
de lignes sur la matrice A jusqu’à la rendre triangulaire supérieure, en faisant les mêmes opérations sur Y . Réécrivez
l’algorithme du pivot partiel sous forme de fonction gauss(A,Y) tel que vu en cours, si possible sans regarder le cours.
1. que vous obtiendrez par Google, ou avec l’adresse suivante : http://pythontutor.com/visualize.html#mode=edit.
2. http://perso.numericable.fr/jules.svartz/prepa/index.html
Svartz
Page 1/3
2014/2015
MPSI 831, PCSI 833
Lycée Massena
Question 3. Vérifiez que tout fonctionne en résolvant le système AX = Y suivant :
>>> A=[[1, 1, 1], [1, -1, 2], [2, -1, 1]]
>>> Y=[[2], [9], [7]]
>>> X=gauss(A,Y)
>>> print(X)
[[1.0], [-2.0], [3.0]]
Question 4. Éxécutez le code suivant :
from random import randint
n=5
A=[[0]*n for i in range(n)]
for i in range(n):
for j in range(n):
A[i][j]=randint(0,50)
Y=[[sum(A[i])] for i in range(n)]
X=gauss(A,Y)
print(X)
La solution est sensée être un vecteur colonne constitué de 1. Qu’observez-vous ?
Question 5. Écrire une fonction matrice_hilbert(n) prenant en entrée un entier naturel strictement positif, et
1
retournant la matrice Hn de taille n × n, dont le coefficient en case (i, j) est i+j+1
(lignes et colonnes sont indexées à
partir de 0).
Question 6. On considère ici n = 5. Résoudre le système H5 X


−7.7
 −6 


 −2.1 
et


 −0.5 
0
= Y , avec pour Y les deux vecteurs ci-desous :


−7.7
 −6 


 −2.1 


 −0.4 
0
On voit donc qu’une petite variation sur Y peut induire une grande variation sur X, c’est parce que la matrice H5 est
très mal conditionnée.
Module fractions et résolution sur Q. En travaillant sur les rationnels, on n’introduit pas d’erreurs d’approximation, et on résout nos systèmes de manière exacte. L’algorithme du pivot fonctionne aussi sur les rationnels, si vous
l’avez correctement écrit. Pour travailler sur Q, il faut importer le module fractions. Pour construire un rationnel,
on donne son numérateur et son dénominateur de la façon suivante :
import fractions
r=fractions.Fraction(p,q)
Question 7. Écrire une fonction matrice_hilbert_Q(n) créant aussi une matrice de Hilbert, mais dont les éléments
sont sous forme de fractions.
Résolution de systèmes linéaires avec Numpy. En transformant nos matrices en tableaux Numpy (comme ceux
vus dans le TP sur les images), on peut utiliser les fonctions Numpy pour résoudre des systèmes. La syntaxe est la
suivante (c’est l’exemple du cours) :
import numpy as np
A=[[1, 1, 1], [1, -1, 2], [2, -1, 1]]
Y=[[2], [9], [7]]
X=np.linalg.solve(A,Y)
Svartz
Page 2/3
2014/2015
MPSI 831, PCSI 833
Lycée Massena
En fait, la conversion en tableaux Numpy se fait automatiquement à l’utilisation de np.linalg.solve. On va
maintenant travailler avec la famille des matrices Hn qui sont connues pour être très mal conditionnées c’est à dire
que les résultats fournis par l’algorithme du pivot sont très vite très éloignés de la réalité, et Numpy ne fera pas mieux
que nous.
Question 8. Résoudre, pour tout n entre 2 et 20, le système Hn X = Y , où Hn est la matrice de Hilbert d’ordre n et
Y le vecteur ayant un 1 dans sa dernière composante, et des 0 ailleurs. On affichera à l’écran la dernière composante du
vecteur X. On fera trois résolutions : sur les rationnels avec notre algorithme du pivot (pour l’affichage, on convertira
cette dernière composante en entier, c’en est un), sur les flottants avec notre algorithme, et sur les flottants avec
Numpy. Comparez.
Question 9. Comparez le temps d’exécution pour des matrices de taille ' 100 du pivot avec des rationnels et du
pivot sur flottants sur la matrice Hn . Le temps d’exécution de l’algorithme avec des rationnels n’est plus en O(n3 ),
pourquoi ?
2
Une application physique : système de masses sur un fil tendu
On considère un fil horizontal, attachés aux extrémités, sur lequel sont régulièrement réparties N masses, numérotées
de 0 à N − 1. On applique sur chaque masse une force perpendiculaire au fil, également dans le plan horizontal, comme
le montre le schéma suivant (N = 3 ici).
F~2
F~0
~ey
η2
η0
~ex
η1
F~1
On cherche à étudier les petits déplacements ηi des masses (largement exagérés sur la figure ci-dessus) dans
la direction perpendiculaire à l’axe du fil. Pour de petits déplacements, la longueur du fil entre deux masses reste
constante égale à `, et la tension du fil est également constante. Au premier ordre, l’équation du mouvement de la
masse numéro i dans la direction ~ey s’écrit :
..
mi ηi = fi −
T
T
(ηi − ηi−1 ) + (ηi+1 − ηi )
`
`
où ` est la longueur du fil entre deux masses (ou une masse et l’extrémité du fil), T est la tension du fil et F~i = fi~ey ,
avec la convention η−1 = ηN = 0.
Question 10. En convenant que T` = 1 N.m−1 pour simplifier, donner le système linéaire vérifié par (η0 , . . . , ηN −1 )
lorsque le système est à l’équilibre (les forces F~i sont de plus supposées constantes). Ce genre de système tridiagonal 3
intervient fréquemment en mécanique.
Question 11. Peut-on améliorer l’algorithme du pivot pour résoudre plus rapidement un tel système tridiagonal, c’est
à dire avec une complexité plus petite que O(n3 ) ? (On ne demande pas de le coder).
Question 12. Reprendre l’étude de la section 1 avec les matrices obtenues pour ce système, et le vecteur Y dont
les composantes sont toutes égales à 1 (pour l’affichage, on passera ce coup-ci le rationnel en flottant). Vérifiez que
l’algorithme du pivot se passe plutôt bien pour ces matrices là.
3
Reprise du TP8
Reprenez les manipulations de fichiers du TP8. Un lien vers le TP et la liste d’élèves est disponible sur la page
web.
3. C’est à dire que les seuls coefficients non nuls sont sur la diagonale principale et les deux diagonales situées en dessous et au dessus.
Svartz
Page 3/3
2014/2015

Download Report