Programmation vectorielle

Le but des exercices est

• d'avoir un programme qui donne la bonne réponse
• qui soit le plus rapide possible (et pour cela on utilise massivement Numpy)

En règle général si vous avez des for imbriqués c'est mauvais signe.

import numpy as np

np.set_printoptions(precision=10, linewidth=150, suppress=True)

Méthode du pivot de Gauss partiel

L'ennoncé est dans le cours. On verifiera sur le cas qui génère des erreurs d'arrondis.

Factorisation de Choleski

Il s'agit de décomposer A en $A = B\, B^T$ avec B une matrice triangulaire inférieure. Cela n'est possible que si la matrice A est symétrique et définie positive (c'est d'ailleurs un facon de vérifier qu'une matrice est définie positive).

Écrire l'algorithme de Choleski qui prend A et retourne B (pour deviner l'algorithme, essayez de trouver les coefficients de B à partir des coefficients de A sur une matrice A 4x4).

Rappel : pas de boucles for imbriquées mais des opérations vectorielles et matricielles (sur des sous-matrices).

Créer une matrice symétrique définie positive et vérifier que sonprogramme marche bien.

Amméliorer Jacobi

Lorsqu'on écrit une itération de la méthode de Jacobi avec l'ensemble des coefficients, on constate que si on calcule la nouvelle valeur de x élément par élement alors lorsqu'on veut mettre à jour x[1], on connait déjà x[0]. Idem lorsqu'on met à jour x[2] on connait déjà x[0] et x[1], etc.

L'idée de la version amméliorée de Jacobi est d'utiliser la nouvelle valeur de x[0] et non pas l'ancienne comme c'est le cas dans l'algorithme du cours. Ainsi en utilisant des valeurs mise à jour on peut espérer converger plus vite.

Dans ce cas, chaque itération demande un calcul ligne par ligne et donc une boucle for dans la boucle for sur les itérations.

Test d'arrêt

On ajoutera un argument error à la fonction qui indique la précision désirée du résultat. Malheureusement pour connaitre la précision de notre calcul il faut connaitre la solution. On pourrait aussi calculer ${\bf A \, x}^t - {\bf b}$ mais ce n'est pas non plus l'écart entre ${\bf x}^t$ et la solution (en plus c'est un calcul en n² donc lourd).

Aussi on va regarder quand la convergence ralentit et donc on s'arrête lorsque $||{\bf x}^{t+1} - {\bf x}^t|| < \texttt{error}$.

Bonus

Réécrire la méthode sans la boucle for mais en prenant bien le nouvelles valeurs de x. Pour cela on va découper la matrice A en deux matrices triangulaires.

# les méthodes de Numpy triu et tril seront utiles