Algoritmus tridiagonální matice - Tridiagonal matrix algorithm

v numerická lineární algebra, algoritmus tridiagonální matice, také známý jako Thomasův algoritmus (pojmenoval podle Llewellyn Thomas ), je zjednodušená forma Gaussova eliminace které lze použít k řešení tridiagonální soustavy rovnic. Trojúhelníkový systém pro n neznámé mohou být zapsány jako

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

kde ${displaystyle a_ {1} = 0,}$ a ${displaystyle c_ {n} = 0,}$ .

{displaystyle {egin {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}} {egin {bmatrix} {x_ {1}} {x_ {2}} {x_ {3}} vdots {x_ {n}} end {bmatrix}} = {egin {bmatrix} {d_ {1}} {d_ {2}} {d_ {3}} vdots {d_ {n}} end {bmatrix}}.}

U těchto systémů lze řešení získat v ${displaystyle O (n)}$ operace místo ${displaystyle O (n ^ {3})}$ požadováno Gaussova eliminace. První zametání vylučuje ${displaystyle a_ {i}}$ a poté (zkrácená) zpětná substituce vytvoří řešení. Příklady takových matic běžně vyplývají z diskretizace 1D Poissonova rovnice a přírodní kubický spline interpolace.

Thomasův algoritmus není stabilní obecně, ale je tomu tak v několika zvláštních případech, například když je matice diagonálně dominantní (buď podle řádků nebo sloupců) nebo symetrický pozitivní určitý;^[1]^[2] pro přesnější charakterizaci stability Thomasova algoritmu viz Higham Theorem 9.12.^[3] Je-li obecně vyžadována stabilita, Gaussova eliminace s částečné otočení Místo toho se doporučuje (GEPP).^[2]

Metoda

Průběh dopředu sestává z výpočtu nových koeficientů následujícím způsobem, označujících nové koeficienty s prvočísly:

{displaystyle c '_ {i} = {egin {cases} {egin {array} {lcl} {cfrac {c_ {i}} {b_ {i}}} &; & i = 1 {cfrac {c_ {i} } {b_ {i} -a_ {i} c '_ {i-1}}} &; & i = 2,3, tečky, n-1 end {pole}} konec {případů}},}

a

{displaystyle d '_ {i} = {egin {cases} {egin {array} {lcl} {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, tečky, n. end {pole }} konec {případů}},}

Řešení se poté získá zpětnou substitucí:

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

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

Výše uvedená metoda nemění původní vektory koeficientů, ale musí také sledovat nové koeficienty. Pokud lze vektory koeficientů upravit, pak algoritmus s menším vedením účetnictví je:

Pro

{displaystyle i = 2,3, tečky, n,}

dělat

{displaystyle {egin {array} {lcl} w = {cfrac {a_ {i}} {b_ {i-1}}} b_ {i}: = b_ {i} -wc_ {i-1} d_ { i}: = d_ {i} -wd_ {i-1} end {pole}}}

následuje zpětné střídání

{displaystyle {egin {array} {lcl} x_ {n} = {cfrac {d_ {n}} {b_ {n}}} x_ {i} = {cfrac {d_ {i} -c_ {i} x_ { i + 1}} {b_ {i}}}, quad {ext {for}} i = n-1, n-2, tečky, 1end {pole}}}

Níže je uvedena implementace v podprogramu VBA bez zachování vektorů koeficientů

Sub TriDiagonal_Matrix_Algorithm(N%, A#(), B #(), C#(), D #(), X#())    Ztlumit i%, W #    Pro i = 2 Na N        Ž = A(i) / B(i - 1)        B(i) = B(i) - Ž * C(i - 1)        D(i) = D(i) - Ž * D(i - 1)    další i    X(N) = D(N) / B(N)    Pro i = N - 1 Na 1 Krok -1        X(i) = (D(i) - C(i) * X(i + 1)) / B(i)    další iKonec Sub

Derivace

Odvození algoritmu tridiagonální matice je zvláštním případem Gaussova eliminace.

Předpokládejme, že neznámí jsou ${displaystyle x_ {1}, ldots, x_ {n}}$ , a že rovnice, které mají být vyřešeny, jsou:

{displaystyle {egin {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}}}

Zvažte úpravu druhé ( ${displaystyle i = 2}$ ) rovnice s první rovnicí následovně:

{displaystyle ({mbox {rovnice 2}}) cdot b_ {1} - ({mbox {rovnice 1}}) cdot a_ {2}}

který by dal:

{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}.}

Všimněte si, že ${displaystyle x_ {1}}$ byl vyloučen z druhé rovnice. Podobnou taktiku využijeme u upraveno druhá rovnice na třetí rovnici poskytuje:

{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}.,}

Tentokrát ${displaystyle x_ {2}}$ byl vyloučen. Pokud se tento postup opakuje až do ${displaystyle n ^ {th}}$ řádek; (upravený) ${displaystyle n ^ {th}}$ rovnice bude zahrnovat pouze jednu neznámou, ${displaystyle x_ {n}}$ . To může být vyřešeno a poté použito k vyřešení ${displaystyle (n-1) ^ {th}}$ rovnice atd., dokud nebudou vyřešeny všechny neznámé.

Je zřejmé, že koeficienty na upravených rovnicích jsou stále komplikovanější, pokud jsou výslovně uvedeny. Prozkoumáním postupu lze místo toho definovat rekurzivně upravené koeficienty (označené tildes):

{displaystyle {ilde {a}} _ {i} = 0,}

{displaystyle {ilde {b}} _ {1} = b_ {1},}

{displaystyle {ilde {b}} _ {i} = b_ {i} {ilde {b}} _ {i-1} - {ilde {c}} _ {i-1} a_ {i},}

{displaystyle {ilde {c}} _ {1} = c_ {1},}

{displaystyle {ilde {c}} _ {i} = c_ {i} {ilde {b}} _ {i-1},}

{displaystyle {ilde {d}} _ {1} = d_ {1},}

{displaystyle {ilde {d}} _ {i} = d_ {i} {ilde {b}} _ {i-1} - {ilde {d}} _ {i-1} a_ {i}.,}

Chcete-li dále urychlit proces řešení, ${displaystyle {ilde {b}} _ {i}}$ lze rozdělit (pokud neexistuje dělení nulovým rizikem), novější upravené koeficienty, každý označený prvočíslem, budou:

{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}}} .,}

Získáte následující systém se stejnými neznámými a koeficienty definovanými v podmínkách původních výše:

{displaystyle {egin {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 {pole}},}

Poslední rovnice zahrnuje pouze jednu neznámou. Jeho řešení zase sníží další poslední rovnici na jednu neznámou, takže lze tuto zpětnou substituci použít k nalezení všech neznámých:

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

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

Varianty

V některých situacích, zejména těch, které se týkají periodické okrajové podmínky, bude možná potřeba vyřešit mírně narušenou formu trojdiagonálního systému:

{displaystyle {egin {alignedat} {4} & 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 i = 2, ldots, n-1 & a_ {n} x_ {n-1} && + b_ {n} x_ {n} && + c_ {n} x_ {1} && = d_ {n} ,. end {alignedat}}}

V tomto případě můžeme využít Sherman-Morrisonův vzorec vyhnout se dalším operacím Gaussovy eliminace a stále používat Thomasův algoritmus. Metoda vyžaduje řešení upravené necyklické verze systému pro vstupní i řídký korekční vektor a následnou kombinaci řešení. To lze provést efektivně, pokud jsou obě řešení počítána najednou, protože lze sdílet přední část algoritmu čisté tridiagonální matice.

V jiných situacích může být soustava rovnic blokovat tridiagonální (vidět bloková matice ), s menšími dílčími maticemi uspořádanými jako jednotlivé prvky ve výše uvedeném maticovém systému (např. 2D Poissonův problém ). Pro tyto situace byly vyvinuty zjednodušené formy Gaussovy eliminace.^[4]

Učebnice Numerická matematika Quarteroni, Sacco a Saleri, uvádí upravenou verzi algoritmu, která se vyhýbá některým divizím (místo toho používá násobení), což je výhodné pro některé počítačové architektury.

Pro mnoho vektorových a paralelních architektur, včetně GPU, byly publikovány paralelní tridiagonální řešiče^[5]^[6]

Pro rozsáhlou léčbu paralelních trojdiagonálních a blokových tridiagonálních řešičů viz ^[7]

Reference

^ Pradip Niyogi (2006). Úvod do výpočetní dynamiky tekutin. Pearson Education India. p. 76. ISBN 978-81-7758-764-7.
^ ^A ^b Biswa Nath Datta (2010). Numerická lineární algebra a aplikace, druhé vydání. SIAM. p. 162. ISBN 978-0-89871-765-5.
^ Nicholas J. Higham (2002). Přesnost a stabilita numerických algoritmů: Druhé vydání. SIAM. p. 175. ISBN 978-0-89871-802-7.
^ Quarteroni, Alfio; Sacco, Riccardo; Saleri, Fausto (2007). „Oddíl 3.8“. Numerická matematika. Springer, New York. ISBN 978-3-540-34658-6.
^ Chang, L.-W .; Hwu, W, -M. (2014). Msgstr "Průvodce implementací tridiagonálních řešičů na GPU". V V. Kidratenko (ed.). Numerické výpočty s GPU. Springer. ISBN 978-3-319-06548-9.
^ Venetis, I.E .; Kouris, A .; Sobczyk, A .; Gallopoulos, E .; Sameh, A. (2015). "Přímý třírozměrný řešič založený na rotacích Givens pro architektury GPU". Parallel Computing. 49: 101–116.
^ Gallopoulos, E .; Philippe, B .; Sameh, A.H. (2016). „Kapitola 5“. Rovnoběžnost v maticových výpočtech. Springer. ISBN 978-94-017-7188-7.

Conte, S.D. & deBoor, C. (1972). Základní numerická analýza. McGraw-Hill, New York. ISBN 0070124469.
Tento článek včlení text z článku Tridiagonal_matrix_algorithm _-_ TDMA_ (Thomas_algorithm) na CFD-Wiki to je pod GFDL licence.

Stiskněte, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). „Oddíl 2.4“. Numerické recepty: Umění vědecké práce na počítači (3. vyd.). New York: Cambridge University Press. ISBN 978-0-521-88068-8.

[Niyogi2006-1] Pradip Niyogi (2006). Úvod do výpočetní dynamiky tekutin. Pearson Education India. p. 76. ISBN 978-81-7758-764-7.

[Datta2010-2] A ^b Biswa Nath Datta (2010). Numerická lineární algebra a aplikace, druhé vydání. SIAM. p. 162. ISBN 978-0-89871-765-5.

[Higham2002-3] Nicholas J. Higham (2002). Přesnost a stabilita numerických algoritmů: Druhé vydání. SIAM. p. 175. ISBN 978-0-89871-802-7.

[Quarteroni2007-4] Quarteroni, Alfio; Sacco, Riccardo; Saleri, Fausto (2007). „Oddíl 3.8“. Numerická matematika. Springer, New York. ISBN 978-3-540-34658-6.

[ChangHwu14-5] Chang, L.-W .; Hwu, W, -M. (2014). Msgstr "Průvodce implementací tridiagonálních řešičů na GPU". V V. Kidratenko (ed.). Numerické výpočty s GPU. Springer. ISBN 978-3-319-06548-9.

[Venetis-6] Venetis, I.E .; Kouris, A .; Sobczyk, A .; Gallopoulos, E .; Sameh, A. (2015). "Přímý třírozměrný řešič založený na rotacích Givens pro architektury GPU". Parallel Computing. 49: 101–116.

[GPS.16-7] Gallopoulos, E .; Philippe, B .; Sameh, A.H. (2016). „Kapitola 5“. Rovnoběžnost v maticových výpočtech. Springer. ISBN 978-94-017-7188-7.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

Numerická lineární algebra
Klíčové koncepty	Plovoucí bod Numerická stabilita
Problémy	Systém lineárních rovnic Maticové rozklady Násobení matic (algoritmy ) Rozdělení matice Řídké problémy
Hardware	Mezipaměť CPU TLB Algoritmus zapomínající na mezipaměť SIMD Multiprocesing
Software	MATLAB Základní podprogramy lineární algebry (BLAS) LAPACK Specializované knihovny Software pro všeobecné účely