Algorithme de matrice tridiagonale

En algèbre linéaire numérique , l' algorithme de matrice tridiagonale , également connu sous le nom d' algorithme de Thomas (du nom de Llewellyn Thomas ), est une forme simplifiée d' élimination gaussienne qui peut être utilisée pour résoudre des systèmes d'équations tridiagonales . Un système tridiagonal pour n inconnues peut s'écrire comme

où et .

Pour de tels systèmes, la solution peut être obtenue par des opérations au lieu de l'élimination gaussienne requise . Un premier balayage élimine les , puis une substitution arrière (abrégée) produit la solution. Des exemples de telles matrices proviennent généralement de la discrétisation de l'équation de Poisson 1D et de l'interpolation par spline cubique naturelle .

L'algorithme de Thomas n'est pas stable en général, mais il l'est dans plusieurs cas particuliers, comme lorsque la matrice est diagonalement dominante (soit par des lignes, soit par des colonnes) ou symétrique définie positive ; pour une caractérisation plus précise de la stabilité de l'algorithme de Thomas, voir le théorème de Higham 9.12. Si la stabilité est requise dans le cas général, l'élimination gaussienne avec pivotement partiel (GEPP) est recommandée à la place.

Le balayage vers l'avant consiste à calculer de nouveaux coefficients comme suit, en désignant les nouveaux coefficients par des nombres premiers :

et

La solution est alors obtenue par substitution arrière :

La méthode ci-dessus ne modifie pas les vecteurs de coefficients d'origine, mais doit également garder une trace des nouveaux coefficients. Si les vecteurs de coefficients peuvent être modifiés, alors un algorithme avec moins de comptabilité est :

Pour faire

suivi du remplacement arrière

L'implémentation sous forme de fonction C , qui utilise l'espace de travail pour éviter de modifier ses entrées pour ac, ce qui permet de les réutiliser :

void thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double scratch [ restrict X ]) { /*  résout Ax = d, où A est une matrice tridiagonale composée des vecteurs a, b, c  X = nombre d'équations  x[] = contient initialement l'entrée v, et renvoie x. indexé de [0, ..., X - 1]  a[] = sous-diagonale, indexée de [1, ..., X - 1]  b[] = diagonale principale, indexée de [0, ..., X - 1]  c[] = superdiagonale, indexée de [0, ..., X - 2]  scratch[] = espace de scratch de longueur X, fourni par l'appelant, permettant à a, b, c d'être const  non exécuté dans cet exemple : élimination manuelle coûteuse des sous-expressions communes  */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; x [ 0 ] = x [ 0 ] / b [ 0 ];
/* boucle de 1 à X - 1 inclus */ 
for ( int ix = 1 ; ix < X ; ix ++ ) { if ( ix < X -1 ){ scratch [ ix ] = c [ ix ] / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); } x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); }
/* boucle de X - 2 à 0 inclus */ 
for ( int ix = X - 2 ; ix >= 0 ; ix -- ) x [ ix ] -= scratch [ ix ] * x [ ix + 1 ]; }

La dérivation de l'algorithme de matrice tridiagonale est un cas particulier d' élimination gaussienne .

Supposons que les inconnues soient , et que les équations à résoudre soient :

Envisagez de modifier la deuxième ( ) équation avec la première équation comme suit :

ce qui donnerait :

Notez que cela a été éliminé de la deuxième équation. En utilisant une tactique similaire avec la deuxième équation modifiée sur la troisième équation, on obtient :

Cette fois, on a éliminé la règle. Si l'on répète cette procédure jusqu'à la ligne, l' équation (modifiée) ne comportera qu'une seule inconnue, . Cette dernière peut être résolue pour puis utilisée pour résoudre l' équation, et ainsi de suite jusqu'à ce que toutes les inconnues soient résolues.

Il est clair que les coefficients des équations modifiées deviennent de plus en plus compliqués s'ils sont énoncés explicitement. En examinant la procédure, les coefficients modifiés (notés avec des tildes) peuvent plutôt être définis de manière récursive :

Pour accélérer encore le processus de résolution, on peut diviser (s'il n'y a pas de division par risque zéro), les nouveaux coefficients modifiés, chacun noté avec un prime, seront :

Cela donne le système suivant avec les mêmes inconnues et coefficients définis en termes de ceux d'origine ci-dessus :

La dernière équation ne comporte qu'une seule inconnue. Sa résolution réduit à son tour l'équation suivante à une seule inconnue, de sorte que cette substitution inverse peut être utilisée pour trouver toutes les inconnues :

Dans certaines situations, notamment celles impliquant des conditions aux limites périodiques , il peut être nécessaire de résoudre une forme légèrement perturbée du système tridiagonal :

Dans ce cas, nous pouvons utiliser la formule de Sherman-Morrison pour éviter les opérations supplémentaires d'élimination gaussienne et continuer à utiliser l'algorithme de Thomas. La méthode nécessite de résoudre une version non cyclique modifiée du système à la fois pour l'entrée et un vecteur correctif clairsemé, puis de combiner les solutions. Cela peut être fait efficacement si les deux solutions sont calculées en même temps, car la partie avant de l'algorithme de matrice tridiagonale pure peut être partagée.

Si nous indiquons par :

Le système à résoudre est alors :

Dans ce cas les coefficients et sont, en général, non nuls, donc leur présence ne permet pas d'appliquer directement l'algorithme de Thomas. On peut donc considérer et comme suit : Où est un paramètre à choisir. La matrice peut être reconstruite comme . La solution est alors obtenue de la façon suivante : on résout d'abord deux systèmes d'équations tridiagonales en appliquant l'algorithme de Thomas :

Ensuite, nous reconstruisons la solution en utilisant la formule de Shermann-Morrison :

L'implémentation sous forme de fonction C , qui utilise l'espace de travail pour éviter de modifier ses entrées pour ac, ce qui permet de les réutiliser :

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double u [ restrict X ]) { /*  résout Ax = v, où A est une matrice tridiagonale cyclique composée de vecteurs a, b, c  X = nombre d'équations  x[] = contient initialement l'entrée v, et renvoie x. indexé de [0, ..., X - 1]  a[] = sous-diagonale, régulièrement indexée de [1, ..., X - 1], a[0] est le coin inférieur gauche  b[] = diagonale principale, indexée de [0, ..., X - 1]  c[] = superdiagonale, régulièrement indexée de [0, ..., X - 2], c[X - 1] est le coin supérieur droit  cmod[], u[] = vecteurs scratch chacun de longueur X  */
/* coins inférieur gauche et supérieur droit du système tridiagonal cyclique respectivement */ 
const double alpha = a [ 0 ]; const double beta = c [ X - 1 ];
/* arbitraire, mais choisi de telle sorte que la division par zéro soit évitée */ 
const double gamma = - b [ 0 ];
cmod [ 0 ] = c [ 0 ] / ( b [ 0 ] - gamma ); u [ 0 ] = gamma / ( b [ 0 ] - gamma ); x [ 0 ] /= ( b [ 0 ] - gamma );
/* boucle de 1 à X - 2 inclus */ 
for ( int ix = 1 ; ix + 1 < X ; ix ++ ) { const double m = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ]); cmod [ ix ] = c [ ix ] * m ; u [ ix ] = ( 0.0f - a [ ix ] * u [ ix - 1 ]) * m ; x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ; }
/* gérer X - 1 */ 
const double m = 1.0 / ( b [ X - 1 ] - alpha * bêta / gamma - a [ X - 1 ] * cmod [ X - 2 ]); u [ X - 1 ] = ( alpha - a [ X - 1 ] * u [ X - 2 ]) * m ; x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;
/* boucle de X - 2 à 0 inclus */ 
for ( int ix = X - 2 ; ix >= 0 ; ix -- ) { u [ ix ] -= cmod [ ix ] * u [ ix + 1 ]; x [ ix ] -= cmod [ ix ] * x [ ix + 1 ]; }
const double fait = ( x [ 0 ] + x [ X - 1 ] * alpha / gamma ) / ( 1.0 + u [ 0 ] + u [ X - 1 ] * alpha / gamma );
/* boucle de 0 à X - 1 inclus */ 
for ( int ix = 0 ; ix < X ; ix ++ ) x [ ix ] -= fact * u [ ix ]; }

Il existe également une autre façon de résoudre la forme légèrement perturbée du système tridiagonal considéré ci-dessus. Considérons deux systèmes linéaires auxiliaires de dimension :

Pour plus de commodité, nous définissons en plus et . Nous pouvons maintenant trouver les solutions et appliquer l'algorithme de Thomas aux deux systèmes tridiagonaux auxiliaires.

La solution peut alors être représentée sous la forme :

En effet, en multipliant chaque équation du second système auxiliaire par , en ajoutant avec l'équation correspondante du premier système auxiliaire et en utilisant la représentation , on voit immédiatement que les équations numérotées à du système original sont satisfaites ; il ne reste plus qu'à satisfaire l'équation numéro . Pour cela, considérons la formule pour et et substituons et dans la première équation du système original. On obtient ainsi une équation scalaire pour :

Ainsi, nous trouvons :

L'implémentation sous forme de fonction C , qui utilise l'espace de travail pour éviter de modifier ses entrées pour ac, ce qui permet de les réutiliser :

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double v [ restrict X ]) { /* résolvez d'abord un système de longueur X - 1 pour deux membres droits, en ignorant ix == 0 */ cmod [ 1 ] = c [ 1 ] / b [ 1 ]; v [ 1 ] = - a [ 1 ] / b [ 1 ]; x [ 1 ] = x [ 1 ] / b [ 1 ];
/* boucle de 2 à X - 1 inclus */ 
for ( int ix = 2 ; ix < X - 1 ; ix ++ ) { const double m = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ]); cmod [ ix ] = c [ ix ] * m ; v [ ix ] = ( 0.0f - a [ ix ] * v [ ix - 1 ]) * m ; x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ; }
/* gérer X - 1 */ 
const double m = 1.0 / ( b [ X - 1 ] - a [ X - 1 ] * cmod [ X - 2 ]); cmod [ X - 1 ] = c [ X - 1 ] * m ; v [ X - 1 ] = ( - c [ 0 ] - a [ X - 1 ] * v [ X - 2 ]) * m ; x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;
/* boucle de X - 2 à 1 inclus */ 
for ( int ix = X - 2 ; ix >= 1 ; ix -- ) { v [ ix ] -= cmod [ ix ] * v [ ix + 1 ]; x [ ix ] -= cmod [ ix ] * x [ ix + 1 ]; }
x [ 0 ] = ( x [ 0 ] - a [ 0 ] * x [ X - 1 ] - c [ 0 ] * x [ 1 ]) / ( b [ 0 ] + a [ 0 ] * v [ X - 1 ] + c [ 0 ] * v [ 1 ]);
/* boucle de 1 à X - 1 inclus */ 
for ( int ix = 1 ; ix < X ; ix ++ ) x [ ix ] += x [ 0 ] * v [ ix ]; }

Dans les deux cas, les systèmes auxiliaires à résoudre sont véritablement tridiagonaux, de sorte que la complexité de calcul globale de la résolution du système reste linéaire par rapport à la dimension du système , c'est-à -dire aux opérations arithmétiques.

Dans d'autres situations, le système d'équations peut être tridiagonal en blocs (voir matrice en blocs ), avec des sous-matrices plus petites disposées comme les éléments individuels dans le système matriciel ci-dessus (par exemple, le problème de Poisson 2D ). Des formes simplifiées d'élimination gaussienne ont été développées pour ces situations.

Le manuel Mathématiques numériques d' Alfio Quarteroni , Sacco et Saleri, répertorie une version modifiée de l'algorithme qui évite certaines divisions (en utilisant à la place des multiplications), ce qui est bénéfique sur certaines architectures informatiques.

Des solveurs tridiagonaux parallèles ont été publiés pour de nombreuses architectures vectorielles et parallèles, y compris les GPU

Pour un traitement approfondi des solveurs tridiagonaux parallèles et tridiagonaux par blocs, voir