Io ho pensato di fare la decomposizione di $A$ in $LU$ riadattando quella che si fa nell'algoritmo di Thomas per le matrici tridiagonali per questa matrice, quindi ho decomposto in questo modo:
In questo modo ottengo che (riporto a destra il costo di ogni commando):
$\gamma_2=c_2-\beta_1d_1$ ($2$ flops)
$\epsilon=d_2/\alpha_(n-2)$ ($1$ flops)
$\beta_i=b_i/\alpha_(i-1)$ con $i=2,...,n-1$ ($n-2$ flops)
$\beta_n=(b_n-\epsilonc_(n-2))/\alpha_(n-1)$ ($3$ flops)
$\alpha_(1)=a_1$
$\alpha_(i)=a_i-\beta_ic_(i-1)$ con $i=2,4,5,...,n-1$ ($2(n-2)$ flops)
$\alpha_(3)=a_3-\beta_3\gamma_2$ ($2$ flops)
In totale abbiamo $3n+2$ flops.
Quindi ci rimane risolvere il sistema lineare come $LUx=b$ ovvero $Ly=b$ e $Ux=y$.
Per quanto riguarda $Ly=b$ abbiamo:
$y_1=b_1$ ($0$ flops)
$y_i=b_i-\beta_iy_(i-1)$ con $i=2,...,n-1$ ($2(n-2)$ flops)
$y_n=b_n-\beta_ny_(n-1)-\epsilony_(n-2)$ ($4$ flops)
In totale abbiamo $2n$ flops.
Invece per $Ux=y$ abbiamo:
$x_n=y_n/\alpha_n$ ($1$ flops)
$x_i=(y_i-c_ix_(i+1))/\alpha_i$ con $i=n-1,...,3$ ($3(n-3)$ flops)
$x_2=(y_2-\gamma_2x_3)/\alpha_2$ ($3$ flops)
$x_1= (y_1-c_1x_2-d_1x_3)/\alpha_1$ ($5$ flops)
In totale abbiamo $3n$ flops.
Il costo complessivo di tutto questo algoritmo è quindi di $8n+2$ flops (molto vicino a quello dell'algoritmo di Thomas). Non so se questo sia veramente il metodo migliore per risolvere questo sistema, se qualcuno saprebbe far di meglio non esisti a scrivere (forse usando la formula di Sherman-Morrison?).
