Tridiagonale-matrix-algoritme

Het tridiagonale-matrix-algoritme, kortweg TDMA genoemd, en ook bekend als het Thomas-algoritme, is een numerieke methode om een vierkant stelsel van lineaire vergelijkingen op te lossen dat wordt beschreven door een tridiagonale matrix. Dit is een matrix waarbij de elementen buiten de diagonaal en de twee nevendiagonalen alle gelijk aan nul zijn, zoals in onderstaande matrix. De methode is niets anders dan een toepassing van de Gauss-eliminatiemethode, met als doel de elementen onder de hoofddiagonaal nul te maken, waardoor het stelsel een geschikte vorm krijgt om te worden opgelost door achterwaartse substitutie.

Een tridiagonaal stelsel van ${\displaystyle n}$ vergelijkingen en ${\displaystyle n}$ onbekenden heeft als vorm:

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

waarin ${\displaystyle a_{1}=0}$ en ${\displaystyle c_{n}=0}$. In matrixvorm wordt dit:

${\displaystyle {\begin{bmatrix}b_{1}&c_{1}&0&\ldots &\ldots &0\\a_{2}&b_{2}&c_{2}&\ddots &\ddots &\vdots \\0&a_{3}&b_{3}&c_{3}&\ddots &\vdots \\\vdots &\ddots &\ddots &\ddots &\ddots &0\\0&\ldots &0&a_{n-1}&b_{n-1}&c_{n-1}\\0&\ldots &0&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}}}$

De elementen ${\displaystyle b_{i}}$ vormen hierbij samen de hoofdiagonaal of kortweg diagonaal. De elementen ${\displaystyle a_{i}}$ en ${\displaystyle c_{i}}$ bevinden zich op de zogenaamde nevendiagonalen, meer bepaald respectievelijk op de onderdiagonaal en de bovendiagonaal.

De oplossing van dergelijke stelsels kan worden uitgevoerd in ${\displaystyle O(n)}$ bewerkingen in plaats van de normale ${\displaystyle O(n^{3})}$ bewerkingen van de Gauss-methode. In een eerste stap worden de coëfficiënten ${\displaystyle a_{i}}$ weggewerkt, waarna via achterwaartse substitutie de oplossing wordt gevonden. Dit soort stelsels komt bijvoorbeeld voor bij de numerieke oplossing van de diffusievergelijking en bij het berekenen van natuurlijke splines.

In de eerste stap worden de coëfficiënten van het stelsel omgevormd, opdat alle elementen van de coëfficiëntenmatrix op de onderdiagonaal nul worden. De nieuwe coëfficiënten worden in onderstaande vergelijkingen met een accent aangeduid:

${\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}}\,}$

en:

${\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}}}$

Na deze voorwaartse stap volgt als tweede stap een achterwaartse substitutie die de oplossing levert:

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

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

Het Engelstalige artikel over dit onderwerp bevat voorbeelden van implementaties in de programmeertalen C, Matlab en Fortran. Bij het programmeren dient men ermee rekening te houden dat de nummering van de elementen van een vector of een matrix in sommige programmeertalen loopt van 1 tot ${\displaystyle n,}$ en in andere van 0 tot ${\displaystyle n-1.}$ In het laatste geval dienen de uitdrukkingen in dit artikel aangepast te worden.

De afleiding van het tridiagonaal-matrix-algoritme is gebaseerd op Gauss-eliminatie, waarbij de elementen van de onderdiagonaal nul worden gemaakt door het uitvoeren van elementaire rijoperaties. Stel dat de onbekenden van het stelsel ${\displaystyle x_{1},\ldots ,x_{n}}$ zijn en de vergelijkingen:

${\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}}}$

Dan wordt de tweede vergelijking door middel van de eerste als volgt gewijzigd:

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

Het resultaat van deze elementaire rijoperatie is:

${\displaystyle (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}}$

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

Op deze manier is ${\displaystyle x_{1}}$ verwijderd uit de tweede vergelijking. Met behulp van de tweede vergelijking wordt dan de derde op analoge manier behandeld, zodat:

${\displaystyle (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}}$

${\displaystyle (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}}$

Dit systeem van rijoperaties wordt nu verder toegepast doorheen de matrix, waardoor alle elementen op de onderdiagonaal nul worden. Als gevolg hiervan bevat de n-de rij slechts één onbekende, namelijk ${\displaystyle x_{n}}$, die dus direct kan gevonden worden. Vervolgens worden de andere onbekenden in volgorde van onder naar boven bepaald. Bij elke stap bevat de te behandelen vergelijking immers nog slechts één onbekende. De laatste onbekende die wordt bepaald is dus ${\displaystyle x_{1}.}$

Indien de coëfficiënten van de gewijzigde vergelijkingen expliciet worden opgeschreven worden ze steeds ingewikkelder. Het is daarom beter de gewijzigde coëfficiënten, hieronder aangeduid met een tilde, recursief te definiëren.

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

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

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

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

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

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

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

Om de oplossing nog vlotter te laten verlopen kan men ook nog ${\displaystyle {\tilde {b}}_{i}}$ wegdelen (indien niet deze niet nul zijn). De aldus gewijzigde coëfficiënten, aangeduid met een asterisk, zijn:

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

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

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

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

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

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

Dit geeft het volgende stelsel, nog steeds in de oorspronkelijke onbekenden:

${\displaystyle {\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}}}$

De oplossingen zijn, startend vanaf ${\displaystyle x_{n}}$ en zo door achterwaartse substitutie naar ${\displaystyle x_{1}}$:

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

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

In bepaalde situaties, zoals in het geval van periodieke randvoorwaarden, komt een lichtjes andere vorm van een tridiagonaal stelsel voor:

${\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 dat geval kan de formule van Sherlan-Morrison gebruikt worden om bijkomende bewerkingen tijdens de elementaire rijoperaties te vermijden. In andere gevallen komt een stelsel voor dat in matrixvorm een blokmatrix geeft. Ook in dit geval bestaan er aangepaste methoden van de Gauss-eliminatiemethode.

• Conte, S.D., and de Boor, C. (1981). Elementary Numerical Analysis: an algorithmic approach. McGraw-Hill, New York, ISBN 0-07-066228-2.

• Example with VBA code

• Dit artikel of een eerdere versie ervan is een (gedeeltelijke) vertaling van het artikel Tridiagonal matrix algorithm op de Engelstalige Wikipedia, dat onder de licentie Creative Commons Naamsvermelding/Gelijk delen valt. Zie de bewerkingsgeschiedenis aldaar.