Esercizio - conduzione del calore in una barra con Jacobi e Gauss-Seidel

Teoria: Metodi iterativi - Jacobi, Gauss-Seidel e convergenzaSi spezza $A=M-N$ con $M$ facile da invertire e si itera $M\mathbf x^{(k+1)}=N\mathbf x^{(k)}+\mathbf b$, cioè $\mathbf x^{(k+1)}=B\mathbf x^{(k)}+\mathbf c$ con $B=M^{-1}N=I-M^{-1}A$. Con $A=D+L+U$ (diagonale, parte strettamente inferiore e superiore): Jacobi $M=D$, $x_i^{(k+1)}=\frac1{a_{ii}}\big(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\big)$; Gauss-Seidel $M=D+L$, usa subito i valori nuovi $x_j^{(k+1)}$, $j<i$. Convergenza per ogni $\mathbf x^{(0)}$ se e solo se $\rho(B)<1$ (D); basta $|B|<1$. Matrice a dominanza diagonale stretta per righe: convergono entrambi; $A$ SPD: converge Gauss-Seidel; per matrici tridiagonali $\rho(B_{GS})=\rho(B_J)^2$. Errore: $\mathbf e^{(k)}=B^k\mathbf e^{(0)}$, $|\mathbf e^{(k)}|\le|B|^k|\mathbf e^{(0)}|$ e $|\mathbf x-\mathbf x^{(k)}|\le\frac{|B|}{1-|B|}|\mathbf x^{(k)}-\mathbf x^{(k-1)}|$. Per guadagnare $d$ cifre servono $k\approx\frac d{-\log_{10}\rho(B)}$ iterazioni.Metodi iterativi - Jacobi, Gauss-Seidel e convergenza →; la stessa matrice è in Esercizio - condizionamento della matrice di rigidezza di una catena di molle.

Testo. Una barra di lunghezza $L=1$ m con generazione interna di calore uniforme ha, a regime, temperatura $T(x)$ che soddisfa $-T''=s$ con $s=8\ \text{K/m}^2$ e T(0)=T(L)=0\,^\circC (scostamento dalla temperatura degli estremi). La soluzione esatta è $T(x)=4x(1-x)$. Si discretizza con $n=3$ nodi interni a passo $h=\frac L4=0{,}25$ m: la formula alle differenze centrate dà $$2T_i-T_{i-1}-T_{i+1}=h^2s=0{,}5\ \text{K}\qquad(i=1,2,3,\ T_0=T_4=0).$$

1. Scrivere il sistema e verificare che la soluzione esatta del sistema è $T_i=T(x_i)$.
2. Eseguire tre iterazioni di Jacobi e tre di Gauss-Seidel da $\mathbf T^{(0)}=\mathbf 0$ e calcolare gli errori.
3. Calcolare $\rho(B_J)$ e $\rho(B_{GS})$ e il numero di iterazioni per ridurre l'errore di un fattore $10^{-3}$.
4. Quante iterazioni servirebbero con $n=100$? Conviene un metodo iterativo?

Soluzione

1. Sistema. $K\mathbf T=0{,}5\,\mathbf 1$ con $K=\begin{pmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{pmatrix}$. I nodi sono $x=0{,}25;\,0{,}5;\,0{,}75$ e $T(x_i)=4x_i(1-x_i)=(0{,}75;\ 1;\ 0{,}75)$ K. Verifica: $2\cdot0{,}75-1=0{,}5$ ✓; $2\cdot1-0{,}75-0{,}75=0{,}5$ ✓; $2\cdot0{,}75-1=0{,}5$ ✓. (Le differenze centrate sono esatte per i polinomi di secondo grado, quindi la soluzione discreta coincide con quella continua nei nodi.) Esplicitando le incognite: $T_1=\frac{0{,}5+T_2}2$, $T_2=\frac{0{,}5+T_1+T_3}2$, $T_3=\frac{0{,}5+T_2}2$.

2. Iterazioni (temperature in K, errore in norma $\infty$ rispetto a $(0{,}75;1;0{,}75)$). Jacobi:

• $k=1$: $T_1=\frac{0{,}5+0}2=0{,}25$; $T_2=\frac{0{,}5}2=0{,}25$; $T_3=0{,}25$. Errore $0{,}75$.
• $k=2$: $T_1=\frac{0{,}5+0{,}25}2=0{,}375$; $T_2=\frac{0{,}5+0{,}25+0{,}25}2=0{,}5$; $T_3=0{,}375$. Errore $0{,}5$.
• $k=3$: $T_1=\frac{0{,}5+0{,}5}2=0{,}5$; $T_2=\frac{0{,}5+0{,}375+0{,}375}2=0{,}625$; $T_3=0{,}5$. Errore $0{,}375$.

Gauss-Seidel:

• $k=1$: $T_1=\frac{0{,}5}2=0{,}25$; $T_2=\frac{0{,}5+0{,}25+0}2=0{,}375$; $T_3=\frac{0{,}5+0{,}375}2=0{,}4375$. Errore $0{,}625$.
• $k=2$: $T_1=\frac{0{,}5+0{,}375}2=0{,}4375$; $T_2=\frac{0{,}5+0{,}4375+0{,}4375}2=0{,}6875$; $T_3=\frac{0{,}5+0{,}6875}2=0{,}59375$. Errore $0{,}3125$.
• $k=3$: $T_1=\frac{0{,}5+0{,}6875}2=0{,}59375$; $T_2=\frac{0{,}5+0{,}59375+0{,}59375}2=0{,}84375$; $T_3=\frac{0{,}5+0{,}84375}2=0{,}671875$. Errore $0{,}15625$.

Dopo 3 iterazioni la temperatura al centro è $0{,}625$ K (Jacobi) o $0{,}844$ K (Gauss-Seidel) contro $1$ K: la convergenza è lenta. Il rapporto fra errori successivi di Gauss-Seidel è esattamente $0{,}5$ da $k=2$ in poi.

3. Raggi spettrali. $B_J=I-\frac12K=\begin{pmatrix}0&\frac12&0\\\frac12&0&\frac12\\0&\frac12&0\end{pmatrix}$ ha autovalori $\cos\frac{j\pi}4$: $\pm\frac{\sqrt2}2,\ 0$, quindi $$\rho(B_J)=\cos\frac\pi4=0{,}7071,\qquad\rho(B_{GS})=\rho(B_J)^2=0{,}5\quad(\text{matrice tridiagonale}).$$ Per ridurre l'errore di $10^{-3}$: Jacobi $k\ge\frac3{-\log_{10}0{,}7071}=\frac3{0{,}1505}=19{,}9\Rightarrow20$ iterazioni; Gauss-Seidel $k\ge\frac3{-\log_{10}0{,}5}=\frac3{0{,}301}=9{,}97\Rightarrow10$ iterazioni. Gauss-Seidel dimezza il numero di iterazioni.

4. Griglia fine. In generale $\rho(B_J)=\cos\frac\pi{n+1}$. Per $n=100$: $\rho_J=0{,}99952$, $\rho_{GS}=0{,}99903$; per guadagnare 6 cifre servono $\frac6{-\log_{10}\rho}=28\,554$ iterazioni di Jacobi e $14\,277$ di Gauss-Seidel. Ogni iterazione costa $\approx5n$ operazioni: in totale $\approx7\cdot10^6$ flop per Gauss-Seidel, contro $\approx8n=800$ flop per risolvere direttamente il sistema tridiagonale (Cholesky o LU per matrici tridiagonali: $O(n)$). Per una matrice a banda stretta conviene il metodo diretto. Gli iterativi servono per matrici sparse di grandi dimensioni non a banda stretta (griglie 2D o 3D), e anche lì con accelerazioni (gradiente coniugato, multigriglia), perché $\rho\to1$ al raffinare della griglia.

Risultati: errori dopo 3 iterazioni: Jacobi $0{,}375$ K, Gauss-Seidel $0{,}156$ K; $\rho_J=0{,}7071$, $\rho_{GS}=0{,}5$; 20 e 10 iterazioni per $10^{-3}$; con $n=100$: $28\,554$ e $14\,277$ iterazioni per $10^{-6}$.

Errori tipici. Usare $\rho_{GS}=\rho_J$ (vale il quadrato); dimenticare che $h^2s$ è il termine noto di ogni equazione e non solo del nodo centrale; concludere che un metodo iterativo è sempre più economico di uno diretto.