Algoritmo di soluzione dei sistemi tridiagonali

In algebra lineare numerica, l'algoritmo di soluzione dei sistemi tridiagonali, anche conosciuto come algoritmo di Thomas (da Llewellyn Thomas), è una forma più efficiente del metodo di eliminazione di Gauss che può essere usato per la soluzione di sistemi di equazioni nella forma $Ax=b$ dove la matrice $A$ è tridiagonale.

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

Il sistema di equazioni si può scrivere più sinteticamente nella seguente maniera:

$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},\,\!$ Dove $a_{1}=c_{n}=0$

Questo genere di matrici compare comunemente nella discretizzazione di equazioni di Poisson 1D e nelle interpolazioni con spline.

Proprietà

Per questi sistemi, è possibile ottenere la soluzione in $O(n)$ operazioni invece che nelle $O(n^{3})$ del metodo di eliminazione di Gauss. Un primo passaggio elimina la necessità delle $a_{i}$ e quindi una backward substitution abbreviata restituisce la soluzione.

L'algoritmo di Thomas non gode di stabilità numerica nel caso generale, ma questa può essere dimostrata per diversi casi speciali, come quando la matrice è diagonalmente dominante (per riga o per colonna) o definita positiva simmetrica. Per una trattazione più precisa della stabilità dell'algoritmo di Thomas, vedi Higham, teorema 9.12. Se, nel caso generale, è richiesta la stabilità, è preferibile l'eliminazione secondo Gauss con partial pivoting (GEPP).

Algoritmo

Il passaggio forward consiste nel modificare i coefficienti $c_{i}$ e $d_{i}$nella seguente maniera (i nuovi coefficienti sono denotati con l'apice):

$c'_{i}={\begin{cases}{\begin{array}{lcl}{\cfrac {c_{i}}{b_{i}}}&;&i=1\\{\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}}&;&i>1\\\end{array}}\end{cases}}\,$

$d'_{i}={\begin{cases}{\begin{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>1\\\end{array}}\end{cases}}\,$

Il secondo e ultimo passaggio è la sostituzione all'indietro (backward substitution):

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

$x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ \ \ i<n$

L'algoritmo in linguaggio VBA è presentato qua in una versione che sovrascrive i vettori dei coefficienti.

Sub TriDiagonal_Matrix_Algorithm(N%, A#(), B#(), C#(), D#(), X#())
 Dim i%, W#
 For i = 2 To N
  W = A(i) / B(i - 1)
  B(i) = B(i) - W * C(i - 1)
  D(i) = D(i) - W * D(i - 1)
 Next i
 X(N) = D(N) / B(N)
 For i = N - 1 To 1 Step -1
  X(i) = (D(i) - C(i) * X(i + 1)) / B(i)
 Next i
End Sub

Dimostrazione

L'algoritmo di Thomas è un caso particolare del metodo di eliminazione di Gauss.

Consideriamo il seguente sistema di equazioni

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

di incognite $x_{1},\ldots ,x_{n}$. Consideriamo la prima e la seconda equazione e modifichiamo la seconda equazione nel modo seguente:

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

Otteniamo quindi:

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

In questo modo abbiamo eliminato il primo elemento $x_{1}$ dalla seconda equazione. Consideriamo quindi la seconda equazione modificata e la terza equazione. Modifichiamo la terza equazione nel modo seguente:

$({\mbox{equazione 3}})\cdot (b_{2}b_{1}-c_{1}a_{2})-({\mbox{equazione 2}})\cdot a_{3}$

Otteniamo quindi:

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

In questo modo abbiamo eliminato il secondo elemento $x_{2}$ dalla terza equazione. La procedura viene ripetuta per ogni equazione: considerata l' $i$-esima equazione questa viene modificata nel modo seguente

$({\mbox{equazione i}})\cdot (coefficiente~di~x_{i-1}~dell'equazione~i-1)-({\mbox{equazione i-1}})\cdot (coefficiente~di~x_{i}~dell'equazione~i)$

In questo modo l'$n$-esima equazione ha come unica incognita $x_{n}$. Risolta questa equazione il valore di $x_{n}$ viene utilizzato per risolvere la $(n-1)$-esima equazione e così via fino alla prima equazione.

Chiaramente i coefficienti delle equazioni modificate si fanno via via più complicati ma, osservando la procedura, si può notare che questi possono essere definiti ricorsivamente:

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

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

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

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

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

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

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

dove abbiamo indicato con ${\tilde {a}}_{i},{\tilde {b}}_{i},{\tilde {c}}_{i},{\tilde {d}}_{i}$ i coefficienti modificati.

Per ottenere l'algoritmo di Thomas modifichiamo ulteriormente i coefficienti dividendo per ${\tilde {b}}_{i}$. In questo modo otteniamo:

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

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

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

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

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

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

Il sistema di equazioni con i coefficienti modificati diventa quindi:

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

Come già precedentemente detto, l'ultima equazione ha come unica incognita $x_{n}$. Risolvendola, possiamo utilizzare il valore di $x_{n}$e ottenere $x_{1},\ldots ,x_{n-1}$ mediante sostituzione all'indietro:

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

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

Varianti

In alcuni casi è necessario risolvere una versione leggermente diversa di sistema tridiagonale:

${\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 questo caso, la formula di Sherman-Morrison permette di evitare le operazioni ulteriori dell'algoritmo di Gauss e usare anche in questo caso l'algoritmo di Thomas. Il metodo richiede che si risolva una versione modificata e non ciclica del sistema, sia per l'input che per un vettore sparso di correzione, per poi combinare le soluzioni. Calcolare le due soluzioni nello stesso momento è particolarmente efficiente, dal momento che la parte forward dell'algoritmo può essere condivisa dalle due operazioni.

Se indichiamo con:

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

Allora il sistema da risolvere è:

$Ax=d$

In questo caso i coefficienti $a_{1}$ e $c_{n}$ sono, in generale, diversi da zero di conseguenza la loro presenza non permette di applicare l'algoritmo di Thomas direttamente.

Definiamo quindi $B\in \mathbb {R} ^{n\times n}$e $u,v\in \mathbb {R} ^{n}$ nel modo seguente:

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

Dove $\gamma \in \mathbb {R}$ è un parametro da scegliere. La matrice $A$ può essere riscritta $A=B+uv^{T}$. La soluzione del sistema viene quindi trovata con la seguente procedura: prima si risolvono due sistemi tridiagonali applicando l'algoritmo di Thomas:

$By=d\qquad \qquad Bq=u$

Successivamente la soluzione viene ricostruita utilizzando la formula di Shermann-Morisson:

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

L'implementazione in Dev-C++ è presentata qua in una versione che sovrascrive i vettori dei coefficienti.

typedef struct{
	double A[n+2]; 
	double B[n+2];
	double C[n+2];
	double D[n+2];
} COEFFICIENTS;

//Apply Thomas Alg., unknowns x[1],...,x[n]
void ThomasAlg(double x[n+1],COEFFICIENTS* coeff){
	double u[n+1]={},v[n+1]={};
	double y[n+1]={},q[n+1]={};
	double* A=coeff->A, *B=coeff->B,*C=coeff->C,*D=coeff->D;
	double Value=0;
	double w;
	int i;
	
	u[1]=gamma;
	u[n]=C[n];
	
	v[1]=1;
	v[n]=A[1]/gamma;
	
	//create matrix B
	A[1]=0;
	B[1]=B[1]-gamma;
	B[n]=B[n]-(C[n]*A[n])/gamma;
	C[n]=0;
	
	for(i=2;i<n+1;i++){
		w=A[i]/B[i-1];	
		B[i]=B[i]-w*C[i-1];
		D[i]=D[i]-w*D[i-1];
		u[i]=u[i]-w*u[i-1];
	}
	y[n]=D[n]/B[n];
	q[n]=u[n]/B[n];
	
	for(i=n-1;i>0;i--){
		y[i]=(D[i]-C[i]*y[i+1])/B[i];
		q[i]=(u[i]-C[i]*q[i+1])/B[i];
	}
	Value=(v[1]*y[1]+v[n]*y[n])/(1+v[1]*q[1]+v[n]*q[n]);
	
	for(i=1;i<n+1;i++){
		x[i]=y[i]-q[i]*Value;
	}
}

C'è anche un altro modo per risolvere il sistema quasi tridiagonale considerato. Consideriamo due sistemi lineari ausiliari di dimensioni $(n-1)\times (n-1)$:

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

Per comodità definiamo $u_{1}=0$ e $v_{1}=1$. Applicando l'algoritmo di Thomas ai due sistemi ausiliari troviamo le soluzioni $\{u_{2},u_{3}\dots ,u_{n}\}$ e $\{v_{2},v_{3}\dots ,v_{n}\}$; la soluzione $\{x_{1},x_{2}\dots ,x_{n}\}$ viene quindi rappresentata nella forma:

$x_{i}=u_{i}+x_{1}v_{i}\qquad i=1,2,\dots ,n$

Infatti, moltiplicando ogni equazione del secondo sistema per $x_{1}$, sommandola con la corrispondente equazione del primo sistema e usando la rappresentazione $x_{i}=u_{i}+x_{1}v_{i}$ vediamo subito che le equazioni del sistema originale sono soddisfatte per ogni $i=2,\ldots ,n$. Per quanto riguarda la prima equazione, consideriamo $x_{2}=u_{2}+x_{1}v_{2}$ e $x_{n}=u_{n}+x_{1}v_{n}$ e sostituiamo le seguenti espressioni nella prima equazione. Otteniamo quindi un'equazione che ha come unica incognita $x_{1}$:

$b_{1}x_{1}+c_{1}(u_{2}+x_{1}v_{2})+a_{1}(u_{n}+x_{1}v_{n})=d_{1}$

Risolvendola, otteniamo:

$x_{1}={\frac {d_{1}-a_{1}u_{n}-c_{1}u_{2}}{b_{1}+a_{1}v_{n}+c_{1}v_{2}}}$

L'implementazione in Dev-C++ è presentata qua in una versione che sovrascrive i vettori dei coefficienti.

typedef struct{
	double A[n+2]; 
	double B[n+2];
	double C[n+2];
	double D[n+2];
} COEFFICIENTS;

//Apply Thomas Alg., unknowns x[1],...,x[n]
void ThomasAlg(double x[n+1],COEFFICIENTS* coeff){
	double u[n+1]={},v[n+1]={};
	double* A=coeff->A, *B=coeff->B,*C=coeff->C,*D=coeff->D;
	double w,F[n+1]={};
	F[2]=-A[2];
	F[n]=-C[n];
	int i;
	u[1]=0;
	v[1]=1;
	for(i=3;i<n+1;i++){
		w=A[i]/B[i-1];	
		B[i]=B[i]-w*C[i-1];
		D[i]=D[i]-w*D[i-1];
		F[i]=F[i]-w*F[i-1];
	}
	u[n]=D[n]/B[n];
	v[n]=F[n]/B[n];
	for(i=n-1;i>1;i--){
		u[i]=(D[i]-C[i]*u[i+1])/B[i];
		v[i]=(F[i]-C[i]*v[i+1])/B[i];
	}
	x[1]=(D[1]-A[1]*u[n]-C[1]*u[2])/(B[1]+A[1]*v[n]+C[1]*v[2]);
	
	for(i=2;i<n+1;i++){
		x[i]=u[i]+x[1]*v[i];
	}
}

In entrambi i casi i sistemi ausiliari da risolvere sono tridiagonali di conseguenza la soluzione si ottiene ancora in $O(n)$ operazioni.

In altri casi, il sistema di equazioni potrebbe essere tridiagonale a blocchi (vedi matrice a blocchi), con sottomatrici più piccole disposte come elementi individuali del sistema matriciale descritto sopra (per esempio nelle equazioni di Poisson 2D). Sono stati sviluppati algoritmi efficienti per condurre l'eliminazione di Gauss anche per questi casi.

Il libro di testo Matematica Numerica di Quarteroni, Sacco e Saleri mostra una versione dell'algoritmo che sostituisce alcune delle divisioni con delle moltiplicazioni, in modo da migliorare l'efficienza dell'algoritmo in alcune architetture di calcolo.

Sono stati pubblicati solutori di sistemi tridiagonali paralleli per diverse architetture, incluse GPU.

Per una trattazione estensiva degli algoritmi di soluzione paralleli di sistemi tridiagonali e tridiagonali a blocchi, vedi Gallopoulos et al..