Algorithme de Thomas

En algèbre linéaire appliquée à la résolution numérique d'équation, l'algorithme de Thomas (du nom de Llewellyn Thomas), est une forme simplifiée du pivot de Gauss qui peut être utilisée pour résoudre des systèmes d'équations tridiagonales. Un système tridiagonal pour n inconnues peut s'écrire

$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},$

où $a_{1}=0$ et $c_{n}=0$.

${\begin{bmatrix}b_{1}&c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}}={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}.$

Pour de tels systèmes, la solution a une complexité algorithmique de $O(n)$ au lieu de $O(n^{3})$ requis par élimination gaussienne. Un premier balayage élimine les coefficients $a_{i}$, puis une retro substitution (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 splines cubiques.

L'algorithme de Thomas n'est pas stable en général, mais il l'est dans plusieurs cas particuliers, notamment lorsque la matrice est diagonalement dominante (soit par lignes, soit par 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 nécessaire dans le cas général, l'élimination gaussienne avec pivotement partiel (GEPP) est alors recommandée.

Méthode

Le balayage vers l'avant consiste en le calcul de nouveaux coefficients comme suit, où les nouveaux coefficients sont notés par des dérivées :

$c'_{i}={\begin{cases}{\cfrac {c_{i}}{b_{i}}},&i=1,\\{\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n-1\end{cases}}$

et

$d'_{i}={\begin{cases}{\cfrac {d_{i}}{b_{i}}},&i=1,\\{\cfrac {d_{i}-a_{i}d'_{i-1}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n.\end{cases}}$

La solution est alors obtenue par substitution à rebours :

$x_{n}=d'_{n},$

$x_{i}=d'_{i}-c'_{i}x_{i+1},\quad i=n-1,n-2,\ldots ,1.$

La méthode ci-dessus ne modifie pas les vecteurs de coefficients d'origine, mais doit stocker les nouveaux coefficients. Si les vecteurs de coefficients doivent être modifiés, alors un algorithme avec moins de gestion de l'information est possible :

Pour $i=2,3,\dots ,n,$ on détermine

$w={\cfrac {a_{i}}{b_{i-1}}},$

$b_{i}:=b_{i}-wc_{i-1},$

$d_{i}:=d_{i}-wd_{i-1},$

suivi du remplacement à rebours

$x_{n}={\cfrac {d_{n}}{b_{n}}},$

$x_{i}={\cfrac {d_{i}-c_{i}x_{i+1}}{b_{i}}}\quad {\text{pour }}i=n-1,n-2,\dots ,1.$

Calcul

L'algorithme matriciel tridiagonal est un cas particulier d'élimination gaussienne.

On suppose que les inconnues soient $x_{1},\ldots ,x_{n}$, et que les équations à résoudre sont :

${\begin{alignedat}{4}&&&b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&&&=d_{n}\,.\end{alignedat}}$

En modifiant la deuxième équation ( $i=2$ ) avec la première équation comme suit :

$({\mbox{équation 2}})\cdot b_{1}-({\mbox{équation 1}})\cdot a_{2}$

ce qui donnerait :

$(b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}=d_{2}b_{1}-d_{1}a_{2}.$

On peut remarquer que $x_{1}$ a été éliminé de la deuxième équation. Utiliser une tactique similaire avec la deuxième équation modifiée sur la troisième équation donne :

$(b_{3}(b_{2}b_{1}-c_{1}a_{2})-c_{2}b_{1}a_{3})x_{3}+c_{3}(b_{2}b_{1}-c_{1}a_{2})x_{4}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}.\,$

Cette fois $x_{2}$ a été éliminé. Si cette procédure est répétée pour l'ensemble des $n$ lignes, alors chaque équation modifiée n'impliquerait qu'une seule inconnue, $x_{n}$. Cette équation peut être résolue puis utilisée pour résoudre l'équation $(n-1)$, et ainsi de suite jusqu’à ce que toutes les inconnues soient obtenues.

On observe que les coefficients des équations modifiées deviennent de plus en plus compliqués s’ils sont exprimés de façon explicite. En examinant la procédure, les coefficients modifiés (notés par des tildes) peuvent plutôt être définis de manière récursive :

${\tilde {a}}_{i}=0\,$

${\tilde {b}}_{1}=b_{1}\,$

${\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,$

${\tilde {c}}_{1}=c_{1}\,$

${\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,$

${\tilde {d}}_{1}=d_{1}\,$

${\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,$

Pour accélérer encore la résolution, les coefficients ${\tilde {b}}_{i}$ peuvent être divisés (s'il n'y a pas de risque de division par zéro), les coefficients modifiés les plus récents, ici notés par une dérivée, seront :

$a'_{i}=0\,$

$b'_{i}=1\,$

$c'_{1}={\frac {c_{1}}{b_{1}}}\,$

$c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,$

$d'_{1}={\frac {d_{1}}{b_{1}}}\,$

$d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,$

Cela donne le système suivant :

${\begin{array}{lcl}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\x_{n}=d'_{n}\qquad &;&\ i=n.\\\end{array}}\,$

La dernière équation ne comportant qu’une seule inconnue, en la résolvant on réduit l'avant-dernière équation à une inconnue, de sorte que cette substitution à rebours puisse être utilisée pour trouver toutes les inconnues :

$x_{n}=d'_{n}\,$

$x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.$

Références

• W. H. Press, S. A. Teukolsky, W. T. Vetterling et B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, New York, Cambridge University Press, 2007 (ISBN 978-0-521-88068-8), « Section 2.4 »