Una matrice tridiagonale

C'è una particolare matrice tridiagonale che viene largamente studiata perché strettamente legata a certi importanti problemi fisici. Si tratta di una matrice che per ogni riga i ha gli elementi a_ii=2 e a_i,i+1=a_i,i-1=-1; i rimanenti elementi sono nulli. Mostriamo in quali applicazioni si ottiene.

Consideriamo un problema che si presenta spesso in applicazioni fisiche, quello della risoluzione (numerica) di problemi differenziali con condizioni ai limiti del tipo

$$y'' = -f(x) \hspace{2cm}(0\leq x\leq1), \hspace{1cm}y(0)=y(1)=0.$$

con f(x) funzione assegnata. Uno dei modi per trovare la soluzione di questa equazione consiste nel considerare l'equazione non sull'intervallo (continuo) [0,1]ma nell'insieme (discreto) di punti x_i = ih, per $i=0,\ldots,n+1$ e h=1/(n+1), e di approssimare l'operazione di ``derivata seconda'' utilizzando solo questi valori x_i. Questo metodo si chiama delle differenze finite.

Per approssimare la derivata seconda nel discreto si parte dalla formula di Taylor calcolata in x_i+h:

y(x_i+h) = y(x_i) + h y'(x_i) + h² y''(x_i)/2 + h³ y'''(x_i)/6 + o(h³).

Calcolando la stessa formula in x_i-h e sommando le due espressioni ottenute, i termini y' e y''' si cancellano e ricaviamo

y(x_i+h)+y(x_i-h) = 2y(x_i) + h² y''(x_i) + o(h³).

Adesso trascuriamo il termine o(h³) e poniamo y_i=y(x_i). Otteniamo così che

$$y''(x_i) \approx { y_{i-1}-2 y_i+y_{i+1} \over h^2} .$$

L'equazione che vogliamo risolvere si trasforma così nel sistema di equazioni

$${ y_{i-1}-2y_i+y_{i+1} \over h^2} = -f(x_i) \hspace{1cm}(i=0,\ldots,n),\hspace{2cm} y_0=y_{n+1}=0,$$

dove le incognite sono $y_1, \ldots,y_n$. Abbiamo cioè ottenuto un sistema lineare con una matrice tridiagonale di ordine n che chiameremo A_n:

$$\left[ \matrix{ 2 & -1 & 0 & \cdots & 0 \cr -1 & 2 & -1 & & \... ...{ f(x_1)\cr f(x_2)\cr f(x_3)\cr \vdots \cr f(x_n) } \right] .$$

Questa matrice (o matrici molto simili) compare talmente spesso da venire usata come matrice ``test'' per verificare l'efficienza e la precisione dei vari metodi per la risoluzione di sistemi lineari (vedi pag. 283 del libro di testo per un esempio in due dimensioni). Ovviamente, per migliorare l'efficienza, i vari metodi vanno specializzati per tenere conto della particolare struttura della matrice (gli elementi non nulli sono solo 3n-2, e quindi molte locazioni di memoria e operazioni di calcolo possono essere risparmiate!).

Sono molte le proprietà di A_n che possono essere dedotte teoricamente. Si può ad esempio dimostrare che gli autovalori sono

$$\lambda_j=2+2\cos{j\pi\over n+1},\hspace{1cm}j=1,\ldots,n,$$

e da essi possiamo ricavare il numero di condizionamento

$$\mu_2( A_n) = {\lambda_{\max}\over\lambda_{\min}} = {2+2\cos... ...i\over n+1} \over 2-2\cos{\pi\over n+1}} \approx {n^2\over 2}$$

sfruttando lo sviluppo in serie di $\cos x=1-x^2/2+o(x^2)$.

È un esercizio interessante cercare di dedurre teoricamente le proprietà dei vari metodi quando sono applicati a questa matrice, e verificare sperimentalmente il loro comportamento.

La decomposizione in forma A_n=L_nU_n di questo tipo di matrici è nota:

$$L_n= \left[ \matrix{ 1 & & & & \cr -{1\over 2} & 1 & & & \... ...& \cr & & & \ddots & -1 \cr & & & & {n+1\over n} } \right]$$

da cui si ricava immediatamente che $\det A_n=\prod_i u_{ii} = n+1$. Per quanto riguarda la matrice $\hat{L}_n$ della decomposizione $\hat{L}\hat{L}^T$ di Cholesky, si ha

$$\hat{L}_n= \left[ \matrix{ \sqrt{2} & & & & & \cr -\sqrt{1... ...\cr & & &&-\sqrt{n-1\over n} & \sqrt{n+1\over n} } \right] .$$

Si noti che in entrambe le decomposizioni la struttura di A_n si mantiene: il numero di elementi non nulli in $L_n, U_n, \hat{L}_n$ è sempre dell'ordine di n. Questa proprietà di sparsità si perde invece nelle inverse, che sono sempre ``piene''. Ecco alcuni esempi di inverse per n=4,5:

$$A_4^{-1}={1\over 5} \left[ \matrix{ 4 & 3 & 2 & 1 \cr 3 & 6 ... ...r 0.4 & 0.8 & 1.2 & 0.6 \cr 0.2 & 0.4 & 0.6 & 0.8 } \right]$$

$$A_5^{-1}={1\over 6} \left[ \matrix{ 5 & 4 & 3 & 2 & 1 \cr 4... ... & 0.667 \cr 0.167 & 0.333 & 0.5 & 0.667 & 0.833 } \right] .$$

Per n=4, queste sono le matrici intermedie A^(k) che si ottengono applicando l'eliminazione di Gauss:

$$A^{(2)}= \left[ \matrix{ 2.0 &-1.0 & 0.0 & 0.0 \cr 0.0 & 1.5 ... ...& 0.00 & 1.33 &-1.00 \cr 0.00 & 0.00 &-1.00 & 2.00 } \right]$$

$$A^{(4)}= \left[ \matrix{ 2.00 &-1.00 & 0.00 & 0.00 \cr 0.00 &... ...& 1.33 &-1.00 \cr 0.00 & 0.00 & 0.00 & 1.25 } \right] = U_4 .$$