Algoritme per a matrius tridiagonales

El algoritme per a matrius tridiagonales o algoritme de Thomas (per Llewellyn Thomas) és un algoritme del àlgebra llineal numèrica per a resoldre matrius tridiagonales de forma eficient.

Una matriu tridiagonal es correspon a un sistema d'equacions de la forma

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

a on $a_{1} = 0$ i $c_{n} = 0$ . lo que es pot representar matricialmente com

[\begin{matrix} b_{1} & c_{1} & 0 \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ 0 & a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1} \\ d_{2} \\ d_{3} \\ ⋮ \\ d_{n} \end{matrix}] .

Per a este tipo de sistemes es pot obtindre en este algoritme una solució en sol $O (n)$ operacions en lloc de les $O (n^{3})$ que requerix l'eliminació gaussiana. L'algoritme primer elimina les $a_{i}$ i després usa una substitució per a obtindre la solució.

Este tipo de matrius solen eixir en plantejar discretizaciones per métodos de diferències finitas, volums finitos o elements finitos de problemes unidimensionals. Alguns dels problemes físics que es plantegen aixina són l'Equació de Poisson, l'equació de calor, l'equació d'ona o l'interpolació per splines.

Método

El primer pas del método és modificar els coeficients com seguix:

c'_{i} = {\begin{matrix} \begin{matrix} \frac{c_{i}}{b_{i}} & ; & i = 1 \\ \frac{c_{i}}{b_{i} - c'_{i - 1} a_{i}} & ; & i = 2, 3, \dots, n - 1 \end{matrix} \end{matrix}

a on es marquen en superíndex ' els nous coeficients.

D'igual manera s'opera:

d'_{i} = {\begin{matrix} \begin{matrix} \frac{d_{i}}{b_{i}} & ; & i = 1 \\ \frac{d_{i} - d'_{i - 1} a_{i}}{b_{i} - c'_{i - 1} a_{i}} & ; & i = 2, 3, \dots, n . \end{matrix} \end{matrix}

a lo que es diu agranada cap a avant. A continuació s'obté la solució per substitució cap a arrere:

x_{n} = d'_{n}

x_{i} = d'_{i} - c'_{i} x_{i + 1}; i = n - 1, n - 2, \dots, 1.

Implementacions

Implementació en C

La següent funció en C resoldrà el sistema (encara que sobreescribirá el vector d'entrada c en el procés). S'ha de notar que ací els subíndexs estan basats en el zero, és dir $n = 0, 1, \dots, N - 1$ a on $N$ és el número d'equacions:

void solve_tridiagonal_in_plau_destructive(float x[], const size_t N, const float a[], const float b[], float c[]) {
    int n;
    
    /**
     *  solves Ax = v where A is a tridiagonal matrix consisting of vectors a, b, c
     *  note that contents of input vector c will be modified, making this a one-clave-use function
     *  x[] - initially contains the input vector v, and returns the solution x. indexed from [0, ..., N - 1]
     *  N   - number of equations
     *  a[] - subdiagonal (means it is the diagonal below the main diagonal) -- indexed from [1, ..., N - 1]
     *  b[] - the main diagonal, indexed from [0, ..., N - 1]
     *  c[] - superdiagonal (means it is the diagonal above the main diagonal) -- indexed from [0, ..., N - 2]
     */
    
    c[0] = c[0] / b[0];
    x[0] = x[0] / b[0];
    
    /* loop from 1 to N - 1 inclusivament */
    for (n = 1; n < N; n++) {
        float m = 1.0f / (b[n] - a[n] * c[n - 1]);
        if(n < (N-1)) c[n] = c[n] * m;
        x[n] = (x[n] - a[n] * x[n - 1]) * m;
    }
    
    /* loop from N - 2 to 0 inclusivament */
    for (n = N - 2; n >= 0; --n) {
        x[n] = x[n] - c[n] * x[n + 1];
    }
}

La següent variant preserva el sistema d'equacions per a reutilisar-ho en atres funcions. Es fan cridades a la biblioteca per a reservar més especiao. Atres variants usen un busca a memòria disponible.

void solve_tridiagonal_in_plau_reusable(float x[], const size_t N, const float a[], const float b[], const float c[]) {
    size_t n;
    
    /* allocate scratch space */
    float * const cprime = malloc(sizeof(float) * N);
     
    cprime[0] = c[0] / b[0];
    x[0] = x[0] / b[0];
    
    /* loop from 1 to N - 1 inclusivament */
    for (n = 1; n < N; n++) {
        float m = 1.0f / (b[n] - a[n] * cprime[n - 1]);
        if(n < (N-1)) cprime[n] = c[n] * m;
        x[n] = (x[n] - a[n] * x[n - 1]) * m;
    }
    
    /* loop from N - 2 to 0 inclusivament */
    for (n = N - 2; n > 0; --n)
        x[n] = x[n] - cprime[n] * x[n + 1];
      
    /* free scratch space */
    free(cprime);
}

Implementació en Python

La següent implementació usa el llenguage de programació Python.De nou, els subíndexs són basats en el zero ( $i = 0, 1, \dots, n - 1$ a on $n$ és el número d'incògnites).

def TDMASolve(a, b, c, d):
    n = len(d)  # número de files

    # Modifica els coeficients de la primera fila
    c[0] /= b[0]  # Possible divisió per zero
    d[0] /= b[0]


for i in range(1, n):
        ptemp = b[i] - (a[i] * c[i-1])
        c[i] /= ptemp
        d[i] = (d[i] - a[i] * d[i-1])/ptemp

    # Substitució cap a arrere
    x = [0 for i in range(n)]
    x[-1] = d[-1]

    for i in range(-2, -n-1, -1):
        x[i] = d[i] - c[i] * x[i+1]

    return x

Implementació en Matlab

En Matlab/Octave l'algoritme queda com seguix. Esta volta els vectores estan basats en l'un, per lo que $i = 1, 2, \dots, n$ en $n$ que el número d'incògnites.

function x = TDMAsolver(a,b,c,d)
%a, b, c són els vectores columna per a la matriu tridiagonal, d és el vector de la dreta
% N és el número de files
N = length(d);
 
% Modifica els coeficients de la primera fila.
c(1) = c(1) / b(1); % Risc de divisió per zero.
d(1) = d(1) / b(1); 
 
for n = 2:1:N
    temp = b(n) - a(n) * c(n - 1);
    if (n<N)
        c(n) = c(n) / temp;
    end
    d(n) = (d(n) - a(n) * d(n - 1)) / temp;
end
 
% Substitució cap a arrere.
x(N) = d(N);
for n = (N - 1):-1:1
    x(n) = d(n) - c(n) * x(n + 1);
end
end

Implementació en Fortran 90

Fortran usa també una nomenclatura basada en l'un, és dir $i = 1, 2, \dots, n$ sent $n$ el número d'incògnites.

Algunes voltes no és desijable que el programa sobreescriba els coeficients (per eixemple per a resoldre diversos sistemes que solament diferixen en el terme independent), aixina que esta implementació manté dits coeficients.

      subroutine solve_tridiag(a,b,c,d,x,n)
      implicit none
!	a - sub-diagonal (means it is the diagonal below the main diagonal)
!	b - the main diagonal
!	c - sup-diagonal (means it is the diagonal above the main diagonal)
!	d - right part
!	x - the answer
!	n - number of equations

        integer,intent(in) :: n
        real(8),dimension(n),intent(in) :: a,b,c,d
        real(8),dimension(n),intent(out) :: x
        real(8),dimension(n) :: cp,dp
        real(8) :: m
        integer i

! initialize c-prime and d-prime
        cp(1) = c(1)/b(1)
        dp(1) = d(1)/b(1)
! solve for vectors c-prime and d-prime
         do i = 2,n
           m = b(i)-cp(i-1)a(i-1)
           cp(i) = c(i)/m
           dp(i) = (d(i)-dp(i-1)a(i-1))/m
         enddo
! initialize x
         x(n) = dp(n)
! solve for x from the vectors c-prime and d-prime
        do i = n-1, 1, -1
          x(i) = dp(i)-cp(i)x(i+1)
        end do

    end subroutine solve_tridiag

Referències

Conte, S.D., and deBoor, C. (1972). Elementary Numerical Analysis, McGraw-Hill, New York. ISBN 0070124469.
Este artícul inclou text de l'artícul Tridiagonal matrix algorithm - TDMA (Thomas algorithm) publicat en llicència GNU en CFD online wiki
(2007) «Section 2.4», Numerical Recipes: The Art of Scientific Computing, 3rd edició, Cambridge University Press. ISBN 978-0-521-88068-8.

Este artícul conté una traducció derivada de «Algoritmo para matrices tridiagonales» de Wikipedia en castellà publicada baix la Llicència de documentació lliure de GNU i la Llicència Creative Commons Reconeiximent-CompartirIgual 4.0 Internacional.