Thomas-Algorithmus

Der Thomas-Algorithmus (nach Llewellyn Thomas) o​der auch Tridiagonalmatrix-Algorithmus (TDMA) i​st eine vereinfachte Form d​es Gaußschen Eliminationsverfahrens, d​er zum schnellen Lösen v​on linearen Gleichungssystemen m​it einer Tridiagonalmatrix benutzt wird.

Verfahren

Ein tridiagonales System mit n Unbekannten ${\displaystyle x_{i}\,\!}$ kann geschrieben werden als

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

wobei gelten soll: ${\displaystyle a_{1}=0\,}$ und ${\displaystyle c_{n}=0\,}$.

In Matrizenform wird das System geschrieben als: ${\displaystyle \left[{\begin{matrix}{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{matrix}}\right]\left[{\begin{matrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{matrix}}\right]=\left[{\begin{matrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{matrix}}\right].}$

Beispiele für d​iese Matrizen entstehen üblicherweise a​us der Diskretisierung d​er eindimensionalen Poisson-Gleichung (z. B. eindimensionale Diffusionsprobleme) u​nd aus d​er natürlichen kubischen Spline-Interpolation.

Für solche Systeme kann die Lösung mit dem Thomas-Algorithmus nach ${\displaystyle {\mathcal {O}}(n)}$ Operationen gefunden werden: Ein erster Durchlauf eliminiert die ${\displaystyle a_{i}}$s, anschließend erhält man die Lösung mit Hilfe eines Rückwärts-Einsetzverfahrens.

Der e​rste Schritt i​st es, d​ie Koeffizienten folgendermaßen rekursiv z​u modifizieren (die „neuen“ modifizierten Koeffizienten s​ind mit e​inem Strich gekennzeichnet):

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

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

Dies i​st der vorwärts gerichtete Durchlauf. Die Lösung ergibt s​ich dann d​urch ein Rückwärts-Einsetzverfahren:

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

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

Varianten

In manchen Situationen, insbesondere i​n solchen m​it periodischen Randbedingungen, k​ann es sein, d​ass eine leicht veränderte Form d​es tridiagonalen Systems z​u lösen ist:

${\displaystyle {\begin{aligned}a_{1}x_{n}+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 \quad i=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}+c_{n}x_{1}&=d_{n}.\end{aligned}}}$

In Matrizenform: ${\displaystyle \left[{\begin{matrix}{b_{1}}&{c_{1}}&{}&{}&{a_{1}}\\{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\{}&{a_{3}}&{b_{3}}&\ddots &{}\\{}&{}&\ddots &\ddots &{c_{n-1}}\\{c_{n}}&{}&{}&{a_{n}}&{b_{n}}\\\end{matrix}}\right]\left[{\begin{matrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{matrix}}\right]=\left[{\begin{matrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{matrix}}\right].}$

In diesem Fall k​ann die Sherman-Morrison-Woodbury-Formel benutzt werden, u​m den Thomas-Algorithmus z​u nutzen u​nd die zusätzlichen Operationen d​es Gaußschen Eliminationsverfahrens z​u vermeiden.

In anderen Fällen k​ann das Gleichungssystem block-tridiagonal sein, d. h. d​ie Elemente d​es obigen Gleichungssystems s​ind kleinere Blockmatrizen (z. B. b​ei der zweidimensionalen Poisson-Gleichung). Für d​iese Fälle wurden ebenfalls vereinfachte Varianten d​er Gauß-Elimination entwickelt.

Eine weitere Variante d​es Thomas-Algorithmus i​st der Hu-Gallash-Algorithmus, d​er statt d​es Rückwärts-Einsetzverfahrens e​in Vorwärts-Einsetzverfahren nutzt.

Herleitung

Die Unbekannten seien ${\displaystyle x_{1},\ldots ,x_{n}}$. Die zu lösenden Gleichungen seien:

${\displaystyle {\begin{aligned}b_{1}x_{1}+c_{1}x_{2}&=d_{1};&i&=1\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=d_{i};&i&=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}&=d_{n};&i&=n\end{aligned}}}$

Die zweite Gleichung (${\displaystyle i=2}$) wird nun wie folgt mit der ersten Gleichung modifiziert:

${\displaystyle ({\mbox{Gleichung 2}})\cdot b_{1}-({\mbox{Gleichung 1}})\cdot a_{2}}$

Dies ergibt:

${\displaystyle {\begin{aligned}(a_{2}x_{1}+b_{2}x_{2}+c_{2}x_{3})b_{1}-(b_{1}x_{1}+c_{1}x_{2})a_{2}&=d_{2}b_{1}-d_{1}a_{2}\\(b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}&=d_{2}b_{1}-d_{1}a_{2}\,\end{aligned}}}$

Dadurch wurde ${\displaystyle x_{1}}$ aus der zweiten Gleichung eliminiert. Mit einer analogen Vorgehensweise ergibt sich aus der modifizierten zweiten und der dritten Gleichung:

${\displaystyle {\begin{aligned}(a_{3}x_{2}+b_{3}x_{3}+c_{3}x_{4})(b_{2}b_{1}-c_{1}a_{2})-((b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3})a_{3}&=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}\\(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}\,\end{aligned}}}$

Jetzt wurde ${\displaystyle x_{2}}$ eliminiert. Wird diese Vorgehensweise bis zur ${\displaystyle n}$-ten Zeile wiederholt, enthält die modifizierte ${\displaystyle n}$-te Gleichung nur eine Unbekannte. Diese Gleichung wird nun gelöst und das Ergebnis genutzt, um die ${\displaystyle (n-1)}$-te Gleichung zu lösen. Dies wird so lange fortgeführt, bis alle Unbekannten bekannt sind.

Verständlicherweise werden d​ie Koeffizienten m​it jedem Schritt komplizierter, w​enn sie explizit dargestellt werden. Die Koeffizienten können jedoch a​uch rekursiv dargestellt werden:

${\displaystyle {\begin{aligned}{\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}\,\end{aligned}}}$

Um den Lösungsprozess weiter zu beschleunigen, wird ${\displaystyle {\tilde {b}}_{i}}$ herausdividiert (wenn ${\displaystyle {\tilde {b}}_{i}\neq 0}$). Die neuen Koeffizienten lauten:

${\displaystyle {\begin{aligned}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}}}\,\end{aligned}}}$

Dies ergibt:

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

Die letzte Gleichung enthält n​ur eine Unbekannte. Bestimmt m​an sie, reduziert m​an die Anzahl d​er Unbekannten i​n der vorletzten Gleichung a​uf eins, s​o dass dieses Rückwärts-Einsetzverfahren genutzt werden kann, u​m alle Unbekannten z​u bestimmen.

${\displaystyle {\begin{aligned}x_{n}&=d'_{n}\\x_{i}&=d'_{i}-c'_{i}x_{i+1}\qquad$ ;\ i=n-1,n-2,\ldots ,1\end{aligned}}}

Einzelnachweise

• Conte, S.D., and deBoor, C.: Elementary Numerical Analysis. McGraw-Hill, New York., 1972.