Algoritmo de matriz tridiagonal

Reducción mejorada para matrices específicas

En álgebra lineal numérica , el algoritmo de la matriz tridiagonal , también conocido como algoritmo de Thomas (llamado así por Llewellyn Thomas ), es una forma simplificada de eliminación gaussiana que se puede utilizar para resolver sistemas de ecuaciones tridiagonales . Un sistema tridiagonal para n incógnitas se puede escribir como

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

donde y . ${\displaystyle a_{1}=0}$ ${\displaystyle c_{n}=0}$

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

Para tales sistemas, la solución se puede obtener en operaciones en lugar de la eliminación gaussiana requerida . Un primer barrido elimina los , y luego una sustitución hacia atrás (abreviada) produce la solución. Los ejemplos de tales matrices surgen comúnmente de la discretización de la ecuación de Poisson 1D y la interpolación spline cúbica natural . ${\displaystyle O(n)}$ ${\displaystyle O(n^{3})}$ ${\displaystyle a_{i}}$

El algoritmo de Thomas no es estable en general, pero sí lo es en varios casos especiales, como cuando la matriz es diagonalmente dominante (ya sea por filas o columnas) o simétrica positiva definida ; ^{[1] [2]} para una caracterización más precisa de la estabilidad del algoritmo de Thomas, véase el Teorema de Higham 9.12. ^[3] Si se requiere estabilidad en el caso general, se recomienda en su lugar la eliminación gaussiana con pivoteo parcial (GEPP). ^[2]

Método

El barrido hacia adelante consiste en el cálculo de nuevos coeficientes de la siguiente manera, denotando los nuevos coeficientes con primos:

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

y

${\displaystyle 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 solución se obtiene entonces por sustitución hacia atrás:

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

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

El método anterior no modifica los vectores de coeficientes originales, pero también debe llevar un registro de los nuevos coeficientes. Si se pueden modificar los vectores de coeficientes, entonces un algoritmo con menos contabilidad sería:

Para hacer ${\displaystyle i=2,3,\dots ,n,}$

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

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

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

Seguido de la sustitución hacia atrás

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

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

La implementación como una función C , que utiliza espacio de trabajo para evitar modificar sus entradas para ac, lo que permite reutilizarlas:

void thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double scratch [ restrict X ]) { /*  resuelve Ax = d, donde A es una matriz tridiagonal que consiste en los vectores a, b, c  X = número de ecuaciones  x[] = inicialmente contiene la entrada v, y devuelve x. indexado desde [0, ..., X - 1]  a[] = subdiagonal, indexado desde [1, ..., X - 1]  b[] = diagonal principal, indexado desde [0, ..., X - 1]  c[] = superdiagonal, indexado desde [0, ..., X - 2]  scratch[] = espacio scratch de longitud X, proporcionado por el llamador, permitiendo que a, b, c sean constantes  no realizado en este ejemplo: eliminación manual y costosa de subexpresiones comunes  */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; x [ 0 ] = x [ 0 ] / b [ 0 ];                                  /* bucle de 1 a X - 1 inclusive */ for ( int ix = 1 ; ix < X ; ix ++ ) { if ( ix < X -1 ){ scratch [ ix ] = c [ ix ] / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); } x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); }                                             /* bucle desde X - 2 hasta 0 inclusive */ for ( int ix = X - 2 ; ix >= 0 ; ix -- ) x [ ix ] -= scratch [ ix ] * x [ ix + 1 ]; }                  

Derivación

La derivación del algoritmo de la matriz tridiagonal es un caso especial de eliminación gaussiana .

Supongamos que las incógnitas son , y que las ecuaciones a resolver son: ${\displaystyle x_{1},\ldots ,x_{n}}$

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

Consideremos modificar la segunda ecuación ( ) con la primera ecuación de la siguiente manera: ${\displaystyle i=2}$

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

Lo que daría:

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

Tenga en cuenta que se ha eliminado de la segunda ecuación. Si se utiliza una táctica similar con la segunda ecuación modificada en la tercera ecuación, se obtiene: ${\displaystyle x_{1}}$

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

Esta vez se eliminó. Si se repite este procedimiento hasta la fila, la ecuación (modificada) implicará solo una incógnita, . Esto se puede resolver y luego se puede usar para resolver la ecuación, y así sucesivamente hasta que se resuelvan todas las incógnitas. ${\displaystyle x_{2}}$ ${\displaystyle n^{th}}$ ${\displaystyle n^{th}}$ ${\displaystyle x_{n}}$ ${\displaystyle (n-1)^{th}}$

Es evidente que los coeficientes de las ecuaciones modificadas se vuelven cada vez más complicados si se los enuncia explícitamente. Si se examina el procedimiento, los coeficientes modificados (anotados con tildes) pueden definirse recursivamente:

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

Para acelerar aún más el proceso de solución, se puede dividir (si no hay división por riesgo cero), los coeficientes modificados más nuevos, cada uno anotado con un primo, serán: ${\displaystyle {\tilde {b}}_{i}}$

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

Esto da el siguiente sistema con las mismas incógnitas y coeficientes definidos en términos de los originales anteriores:

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

La última ecuación implica solo una incógnita. Al resolverla, la siguiente ecuación queda reducida a una incógnita, de modo que se puede utilizar esta sustitución hacia atrás para hallar todas las incógnitas:

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

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

Variantes

En algunas situaciones, particularmente aquellas que involucran condiciones de contorno periódicas , puede ser necesario resolver una forma ligeramente perturbada del sistema tridiagonal:

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

En este caso, podemos utilizar la fórmula de Sherman-Morrison para evitar las operaciones adicionales de eliminación gaussiana y seguir utilizando el algoritmo de Thomas. El método requiere resolver una versión no cíclica modificada del sistema tanto para la entrada como para un vector correctivo disperso, y luego combinar las soluciones. Esto se puede hacer de manera eficiente si ambas soluciones se calculan a la vez, ya que se puede compartir la parte directa del algoritmo de matriz tridiagonal pura.

Si indicamos por: $${\displaystyle A={\begin{bmatrix}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{bmatrix}},x={\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}},d={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}}$$

Entonces el sistema a resolver es: $${\displaystyle Ax=d}$$

En este caso los coeficientes y son, en general, distintos de cero, por lo que su presencia no permite aplicar directamente el algoritmo de Thomas. Por tanto, podemos considerar y como sigue: Donde es un parámetro a elegir. La matriz puede reconstruirse como . La solución se obtiene entonces de la siguiente manera: ^[4] primero resolvemos dos sistemas de ecuaciones tridiagonales aplicando el algoritmo de Thomas: ${\displaystyle a_{1}}$ ${\displaystyle c_{n}}$ ${\displaystyle B\in \mathbb {R} ^{n\times n}}$ ${\displaystyle u,v\in \mathbb {R} ^{n}}$ $${\displaystyle B={\begin{bmatrix}b_{1}-\gamma &c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}-{\frac {c_{n}a_{1}}{\gamma }}\end{bmatrix}},u={\begin{bmatrix}\gamma \\0\\0\\\vdots \\c_{n}\end{bmatrix}},v={\begin{bmatrix}1\\0\\0\\\vdots \\a_{1}/\gamma \end{bmatrix}}.}$$ ${\displaystyle \gamma \in \mathbb {R} }$ ${\displaystyle A}$ ${\displaystyle A=B+uv^{\mathsf {T}}}$ $${\displaystyle By=d\qquad \qquad Bq=u}$$

Luego reconstruimos la solución utilizando la fórmula de Shermann-Morrison : ${\displaystyle x}$ $${\displaystyle {\begin{aligned}x&=A^{-1}d=(B+uv^{T})^{-1}d=B^{-1}d-{\frac {B^{-1}uv^{T}B^{-1}}{1+v^{T}B^{-1}u}}d=y-{\frac {qv^{T}y}{1+v^{T}q}}\end{aligned}}}$$

La implementación como una función C , que utiliza espacio de trabajo para evitar modificar sus entradas para ac, lo que permite reutilizarlas:

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double u [ restrict X ]) { /*  resuelve Ax = v, donde A es una matriz tridiagonal cíclica que consiste en los vectores a, b, c  X = número de ecuaciones  x[] = inicialmente contiene la entrada v, y devuelve x. indexado desde [0, ..., X - 1]  a[] = subdiagonal, indexada regularmente desde [1, ..., X - 1], a[0] es la esquina inferior izquierda  b[] = diagonal principal, indexada desde [0, ..., X - 1]  c[] = superdiagonal, indexada regularmente desde [0, ..., X - 2], c[X - 1] es la esquina superior derecha  cmod[], u[] = vectores de referencia, cada uno de longitud X  */                           /* esquinas inferior izquierda y superior derecha del sistema tridiagonal cíclico respectivamente */ const double alpha = a [ 0 ]; const double beta = c [ X - 1 ];             /* arbitrario, pero elegido de tal manera que se evita la división por cero */ const double gamma = - b [ 0 ];      cmod [ 0 ] = alfa / ( b [ 0 ] -gamma ) ; u [ 0 ] = gamma / ( b [ 0 ] -gamma ); x [ 0 ] / = ( b [ 0 ] -gamma ) ;                   /* bucle de 1 a X - 2 inclusive */ for ( int ix = 1 ; ix + 1 < X ; ix ++ ) { const double m = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ]); cmod [ ix ] = c [ ix ] * m ; u [ ix ] = ( 0.0f - a [ ix ] * u [ ix - 1 ]) * m ; x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ; }                                                      /* manejar X - 1 */ const double m = 1.0 / ( b [ X - 1 ] - alfa * beta / gamma - beta * cmod [ X - 2 ]); u [ X - 1 ] = ( alfa - a [ X - 1 ] * u [ X - 2 ]) * m ; x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;                                                      /* bucle desde X - 2 hasta 0 inclusive */ for ( int ix = X - 2 ; ix >= 0 ; ix -- ) { u [ ix ] -= cmod [ ix ] * u [ ix + 1 ]; x [ ix ] -= cmod [ ix ] * x [ ix + 1 ]; }                            constante doble hecho = ( x [ 0 ] + x [ X - 1 ] * beta / gamma ) / ( 1.0 + u [ 0 ] + u [ X - 1 ] * beta / gamma );                         /* bucle de 0 a X - 1 inclusive */ for ( int ix = 0 ; ix < X ; ix ++ ) x [ ix ] -= fact * u [ ix ]; }              

Existe también otra forma de resolver la forma ligeramente perturbada del sistema tridiagonal considerado anteriormente. ^[5] Consideremos dos sistemas lineales auxiliares de dimensión : ${\displaystyle (n-1)\times (n-1)}$ $${\displaystyle {\begin{aligned}\qquad \ \ \ \ \ b_{2}u_{2}+c_{2}u_{3}&=d_{2}\\a_{3}u_{2}+b_{3}u_{3}+c_{3}u_{4}&=d_{3}\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=d_{i}\\\dots \\a_{n}u_{n-1}+b_{n}u_{n}\qquad &=d_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1\qquad \qquad {\begin{aligned}\qquad \ \ \ \ \ b_{2}v_{2}+c_{2}v_{3}&=-a_{2}\\a_{3}v_{2}+b_{3}v_{3}+c_{3}v_{4}&=0\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=0\\\dots \\a_{n}v_{n-1}+b_{n}v_{n}\qquad &=-c_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1}$$

Para mayor comodidad, definimos adicionalmente y . Ahora podemos encontrar las soluciones y aplicar el algoritmo de Thomas al sistema tridiagonal de dos auxiliares. ${\displaystyle u_{1}=0}$ ${\displaystyle v_{1}=1}$ ${\displaystyle \{u_{2},u_{3}\dots ,u_{n}\}}$ ${\displaystyle \{v_{2},v_{3}\dots ,v_{n}\}}$

La solución puede entonces representarse en la forma: ${\displaystyle \{x_{1},x_{2}\dots ,x_{n}\}}$ $${\displaystyle x_{i}=u_{i}+x_{1}v_{i}\qquad i=1,2,\dots ,n}$$

En efecto, multiplicando cada ecuación del segundo sistema auxiliar por , sumando con la ecuación correspondiente del primer sistema auxiliar y utilizando la representación , vemos inmediatamente que se satisfacen las ecuaciones número a del sistema original; sólo queda satisfacer la ecuación número . Para ello, considere la fórmula para y y sustituya y en la primera ecuación del sistema original. Esto produce una ecuación escalar para : ${\displaystyle x_{1}}$ ${\displaystyle x_{i}=u_{i}+x_{1}v_{i}}$ ${\displaystyle 2}$ ${\displaystyle n}$ ${\displaystyle 1}$ ${\displaystyle i=2}$ ${\displaystyle i=n}$ ${\displaystyle x_{2}=u_{2}+x_{1}v_{2}}$ ${\displaystyle x_{n}=u_{n}+x_{1}v_{n}}$ ${\displaystyle x_{1}}$ $${\displaystyle b_{1}x_{1}+c_{1}(u_{2}+x_{1}v_{2})+a_{1}(u_{n}+x_{1}v_{n})=d_{1}}$$

Así, encontramos: $${\displaystyle x_{1}={\frac {d_{1}-a_{1}u_{n}-c_{1}u_{2}}{b_{1}+a_{1}v_{n}+c_{1}v_{2}}}}$$

La implementación como una función C , que utiliza espacio de trabajo para evitar modificar sus entradas para ac, lo que permite reutilizarlas:

void cyclic_thomas ( const int X , double x [ restringe X ], const double a [ restringe X ], const double b [ restringe X ], const double c [ restringe X ], double cmod [ restringe X ], double v [ restringe X ]) { /* primero resuelve un sistema de longitud X - 1 para dos lados derechos, ignorando ix == 0 */ cmod [ 1 ] = c [ 1 ] / b [ 1 ]; v [ 1 ] = - a [ 1 ] / b [ 1 ]; x [ 1 ] = x [ 1 ] / b [ 1 ];                                          /* bucle de 2 a X - 1 inclusive */ for ( int ix = 2 ; ix < X - 1 ; ix ++ ) { const double m = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ]); cmod [ ix ] = c [ ix ] * m ; v [ ix ] = ( 0.0f - a [ ix ] * v [ ix - 1 ]) * m ; x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ; }                                                      /* manejar X - 1 */ const double m = 1.0 / ( b [ X - 1 ] - a [ X - 1 ] * cmod [ X - 2 ]); cmod [ X - 1 ] = c [ X - 1 ] * m ; v [ X - 1 ] = ( - c [ 0 ] - a [ X - 1 ] * v [ X - 2 ]) * m ; x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;                                                           /* bucle desde X - 2 hasta 1 inclusive */ for ( int ix = X - 2 ; ix >= 1 ; ix -- ) { v [ ix ] -= cmod [ ix ] * v [ ix + 1 ]; x [ ix ] -= cmod [ ix ] * x [ ix + 1 ]; }                            x [ 0 ] = ( x [ 0 ] - a [ 0 ] * x [ X - 1 ] - c [ 0 ] * x [ 1 ]) / ( b [ 0 ] + a [ 0 ] * v [ X - 1 ] + c [ 0 ] * v [ 1 ]);                         /* bucle de 1 a X - 1 inclusive */ for ( int ix = 1 ; ix < X ; ix ++ ) x [ ix ] += x [ 0 ] * v [ ix ]; }              

En ambos casos los sistemas auxiliares a resolver son genuinamente tridiagonales, por lo que la complejidad computacional global de la resolución del sistema sigue siendo lineal con respecto a la dimensión del sistema , es decir, las operaciones aritméticas. ${\displaystyle Ax=d}$ ${\displaystyle n}$ ${\displaystyle O(n)}$

En otras situaciones, el sistema de ecuaciones puede ser tridiagonal en bloques (véase matriz de bloques ), con submatrices más pequeñas dispuestas como elementos individuales en el sistema matricial anterior (por ejemplo, el problema de Poisson 2D ). Se han desarrollado formas simplificadas de eliminación gaussiana para estas situaciones. ^[6]

El libro de texto Matemáticas numéricas de Alfio Quarteroni , Sacco y Saleri enumera una versión modificada del algoritmo que evita algunas de las divisiones (utilizando en su lugar multiplicaciones), lo que resulta beneficioso en algunas arquitecturas de computadoras.

Se han publicado solucionadores tridiagonales paralelos para muchas arquitecturas vectoriales y paralelas, incluidas las GPU ^{[7] [8]}

Para un tratamiento extenso de los solucionadores tridiagonales paralelos y tridiagonales en bloque, consulte ^[9].

Referencias

1. Pradip Niyogi (2006). Introducción a la dinámica de fluidos computacional . Pearson Education India. pág. 76. ISBN 978-81-7758-764-7.
2. Biswa Nath Datta (2010). Álgebra lineal numérica y aplicaciones, segunda edición . SIAM. pág. 162. ISBN 978-0-89871-765-5.
3. Nicholas J. Higham (2002). Precisión y estabilidad de algoritmos numéricos: segunda edición . SIAM. pág. 175. ISBN 978-0-89871-802-7.
4. Batista, Milán; Ibrahim Karawia, Abdel Rahman A. (2009). "El uso de la fórmula de Sherman–Morrison–Woodbury para resolver sistemas de ecuaciones lineales de bloques cíclicos tridiagonales y de bloques cíclicos pentadiagonales". Matemáticas Aplicadas y Computación . 210 (2): 558–563. doi :10.1016/j.amc.2009.01.003. ISSN  0096-3003.
5. Ryaben'kii, Victor S.; Tsynkov, Semyon V. (2 de noviembre de 2006), "Introducción", Una introducción teórica al análisis numérico , Chapman y Hall/CRC, págs. 1–19, doi :10.1201/9781420011166-1, ISBN 978-0-429-14339-7, consultado el 25 de mayo de 2022
6. Quarteroni, Alfio ; Sacco, Ricardo; Saleri, Fausto (2007). "Sección 3.8". Matemáticas Numéricas . Springer, Nueva York. ISBN 978-3-540-34658-6.
7. Chang, L.-W.; Hwu, W,-M. (2014). "Una guía para implementar solucionadores tridiagonales en GPU". En V. Kidratenko (ed.). Cálculos numéricos con GPU . Springer. ISBN 978-3-319-06548-9.{{cite conference}}: CS1 maint: multiple names: authors list (link)
8. Venetis, es decir; Kouris, A.; Sobczyk, A.; Gallopoulos, E.; Sameh, A. (2015). "Un solucionador tridiagonal directo basado en rotaciones de Givens para arquitecturas de GPU". Computación Paralela . 49 : 101-116. doi :10.1016/j.parco.2015.03.008.
9. Gallopoulos, E.; Philippe, B.; Sameh, AH (2016). "Capítulo 5". Paralelismo en cálculos matriciales . Springer. ISBN 978-94-017-7188-7.

• Conte, SD; de Boor, C. (1972). Análisis numérico elemental . McGraw-Hill, Nueva York. ISBN 0070124469.
• Este artículo incorpora texto del artículo Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) en CFD-Wiki que está bajo la licencia GFDL .

• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Sección 2.4". Recetas numéricas: el arte de la computación científica (3.ª ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8.