Esercizio - Eulero esplicito a mano e confronto con la soluzione esatta del raffreddamento

Teoria: Metodi di Eulero esplicito, implicito e Crank-NicolsonProblema di Cauchy $y'=f(t,y)$, $y(t_0)=y_0$, con $f$ lipschitziana in $y$ (costante $L$): ben posto. Nodi $t_n=t_0+nh$, $u_n\approx y(t_n)$. Eulero esplicito $u_{n+1}=u_n+h,f(t_n,u_n)$; Eulero implicito $u_{n+1}=u_n+h,f(t_{n+1},u_{n+1})$ (equazione da risolvere a ogni passo: punto fisso se $hL<1$, Newton, o formula chiusa se $f$ è lineare); Crank-Nicolson $u_{n+1}=u_n+\frac h2[f(t_n,u_n)+f(t_{n+1},u_{n+1})]$ (regola dei trapezi). Errore di troncamento locale $\tau_n=\frac{y(t_{n+1})-y(t_n)}h-\Phi$: $\frac h2y''$ (EE), $-\frac h2y''$ (EI), $-\frac{h^2}{12}y'''$ (CN). Consistenza + stabilità $\Rightarrow$ convergenza: ordine 1 per EE e EI, ordine 2 per CN, cioè errore globale $O(h)$ e $O(h^2)$; per EE $|e_n|\le\frac{Mh}{2L}\big(e^{L(t_n-t_0)}-1\big)$ con $M=\max|y''|$ (D). Dimezzando $h$ l'errore si divide per 2 (EE, EI) o per 4 (CN).Metodi di Eulero esplicito, implicito e Crank-Nicolson →.

Testo. Un pezzo metallico a 90\,^\circC viene lasciato raffreddare in un ambiente a 20\,^\circC. Per la legge di Newton $T'=-k\,(T-20)$ con $k=0{,}4\ \text{min}^{-1}$ ($t$ in minuti, $T$ in $^\circ$C).

1. Eseguire tre passi di Eulero esplicito con $h=1$ min e confrontare con la soluzione esatta $T(t)=20+70\,e^{-0{,}4t}$.
2. Ripetere fino a $t=3$ min con $h=0{,}5$ min e confrontare gli errori.
3. Calcolare la maggiorazione teorica dell'errore globale in $t=3$ per $h=1$.
4. Per quali valori di $h$ lo schema decresce monotonicamente verso 20\,^\circC? Cosa succede per $h=6$ min?

Soluzione

1. Tre passi con $h=1$. Lo schema è $T_{n+1}=T_n+h\,f(t_n,T_n)=T_n-0{,}4\,(T_n-20)$.

• $T_1=90-0{,}4\cdot70=90-28=62$;
• $T_2=62-0{,}4\cdot42=62-16{,}8=45{,}2$;
• $T_3=45{,}2-0{,}4\cdot25{,}2=45{,}2-10{,}08=35{,}12$.

Valori esatti: $T(1)=20+70\cdot0{,}670320=66{,}922$; $T(2)=20+70\cdot0{,}449329=51{,}453$; $T(3)=20+70\cdot0{,}301194=41{,}084$. Errori $T(t_n)-T_n$: $4{,}92;\ 6{,}25;\ 5{,}96$ $^\circ$C. Eulero sottostima la temperatura (il pezzo si raffredda troppo in fretta): la tangente di una funzione convessa decrescente sta sotto la curva.

2. Con $h=0{,}5$. $T_{n+1}-20=(1-0{,}4\cdot0{,}5)(T_n-20)=0{,}8\,(T_n-20)$, quindi $T_n-20=70\cdot0{,}8^n$. Dopo 6 passi ($t=3$): $T_6=20+70\cdot0{,}8^6=20+70\cdot0{,}262144=38{,}350$. Errore $41{,}084-38{,}350=2{,}73$ $^\circ$C. Il rapporto con l'errore di $h=1$ è $\frac{5{,}96}{2{,}73}=2{,}18$: coerente con l'ordine $1$ (atteso $\approx2$).

3. Maggiorazione. $f(t,T)=-0{,}4\,(T-20)$ ha $L=0{,}4$ min$^{-1}$. $T''=0{,}16\cdot70\,e^{-0{,}4t}\le11{,}2$ $^\circ$C/min$^2$, quindi $M=11{,}2$. Il teorema dà $$|e_n|\le\frac{Mh}{2L}\Big(e^{L(t_n-t_0)}-1\Big)=\frac{11{,}2\cdot1}{0{,}8}\big(e^{1{,}2}-1\big)=14\cdot2{,}320=32{,}5\ ^\circ\text{C}.$$ È una maggiorazione corretta ($5{,}96\le32{,}5$) ma molto pessimistica (fattore $5{,}5$): il termine $e^{L(t-t_0)}$ prevede un'amplificazione degli errori che qui non avviene, perché la soluzione decresce verso l'equilibrio e gli errori vengono smorzati.

4. Comportamento in $h$. Lo schema è $T_{n+1}-20=(1-kh)(T_n-20)=(1-0{,}4h)(T_n-20)$.

• Per $0<h<2{,}5$ min il fattore è in $(0,1)$: la temperatura decresce monotonicamente a $20$.
• Per $2{,}5<h<5$ min il fattore è in $(-1,0)$: la temperatura oscilla attorno a $20$ ma tende comunque a $20$ (non fisico).
• Per $h=5$ min il fattore è $-1$ (oscillazione permanente) e per $h>5$ min è minore di $-1$: la soluzione numerica diverge. Con $h=6$: fattore $1-2{,}4=-1{,}4$: $T_1=90-2{,}4\cdot70=-78$, $T_2=-78+2{,}4\cdot98=157{,}2$, $T_3=157{,}2-2{,}4\cdot137{,}2=-172{,}1$ °C: temperature assurde, che crescono in modulo. Il limite $h<\frac2k=5$ min per non divergere è l'oggetto della nota Stabilità assoluta, problemi stiff e sistemi di equazioni differenzialiSul problema test $y'=\lambda y$ ($\mathrm{Re},\lambda<0$, soluzione che decade) i metodi danno $u_{n+1}=R(z),u_n$ con $z=h\lambda$: EE $R=1+z$, EI $R=\frac1{1-z}$, CN $R=\frac{1+z/2}{1-z/2}$. La regione di stabilità assoluta è ${z:|R(z)|<1}$: per EE il disco di centro $-1$ e raggio $1$ (per $\lambda<0$ reale: $h<\frac2{|\lambda|}$), per EI l'esterno del disco di centro $1$ e raggio $1$, per CN il semipiano $\mathrm{Re},z<0$. EI e CN sono A-stabili (nessun limite su $h$ per $\mathrm{Re},\lambda<0$); per $z\to-\infty$ è $R_{EI}\to0$ (smorza) e $R_{CN}\to-1$ (oscilla con smorzamento lento). Un problema è stiff quando ha componenti che decadono a velocità molto diverse (rapporto $\frac{\max|\mathrm{Re}\lambda_i|}{\min|\mathrm{Re}\lambda_i|}$ grande): l'accuratezza permetterebbe $h$ grande, ma la stabilità di un metodo esplicito impone $h<\frac2{\max|\lambda_i|}$. Sistemi $\mathbf y'=A\mathbf y$: EE $\mathbf u_{n+1}=(I+hA)\mathbf u_n$, EI $(I-hA)\mathbf u_{n+1}=\mathbf u_n$, CN $(I-\frac h2A)\mathbf u_{n+1}=(I+\frac h2A)\mathbf u_n$; stabilità se $h\lambda_i$ sta nella regione per ogni autovalore. Oscillatore ($\lambda=\pm i\omega$): $|R_{EE}|=\sqrt{1+\omega^2h^2}>1$ (energia cresce), $|R_{EI}|<1$ (smorza), $|R_{CN}|=1$ (conserva). Equazioni di ordine $m$ si riducono a sistemi del primo ordine.Stabilità assoluta, problemi stiff e sistemi di equazioni differenziali →.

Risultati: $T_1=62$, $T_2=45{,}2$, $T_3=35{,}12$ (errori $4{,}9;\ 6{,}3;\ 6{,}0$); con $h=0{,}5$: $T(3)\approx38{,}35$ (errore $2{,}7$); limite teorico $32{,}5$; stabile solo per $h<5$ min.

Errori tipici. Calcolare il passo con il valore sbagliato $f(t_{n+1},T_{n+1})$; dimenticare il $-20$ (usare $T'=-kT$); confondere «errore diminuisce con $h$» con «stabilità per $h$ grande».