Algoritmo de matriz tridiagonal

En álgebra lineal numérica , el algoritmo de matriz tridiagonal , también conocido como algoritmo de Thomas (llamado así en honor a 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},}$

dónde${\displaystyle a_{1}=0}$y${\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${\displaystyle O(n)}$operaciones en lugar de${\displaystyle O(n^{3})}$requerido por eliminación gaussiana . Un primer barrido elimina el${\displaystyle a_{i}}$y luego una sustitución hacia atrás (abreviada) produce la solución. Ejemplos de tales matrices surgen comúnmente de la discretización de la ecuación de Poisson 1D y la interpolación de splines cúbicos naturales .

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 definida positiva ; ^{[ 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 ]}

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

${\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 mediante 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 los vectores de coeficientes pueden modificarse, entonces un algoritmo con menos gestión de datos es:

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

${\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 temporal para evitar modificar sus entradas para ac, permitiendo que se reutilicen:

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 consta de vectores a, b, c  X = número de ecuaciones  x[] = inicialmente contiene la entrada, d, 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 de trabajo de longitud X, proporcionado por quien llama, que permite que a, b, c sean const  no realizado en este ejemplo: eliminación manual costosa de subexpresiones comunes  */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; x [ 0 ] = x [ 0 ] / b [ 0 ];/* bucle desde 1 hasta 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 ]; }

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

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

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

Considere modificar el segundo (${\displaystyle i=2}$) ecuación con la primera ecuación como sigue:

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

lo que daría como resultado:

${\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${\displaystyle x_{1}}$Se ha eliminado de la segunda ecuación. Utilizando una táctica similar con la segunda ecuación modificada en la tercera ecuación se obtiene:

${\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${\displaystyle x_{2}}$fue eliminado. Si este procedimiento se repite hasta que${\displaystyle n^{th}}$fila; la (modificada)${\displaystyle n^{th}}$La ecuación involucrará solo una incógnita,${\displaystyle x_{n}}$. Esto se puede resolver y luego utilizar para resolver el${\displaystyle (n-1)^{th}}$ecuación, y así sucesivamente hasta que se resuelvan todas las incógnitas.

Evidentemente, los coeficientes de las ecuaciones modificadas se vuelven cada vez más complicados si se expresan explícitamente. Al examinar el procedimiento, los coeficientes modificados (indicados 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,${\displaystyle {\tilde {b}}_{i}}$Se puede dividir (si no hay riesgo de división por cero ), los nuevos coeficientes modificados, cada uno denotado con una prima, serán:

${\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 como resultado 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 contiene una sola incógnita. Resolverla reduce la penúltima ecuación a una sola incógnita, de modo que esta sustitución hacia atrás puede utilizarse 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.}

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 empleando 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 realizar de manera eficiente si ambas soluciones se calculan simultáneamente, ya que la parte directa del algoritmo de matriz tridiagonal pura se puede compartir.

Si lo 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${\displaystyle a_{1}}$y${\displaystyle c_{n}}$son, en general hablando, distintos de cero, por lo que su presencia no permite aplicar directamente el algoritmo de Thomas. Por lo tanto, podemos considerar${\displaystyle B\in \mathbb {R} ^{n\times n}}$y${\displaystyle u,v\in \mathbb {R} ^{n}}$como sigue: $${\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}}.}$$ Dónde${\displaystyle \gamma \in \mathbb {R} }$es un parámetro a elegir. La matriz A se puede reconstruir como${\displaystyle A=B+uv^{\mathsf {T}}}$. La solución se obtiene entonces de la siguiente manera: ^{[ 4 ]} primero resolvemos dos sistemas de ecuaciones tridiagonales aplicando el algoritmo de Thomas: $${\displaystyle By=d\qquad \qquad Bq=u}$$

Luego reconstruimos la solución x utilizando la fórmula de Shermann-Morrison : $${\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 temporal para evitar modificar sus entradas para ac, permitiendo que se reutilicen:

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 consta de 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 prueba, 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 evite la división por cero */const double gamma = - b [ 0 ];cmod [ 0 ] = c [ 0 ] / ( b [ 0 ] - gamma );u [ 0 ] = gamma / ( b [ 0 ] - gamma );x [ 0 ] /= ( b [ 0 ] - gamma );/* bucle desde 1 hasta X - 2 inclusive */para ( 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 ;}/* maneja X - 1 */const double m = 1.0 / ( b [ X - 1 ] - alpha * beta / gamma - a [ X - 1 ] * 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 */para ( int ix = X - 2 ; ix >= 0 ; ix -- ) {u [ ix ] -= cmod [ ix ] * u [ ix + 1 ];x [ ix ] -= cmod [ ix ] * x [ ix + 1 ];}const double fact = ( x [ 0 ] + x [ X - 1 ] * alpha / gamma ) / ( 1.0 + u [ 0 ] + u [ X - 1 ] * alpha / gamma );/* bucle desde 0 hasta X - 1 inclusive */para ( int ix = 0 ; ix < X ; ix ++ )x [ ix ] -= hecho * u [ ix ];}

También existe 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}v_{i-1}+b_{i}v_{i}+c_{i}v_{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${\displaystyle u_{1}=0}$y${\displaystyle v_{1}=1}$Ahora podemos encontrar las soluciones.${\displaystyle \{u_{2},u_{3}\dots ,u_{n}\}}$y${\displaystyle \{v_{2},v_{3}\dots ,v_{n}\}}$Aplicando el algoritmo de Thomas al sistema tridiagonal auxiliar de dos vías.

La solución${\displaystyle \{x_{1},x_{2}\dots ,x_{n}\}}$entonces se puede representar de la forma: $${\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${\displaystyle x_{1}}$, sumando con la ecuación correspondiente del primer sistema auxiliar y utilizando la representación${\displaystyle x_{i}=u_{i}+x_{1}v_{i}}$, inmediatamente vemos que las ecuaciones númeroSe satisfacen las ecuaciones 2 a n del sistema original; solo queda satisfacer la ecuación número1. Para ello, considere la fórmula para${\displaystyle i=2}$y${\displaystyle i=n}$y sustituir${\displaystyle x_{2}=u_{2}+x_{1}v_{2}}$y${\displaystyle x_{n}=u_{n}+x_{1}v_{n}}$en la primera ecuación del sistema original. Esto produce una ecuación escalar para${\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}}$$

Por lo tanto, 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 temporal para evitar modificar sus entradas para ac, permitiendo que se reutilicen:

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 v [ restrict X ]) {/* Primero, resuelve un sistema de longitud X - 1 para dos segundos miembros, ignorando ix == 0 */cmod [ 1 ] = c [ 1 ] / b [ 1 ];v [ 1 ] = - a [ 1 ] / b [ 1 ];x [ 1 ] = x [ 1 ] / b [ 1 ];/* bucle desde 2 hasta X - 1 inclusive */para ( 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 ;}/* maneja 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 */para ( 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 desde 1 hasta X - 1 inclusive */para ( int ix = 1 ; ix < X ; ix ++ )x [ ix ] += x [ 0 ] * v [ ix ];}

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

En otras situaciones, el sistema de ecuaciones puede ser tridiagonal por bloques (véase matriz por 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, incluye una versión modificada del algoritmo que evita algunas de las divisiones (utilizando multiplicaciones en su lugar), lo cual resulta beneficioso en algunas arquitecturas informáticas.

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

Para un tratamiento extenso de solucionadores tridiagonales paralelos y tridiagonales por bloques, consulte ^{[ 9 ].}

• 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» . Numerical Recipes: The Art of Scientific Computing (3.ª  ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8.