Esercizio - massa-molla-smorzatore ridotto al primo ordine e passo di stabilità

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 →, Oscillatore massa-molla-smorzatore e transitoriMassa m, molla k e smorzatore viscoso c con forza esterna F(t): m x'' + c x' + k x = F(t), cioè x'' + (c/m) x' + (k/m) x = F/m. Grandezze: ω0 = √(k/m), smorzamento critico c_cr = 2√(km), ζ = c/c_cr. Libero: c > c_cr due esponenziali negativi (sovrasmorzato, nessuna oscillazione), c = c_cr (c1 + c2 t) e^(−ω0 t), c < c_cr e^(−ct/2m)(c1 cos ωd t + c2 sin ωd t) con ωd = ω0√(1−ζ²); con c > 0 ogni moto libero tende a 0 (transitorio). Forzato F0 cos Ωt con c > 0: x = transitorio + regime permanente X cos(Ωt − φ), X = F0/√((k−mΩ²)² + c²Ω²), tan φ = cΩ/(k−mΩ²); massimo di X per Ω² = k/m − c²/(2m²). Primo ordine m v' + c v = F: v = F/c + (v0 − F/c) e^(−ct/m), costante di tempo τ = m/c.Oscillatore massa-molla-smorzatore e transitori →.

Testo. Il moto libero di un sistema massa-molla-smorzatore con $m=1$ kg, $c=2$ N s/m, $k=101$ N/m è $x''+2x'+101x=0$, con $x(0)=1$ m, $x'(0)=0$.

1. Ridurre l'equazione a un sistema del primo ordine e trovare gli autovalori.
2. Eseguire tre passi di Eulero esplicito con $h=0{,}01$ s e confrontare con la soluzione esatta $x(t)=e^{-t}\big(\cos10t+0{,}1\sin10t\big)$.
3. Trovare il massimo passo $h$ per cui EE è assolutamente stabile. Confrontare con la regola $h<\frac2{|\lambda|}$.
4. Cosa succede con EI e CN per $h=0{,}1$ s?

Soluzione

1. Sistema. Con $v=x'$: $x'=v$, $v'=-101x-2v$, $$\begin{pmatrix}x\\v\end{pmatrix}'=\begin{pmatrix}0&1\\-101&-2\end{pmatrix}\begin{pmatrix}x\\v\end{pmatrix},\qquad\mathbf y(0)=\begin{pmatrix}1\\0\end{pmatrix}.$$ Polinomio caratteristico $\lambda^2+2\lambda+101=0$, quindi $\lambda=\dfrac{-2\pm\sqrt{4-404}}{2}=-1\pm10i$ (sistema sottosmorzato: $\mu=-1$ dà il decadimento $e^{-t}$, $\omega=10$ rad/s l'oscillazione). $|\lambda|=\sqrt{101}=10{,}05$.

2. Tre passi di EE con $h=0{,}01$. $\mathbf u_{n+1}=\mathbf u_n+h\begin{pmatrix}v_n\\-101x_n-2v_n\end{pmatrix}$.

• $n=0$: $x_1=1+0{,}01\cdot0=1$; $v_1=0+0{,}01\,(-101\cdot1-0)=-1{,}01$.
• $n=1$: $x_2=1+0{,}01\cdot(-1{,}01)=0{,}9899$; $v_2=-1{,}01+0{,}01\,(-101\cdot1-2\cdot(-1{,}01))=-1{,}01+0{,}01\cdot(-98{,}98)=-1{,}9998$.
• $n=2$: $x_3=0{,}9899+0{,}01\cdot(-1{,}9998)=0{,}969902$; $v_3=-1{,}9998+0{,}01\,(-101\cdot0{,}9899-2\cdot(-1{,}9998))=-1{,}9998+0{,}01\cdot(-99{,}9799+3{,}9996)=-1{,}9998-0{,}959803=-2{,}959603$.

Valori esatti: $x(0{,}01)=0{,}99499$, $x(0{,}02)=0{,}98013$, $x(0{,}03)=0{,}95578$. Errori: $x_1$ $+5\cdot10^{-3}$, $x_2$ $+9{,}8\cdot10^{-3}$, $x_3$ $+1{,}4\cdot10^{-2}$ (EE sovrastima: la molla non «frena» abbastanza nei primi istanti).

3. Stabilità di EE. Si impone $|1+h\lambda|<1$ con $\lambda=-1+10i$: $$|1+h\lambda|^2=(1-h)^2+(10h)^2=1-2h+101h^2<1\iff h\,(101h-2)<0\iff h<\frac2{101}=0{,}0198\ \text{s}.$$ In generale per $\lambda=\mu\pm i\omega$ (con $\mu<0$): $h<\dfrac{2|\mu|}{\mu^2+\omega^2}=\dfrac{2|\mathrm{Re}\,\lambda|}{|\lambda|^2}$. Controllo numerico del raggio spettrale di $I+hA$: $h=0{,}01$: $0{,}99504$; $h=0{,}0198$: $0{,}999998$; $h=0{,}02$: $1{,}0002$ (instabile); $h=0{,}05$: $1{,}0735$. Confronto con la regola $h<\frac2{|\lambda|}$: darebbe $\frac2{10{,}05}=0{,}199$ s, dieci volte più grande di quello giusto. La regola $\frac2{|\lambda|}$ vale solo per $\lambda$ reale negativo; qui il decadimento è lentissimo rispetto all'oscillazione ($|\mu|=1\ll\omega=10$) e il disco di stabilità di EE è tangente all'asse immaginario, quindi un autovalore quasi immaginario lo lascia quasi subito.

4. EI e CN con $h=0{,}1$ s. I raggi spettrali dei fattori di passo sono $\rho_{EI}=0{,}6727<1$ e $\rho_{CN}=0{,}9231<1$: sono stabili (A-stabili: $\mathrm{Re}\,\lambda<0$). Con $h=0{,}1$, però, ogni passo copre circa $\frac16$ del periodo ($T=\frac{2\pi}{10}=0{,}63$ s, circa $6$ passi per periodo): l'accuratezza è scarsa, EI smorza troppo (decadimento numerico $0{,}67^n$ contro $e^{-0{,}1n}=0{,}905^n$ del vero), CN è più vicina ($0{,}923^n$ contro $0{,}905^n$). La stabilità non basta: servono almeno $\approx20$–$30$ passi per periodo per l'accuratezza.

Risultati: $\lambda=-1\pm10i$; EE con $h=0{,}01$: $x_3=0{,}9699$ (esatto $0{,}9558$); $\mathbf{h_{\max}=\frac{2|\mathrm{Re}\lambda|}{|\lambda|^2}=\frac2{101}=0{,}0198}$ s (non $0{,}199$); EI e CN stabili per ogni $h$.

Errori tipici. Usare $\frac2{|\lambda|}$ con $\lambda$ complesso; dimenticare il termine $-2v$ nell'equazione di $v$; confondere stabilità (si esplode o no) con accuratezza (quanto è vicino).