Modélisation de la propagation de la houle en présence de courants cisaillés et par bathymétrie variable
Méthode de résolution d’un système linéaire
Afin de présenter plus tard et en toute clarté le solveur utilisé dans cette étude, faisons ici un bref rappel de l’état de l’art sur la résolution de systèmes linéaires (une revue détaillée est disponible dans l’étude [153]) et sur ses trois grandes familles de méthodes de résolution. Si on pose p le nombre d’équations et q le nombre d’inconnues de notre problème, nous allons chercher à résoudre le système matrice-vecteur A−→x = −→b avec A ∈ R p∗q la matrice inversible des coe- cients de notre modèle théorique, −→b ∈ R q le vecteur colonne second membre, également imposé par notre modèle, et −→x ∈ R q notre vecteur solution. Concernant les caractéristiques éventuelles de A, on rappelle que A est symétrique si AT = A, qu’elle est positive si ∀ −→x R n , −→x TA−→x ≥ 0 et qu’elle est dénie si pour −→x R n , −→x TA−→x = 0 ⇒ −→x = 0. On dénira également la norme matricielle comme suit : kAk = max kxk=1 kAxk = max kxk6=0 kAxk kxk (B.1) Existence et unicité de la solution Plusieurs cas sont à considérer. Commençons par le cas le plus rare de matrices non carrées: s’il y a plus d’équations que d’inconnues (système surdéterminé), il n’y aura probablement pas de solution. Si au contraire il y a moins d’équations que d’inconnues (système sous-déterminé), il y aura une innité de solutions. Mais le plus fréquent est le cas des matrices carrées, où il y a autant d’équations que d’inconnues – c’est également le cas dans nos modèles théoriques. Si le système associé était homogène, l’inversibilité de la matrice déterminerait le nombre de solution : si elle est inversible, le système aurait une seule solution −→x = A−1−→b , si elle ne l’est pas, il y aurait une innité de solutions. Dans le cas d’un système non homogène (c’est le cas pour nos trois modèles, puisqu’un second membre intervient dans les conditions aux limites), là encore une matrice inversible va imposer une solution unique, tandis que dans le cas d’une matrice non inversible, pour qu’il y ait au moins une solution, il faudra vérifier la condition rang(A) = rang(Ab) (où la matrice Ab est formée par la matrice A à laquelle −→b est accolé). Avant d’entamer la résolution, il est donc important de vérier l’inversibilité du système. Dans notre étude, la matrice est carrée et appartient au corps des nombres complexes C, qui est un corps commutatif comme R, alors il sura de vérier au choix que le déterminant de A est non nul, ou que 0 n’est pas valeur propre de A , ou que Nx est le rang de A. Il est à noter que la transposée de A sera également inversible. Dans la présentation des méthodes générales de résolution, pour ne pas nuire à la clarté du propos, nous nous ramènerons au cas d’un 215 système réel. Le choix de la méthode de résolution va être conditionné par les caractéristiques de la matrice (positivité, symétrie, dénie, son rang, son conditionnement), impliquera divers points de contrôle (stabilité, nécessité d’un préconditionneur, convergence), et orira certaines possibilités (parallélisme).
Méthodes directes
Tout d’abord, les méthodes directes donnent la solution exacte en un nombre ni d’opérations, et peuvent s’appliquer indifféremment aux matrices denses ([73]) ou creuses ([52, 47]). Parmi elles, les formules de Cramer sont particulièrement simples d’utilisation: elles écrivent les solutions comme des rapports de déterminants. Mais le calcul de déterminant peut vite être trop gourmand en calcul. Ces formules sont donc réservées à de très petits systèmes. D’autres méthodes directes reposent sur la transformation de la matrice A, sachant que l’on ne change pas la solution d’un système linéaire lorsque l’on permute deux lignes ou deux colonnes, lorsqu’on fait une combinaison linéaires de lignes. On retrouve ainsi les procédés de triangularisation remontée, en général articulés autour des principes de pivot ou substitution, des variantes de la méthode d’élimination de Gauss: la factorisation LU et en particulier la factorisation de Cholesky, celle de Crout… Les solveurs directs ont l’avantage d’être robustes: ils fournissent la solution arithmétique exacte, mais le nombre d’opérations nécessaires peut être grand. Ils sont donc très e caces sur des systèmes de petites tailles, ou demandant des calculs simples. Il est d’ailleurs possible d’optimiser le stockage pour des systèmes creux (théorie des graphes). Les méthodes solve (Python, Octave, Matlab) utilisent par exemple une factorisation LU. Par contre, pour des systèmes de grande taille, le coût en mémoire, le poids des calculs et les erreurs de troncature viennent souvent déstabiliser ces méthodes. Du point de vue de l’optimisation numérique, la parallélisation est possible mais pas évidente. Il est à noter qu’il existe des méthodes hybrides, qui souhaitent combiner les avantages des méthodes directes et itératives. Elles se basent sur une décomposition de domaine (ce qui implique un maillage multi-grilles), ou des solveurs directs s’occupent des problèmes locaux tandis que des solveurs de krylov équilibrent les interfaces. On retrouve ainsi les méthodes de Schur, FETI-I, … pour n’en citer que quelques unes, mais n’étant pas l’objet de cette étude elles ne seront pas développées ici.
Méthodes itératives stationnaires
De façon générale, les méthodes itératives ([200]) nécessitent très peu de mémoire, elles sont donc particulièrement adaptées à la résolution de très gros problèmes. Pour notre étude, le choix d’une méthode de ce type s’est imposé de lui-même, nous en présenterons le mode de fonctionnement général. Parmi la multitude de méthodes itératives à disposition aujourd’hui, choisir la méthode la plus adaptée à notre problème n’est pas toujours évident. Si elles sont parfois plus diciles à paralléliser, les méthodes de Krylov convergent généralement plus rapidement, et c’est vers l’une d’elles que notre choix s’est porté: la méthode du Gradient Biconjugué, une extension de la méthode de Gradient Conjugué. Nous allons ainsi présenter cette dernière pour plus tard pouvoir restituer le Gradient Biconjugué dans un contexte clair. Le soin que nous allons apporter à la présentation des méthodes type Gradient Conjugué et à leur fonctionnement n’est pas anodin. L’implémentation du solveur Gradient Biconjugué pour système complexe a pris une grande importance dans cette étude, car il a été implémenté hors de toute bibliothèque préconçue, dans une politique de feuille blanche an de se confronter à la science et au problème dans toutes ses nuances. 216 Algorithme et préconditionnement d’une méthode itérative générale Les méthodes itératives stationnaires cherchent à construire une suite de vecteurs −→x k qui converge vers la solution que l’on notera −→x ∗ (pour la différencier mathématiquement de la solution exacte):
Résumé |