Thomas-Algorithmus

Der Thomas-Algorithmus (nach Llewellyn Thomas) oder auch Tridiagonalmatrix-Algorithmus (TDMA) ist eine vereinfachte Form des Gaußschen Eliminationsverfahrens, der zum schnellen Lösen von linearen Gleichungssystemen mit einer Tridiagonalmatrix benutzt wird.

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

a i x i − 1 + b i x i + c i x i + 1 = d i , {\displaystyle a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},\,\!}

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

In Matrizenform wird das System geschrieben als: [ b 1 c 1 0 a 2 b 2 c 2 a 3 b 3 ⋱ ⋱ ⋱ c n − 1 0 a n b n ] [ x 1 x 2 x 3 ⋮ x n ] = [ d 1 d 2 d 3 ⋮ d n ] . {\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 diese Matrizen entstehen üblicherweise aus der Diskretisierung der eindimensionalen Poisson-Gleichung (z. B. eindimensionale Diffusionsprobleme) und aus der natürlichen kubischen Spline-Interpolation.

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

Der erste Schritt ist es, die Koeffizienten folgendermaßen rekursiv zu modifizieren (die „neuen“ modifizierten Koeffizienten sind mit einem Strich gekennzeichnet):

c i ′ = { c 1 b 1 ; i = 1 c i b i − c i − 1 ′ a i ; i = 2 , 3 , … , n − 1 {\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}}\,}

d i ′ = { d 1 b 1 ; i = 1 d i − d i − 1 ′ a i b i − c i − 1 ′ a i ; i = 2 , 3 , … , n {\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 ist der vorwärts gerichtete Durchlauf. Die Lösung ergibt sich dann durch ein Rückwärts-Einsetzverfahren:

x n = d n ′ {\displaystyle x_{n}=d'_{n}\,}

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

In manchen Situationen, insbesondere in solchen mit periodischen Randbedingungen, kann es sein, dass eine leicht veränderte Form des tridiagonalen Systems zu lösen ist:

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 , i = 2 , … , n − 1 a n x n − 1 + b n x n + c n x 1 = d n . {\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: [ b 1 c 1 a 1 a 2 b 2 c 2 a 3 b 3 ⋱ ⋱ ⋱ c n − 1 c n a n b n ] [ x 1 x 2 x 3 ⋮ x n ] = [ d 1 d 2 d 3 ⋮ d n ] . {\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 kann die Sherman-Morrison-Woodbury-Formel benutzt werden, um den Thomas-Algorithmus zu nutzen und die zusätzlichen Operationen des Gaußschen Eliminationsverfahrens zu vermeiden.

In anderen Fällen kann das Gleichungssystem block-tridiagonal sein, d. h. die Elemente des obigen Gleichungssystems sind kleinere Blockmatrizen (z. B. bei der zweidimensionalen Poisson-Gleichung). Für diese Fälle wurden ebenfalls vereinfachte Varianten der Gauß-Elimination entwickelt.

Eine weitere Variante des Thomas-Algorithmus ist der Hu-Gallash-Algorithmus, der statt des Rückwärts-Einsetzverfahrens ein Vorwärts-Einsetzverfahren nutzt.

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

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 , … , n − 1 a n x n − 1 + b n x n = d n ; i = n {\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 ( i = 2 {\displaystyle i=2} ) wird nun wie folgt mit der ersten Gleichung modifiziert:

( Gleichung 2 ) ⋅ b 1 − ( Gleichung 1 ) ⋅ a 2 {\displaystyle ({\mbox{Gleichung 2}})\cdot b_{1}-({\mbox{Gleichung 1}})\cdot a_{2}}

Dies ergibt:

( 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 {\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 x 1 {\displaystyle x_{1}} aus der zweiten Gleichung eliminiert. Mit einer analogen Vorgehensweise ergibt sich aus der modifizierten zweiten und der dritten Gleichung:

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

Verständlicherweise werden die Koeffizienten mit jedem Schritt komplizierter, wenn sie explizit dargestellt werden. Die Koeffizienten können jedoch auch rekursiv dargestellt werden:

a ~ i = 0 b ~ 1 = b 1 b ~ i = b i b ~ i − 1 − c ~ i − 1 a i c ~ 1 = c 1 c ~ i = c i b ~ i − 1 d ~ 1 = d 1 d ~ i = d i b ~ i − 1 − d ~ i − 1 a i {\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 b ~ i {\displaystyle {\tilde {b}}_{i}} herausdividiert (wenn b ~ i ≠ 0 {\displaystyle {\tilde {b}}_{i}\neq 0} ). Die neuen Koeffizienten lauten:

a i ′ = 0 b i ′ = 1 c 1 ′ = c 1 b 1 c i ′ = c i b i − c i − 1 ′ a i d 1 ′ = d 1 b 1 d i ′ = d i − d i − 1 ′ a i b i − c i − 1 ′ a i {\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:

x i + c i ′ x i + 1 = d i ′ ;   i = 1 , … , n − 1 x n = d n ′ ;   i = n {\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 nur eine Unbekannte. Bestimmt man sie, reduziert man die Anzahl der Unbekannten in der vorletzten Gleichung auf eins, so dass dieses Rückwärts-Einsetzverfahren genutzt werden kann, um alle Unbekannten zu bestimmen.

x n = d n ′ x i = d i ′ − c i ′ x i + 1 ;   i = n − 1 , n − 2 , … , 1 {\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}}}

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