Esercizio - passo massimo di Eulero esplicito per un'equazione test e confronto con i metodi impliciti

Teoria: 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 →, 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. La tensione $v$ su un condensatore che si scarica attraverso una resistenza obbedisce a $v'=-80\,v$ ($t$ in secondi, $v$ in volt), con $v(0)=5$ V (costante di tempo $\tau=\frac1{80}=12{,}5$ ms).

1. Ricavare il fattore di amplificazione di Eulero esplicito (EE), Eulero implicito (EI) e Crank-Nicolson (CN) e il passo massimo $h$ per la stabilità di EE.
2. Calcolare i fattori per $h=0{,}01;\ 0{,}02;\ 0{,}025;\ 0{,}03;\ 0{,}05$ s e classificare il comportamento di EE.
3. Con $h=0{,}05$ s calcolare 4 passi dei tre metodi (fino a $t=0{,}2$ s) e confrontare con $v(0{,}2)=5\,e^{-16}=5{,}6\cdot10^{-7}$ V.
4. (Teoria) Con $h=0{,}02$ s si può dire che EE è «accettabile»? Fino a quale $h$ la soluzione di EE è anche monotona?

Soluzione

1. Fattori. Con $z=h\lambda=-80h$: EE $R=1+z=1-80h$; EI $R=\frac1{1-z}=\frac1{1+80h}$; CN $R=\frac{1+z/2}{1-z/2}=\frac{1-40h}{1+40h}$. EE è stabile se $|1-80h|<1$, cioè $0<h<\frac2{80}=0{,}025$ s. EI e CN sono stabili per ogni $h>0$.

2. Tabella dei fattori.

3. Con $h=0{,}05$. $v_n=v_0R^n$, $v_0=5$:

• EE ($R=-3$): $-15;\ 45;\ -135;\ 405$ V. La soluzione numerica esplode e cambia segno a ogni passo: fisicamente assurda (una tensione di $405$ V da un condensatore a $5$ V!).
• EI ($R=0{,}2$): $1;\ 0{,}2;\ 0{,}04;\ 0{,}008$ V.
• CN ($R=-\frac13$): $-1{,}667;\ 0{,}556;\ -0{,}185;\ 0{,}0617$ V. Valore esatto a $t=0{,}2$ s: $5{,}6\cdot10^{-7}$ V. EI e CN sono molto più vicini di EE ($405$ V); CN oscilla ma decade, EI decade monotono (è più piatto del vero, ma non oscilla). L'accuratezza con $h=0{,}05$ è comunque scarsa perché $|\lambda|h=4$ è grande rispetto a un'accuratezza reale ($\tau=0{,}0125$ s è minore di $h$): se si volesse seguire la scarica si userebbe $h\lesssim\tau/5$; ma la stabilità di EI e CN è garantita.

4. Teoria. Con $h=0{,}02$ il fattore è $-0{,}6$: la soluzione decade (stabile) ma oscilla di segno ($5;\ -3;\ 1{,}8;\ -1{,}08;\dots$), mentre la vera è positiva e monotona. Per avere decadimento monotono serve $R\ge0$, cioè $1-80h\ge0\iff h\le\frac1{80}=0{,}0125$ s $=\tau$. Quindi: stabilità $h<2\tau=0{,}025$ s; assenza di oscillazioni spurie $h\le\tau=0{,}0125$ s. Per un'accuratezza ragionevole (errore dell'ordine del 10% per passo) $h\approx\frac\tau5=2{,}5$ ms. Qui la stabilità non è il vincolo (si usa spesso $h\ll2\tau$ per accuratezza), ma lo diventa se il sistema contiene anche una componente lenta e una molto veloce (problema stiff).

Risultati: $h_{\max}^{EE}=\frac2{80}=0{,}025$ s; con $h=0{,}05$ s: EE $405$ V al quarto passo (instabile), EI $0{,}008$ V, CN $0{,}062$ V.

Errori tipici. Confondere $h<\frac2{|\lambda|}$ (stabilità) con $h<\frac1{|\lambda|}$ (assenza di oscillazioni); pensare che «decade» significhi «accurato»; usare $|\lambda|h$ con $\lambda$ positivo.