Résolution de systèmes linéaires / Application à l'équation de la chaleur
Le but de ce projet consiste à implémenter des algorithmes de résolution de systèmes linéaires de grande taille, et à les appliquer à un problème de résolution d'équation aux dérivées partielles. Dans ce devoir, on considère uniquement des systèmes linéaires creux (ne comportant que relativement peu d'éléments non nuls), et on exploite cette propriété pour obtenir des algorithmes plus efficaces.
- la première partie s'intéresse à une méthode de résolution directe appelée factorisation "LU";
- la seconde partie propose de tester les diverses méthodes itératives de résolution;
- enfin, la dernière partie utilise les algorithmes définis dans les deux premières parties afin de résoudre une équation aux dérivées partielles, l'équation de la chaleur.
Résolution de systèmes linéaires par la factorisation LU
La factorisation LU permet de calculer à partir d'une matrice A les facteurs L et U de telle sorte que A = LU avec L matrice triangulaire inférieure avec des 1 sur la diagonale, et U matrice triangulaire supérieure. On écrira dans cette partie le code pour calculer les facteurs L et U dans le cas particulier de matrices bandes. On considère une matrice-bande dont la largeur d de bande est inconnue, par exemple lorsque d = 2 :
On stockera ainsi la partie triangulaire inférieure de A dans les coefficients b[i, j], i indiquant le numéro de la sous-diagonale et j le numéro de colonne. La partie triangulaire supérieure de A est stocké dans les coefficients e[i, j], i indiquant le numéro de la sur-diagonale et j le numéro de la colonne. La diagonale de A est stockée dans un vecteur de coefficients c[i].
L'avantage de ce stockage tient au fait qu'on sait que les facteurs L et U (obtenus sans pivotage) sont aussi des matrices triangulaires bandes de même largeur de bande. Nous allons donc stocker ces facteurs L et U sur le même modèle que A.
Par exemple pour :
L'avantage de ce stockage tient au fait qu'on sait que les facteurs L et U (obtenus sans pivotage) sont aussi des matrices triangulaires bandes de même largeur de bande. Nous allons donc stocker ces facteurs L et U sur le même modèle que A.
Par exemple pour :
Sur cette matrice, nous obtenons la factorisation suivante :
- Écrire des fonctions permettant de convertir une matrice pleine en une matrice-bande. On vous demande donc à partir d'une matrice pleine de calculer la largeur de bande d et de stocker les coefficients b[i, j], e[i, j] et c[i]. De manière réciproque, écrire une fonction de conversion d'une matrice-bande vers une matrice pleine.
- La factorisation LU se déroule de la même manière que le pivot de Gauss, on élimine successivement la première colonne, puis la deuxième jusqu'à obtenir une matrice triangulaire supérieure qui est U. On peut alors écrire l'algorithme suivant (L est ici stockée dans des coefficients l[i, j] sur le même modèle que les coefficients b[i, j] de A) :
For i from 0 to n-2 For k from 0 to d-1 l[k,i] = b[k,i]/c[i] b[k,i] = 0 For j from k-1 downto 0 If (i+k < n) Then b[j,i+k-j] = b[j,i+k-j] - l[k,i] e[k-j-1,i+k-j] End if End for If (i+k+1 < n) Then c[i+k+1] = c[i+k+1] - l[k,i] e[k,i+k+1] End if For j from 0 to d-k-2 If (i+j+k+2 < n) Then e[j,i+j+k+2] = e[j,i+j+k+2] - l[k,i] e[j+k+1,i+j+k+2] End if End for End for End for
Écrire alors la factorisation LU d'une matrice bande A à partir de cet algorithme. A la fin de l'algorithme, les coefficients l[i, j] contiennent la partie sous-diagonale de la matrice L et les coefficients c[i] et e[i, j] représentent la matrice U. Vérifier sur l'exemple simple montré au-dessus que vous obtenez les bons facteurs L et U.
3. Ecrire la résolution du système linéaire Ax = b en utilisant les facteurs L et U précédemment calculés. Cette tâche est aisée puisqu'elle consiste à résoudre deux systèmes triangulaires successifs. Typiquement, résoudre le système Lx = y lorsque d = 1 se fait de la manière suivante :
3. Ecrire la résolution du système linéaire Ax = b en utilisant les facteurs L et U précédemment calculés. Cette tâche est aisée puisqu'elle consiste à résoudre deux systèmes triangulaires successifs. Typiquement, résoudre le système Lx = y lorsque d = 1 se fait de la manière suivante :
Il est interdit de retransformer les matrices creuses en matrices pleines pour faire ce calcul. Vérifier que vous obtenez le même résultat qu'avec linalg.solve.
Résolution de systèmes linéaires par des méthodes itératives
Cette partie du problème porte sur les méthodes itératives par opposition aux méthodes précédentes, qui étaient directes. Ces méthodes permettent d'obtenir une suite convergeant vers la solution du système linéaire, sans jamais être assuré de la convergence en un nombre fini d'étapes.
Pour éviter les boucles infinies, il est demandé d'introduire des compteurs permettant d'arrêter le calcul au bout d'un nombre fixé d'itérations.
On suppose dans cette partie que la matrice considérée est une matrice creuse dont tous les coefficients sont nuls sauf sur certaines diagonales. En conséquence on ne stockera que les diagonales de A.
Par exemple la matrice :
Par exemple la matrice :
va être stockée sous la forme d'un tableau avec uniquement les diagonales non-nulles, ainsi que les indices de ces diagonales :
d'indices : ind = [-1, 0, 2]
L'indice 0 correspond à la diagonale principale, l'indice -1 à la sous-diagonale, et l'indice +1 à la sur-diagonale. Ce sont les mêmes indices utilisés par la fonction diag de Python, pour créer une matrice avec une diagonale non-nulle, par exemple :
A = diag(rand(4), 2) reviendrait à prendre ici v = rand(4); ind = array([2])
Par la suite, on écrira par la matrice A sous la forme suivante (en supposant ici que d est l'indice tel que ind(d) = 0) :
A = diag(rand(4), 2) reviendrait à prendre ici v = rand(4); ind = array([2])
Par la suite, on écrira par la matrice A sous la forme suivante (en supposant ici que d est l'indice tel que ind(d) = 0) :
1. Écrire un algorithme permettant de générer aléatoirement des matrices à diagonale strictement dominante, avec un nombre aléatoire de diagonales non-nulles. On générera directement le tableau v et le tableau d'indices ind.
2. Implémenter la méthode de Jacobi qui consiste à calculer la solution du système Ax = b à partir d'un itéré initial x[0] (souvent le vecteur nul), puis de calculer une suite (x[n]), où x[n+1] est solution du système linéaire suivant (on suppose toujours que que ind(d) = 0) :
2. Implémenter la méthode de Jacobi qui consiste à calculer la solution du système Ax = b à partir d'un itéré initial x[0] (souvent le vecteur nul), puis de calculer une suite (x[n]), où x[n+1] est solution du système linéaire suivant (on suppose toujours que que ind(d) = 0) :
... qu'on peut résoudre immédiatement :
Tester cette méthode sur une matrice de la question précédente.
3. Implémenter la méthode de Gauss-Seidel qui consiste à résoudre le système linéaire suivant :
3. Implémenter la méthode de Gauss-Seidel qui consiste à résoudre le système linéaire suivant :
... qu'on peut résoudre immédiatement :
Tester la méthode, converge-t-elle plus vite que la méthode de Jacobi ?
4. Implémenter la méthode de relaxation qui s'obtient à partir de l'algorithme de Gauss-Seidel.
À l'étape i, on a :
4. Implémenter la méthode de relaxation qui s'obtient à partir de l'algorithme de Gauss-Seidel.
À l'étape i, on a :
Le paramètre w de relaxation appartient à l'intervalle ]0; 2[.
5. Comparer les trois méthodes sur la matrice tridiagonale suivante :
5. Comparer les trois méthodes sur la matrice tridiagonale suivante :
Pour la méthode de relaxation, on prendra :
où :
N désigne ici le nombre de lignes de la matrice A. Dessiner le résidu (la norme de Ax^n - b) en fonction du numéro de l'itération en échelle logarithmique. Quelle méthode est la plus efficace sur ce cas ?
Application à l'équation de la chaleur
Dans cette partie, on s'intéresse à la résolution d'un système linéaire particulier lié à une équation aux dérivées partielles connue : l'équation de la chaleur (en mode stationnaire). Concrètement, il s'agit de modéliser l'espace (plan) par un carré [0; 1] x [0; 1], et de calculer l'évolution de la température en chacun des points de cet espace. Si l'on note T(x, y) la température en un point du carré, et f(x, y) l'apport de chaleur extérieur en ce même point, alors l'équation est la suivante :
Notre but ici consiste à déterminer, étant donnée une fonction f, une fonction T qui soit une solution approchée de cette équation (en admettant qu'une telle solution existe). La méthode de résolution appliquée est appelée méthode des différences finies. Elle consiste à représenter l'espace [0; 1] x [0; 1] par un ensemble de N x N points choisis uniformément dans le carré.
Le domaine est alors vu comme un ensemble fini de points organisés sur un carré, et les fonctions T et f sont des matrices carrées associant à un point la valeur correspondant à ce point. De plus, les dérivées partielles sont approximées par les équations linéaires suivantes :
où h est la distance entre deux points consécutifs sur la même ligne ou la même colonne, c'est-à-dire ici
h = 1/(N+1).
Conditions aux bords de Dirichlet : Pour bien définir le problème précédent, il est nécessaire de se donner des conditions aux limites sur les bords du domaine. Les conditions les plus simples sont les suivantes :
h = 1/(N+1).
Conditions aux bords de Dirichlet : Pour bien définir le problème précédent, il est nécessaire de se donner des conditions aux limites sur les bords du domaine. Les conditions les plus simples sont les suivantes :
1. Vérifier que le problème peut se mettre alors sous la forme d'un système linéaire Ax = b où la matrice A est une matrice tridiagonale par blocs de taille N² x N² et de la forme suivante :
Noter que la solution x est un vecteur de N² lignes, qu'il s'agira de retransformer en une matrice (T[i, j]) : N x N telle que T[i, j] = x[i*N+j]. L'entier N compte le nombre de points par lignes et par colonnes dans le domaine. Écrire un algorithme permettant de générer la matrice A et le vecteur b en fonction du nombre N.
2. Résoudre le système dans le cas d'un radiateur placé au centre du carré (quel est le vecteur image ?).
Tracer l'image du résultat.
3. Résoudre le système dans le cas d'un mur chaud placé au nord. du carré (quel est le vecteur image ?).
Tracer l'image du résultat.
Pour exemple, voilà le type d'image que l'on peut obtenir en résolvant certains de ces systèmes :
2. Résoudre le système dans le cas d'un radiateur placé au centre du carré (quel est le vecteur image ?).
Tracer l'image du résultat.
3. Résoudre le système dans le cas d'un mur chaud placé au nord. du carré (quel est le vecteur image ?).
Tracer l'image du résultat.
Pour exemple, voilà le type d'image que l'on peut obtenir en résolvant certains de ces systèmes :