Salta al contenuto
Note per Studenti Esercizio - sistema 2x2 stiff con autovalori -1 e -100 e limite di stabilità

Esercizio - sistema 2x2 stiff con autovalori -1 e -100 e limite 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 →, Autovalori e autovettoriUn autovettore è un vettore non nullo che una funzione lineare manda in un suo multiplo; si trovano gli autovalori come radici del polinomio caratteristico det(A − λI) e gli autovettori come nucleo di A − λI. Matrici simili hanno gli stessi autovalori.Autovalori e autovettori →.

Testo. Due serbatoi termici accoppiati hanno temperature (scostamenti) y=(y1,y2)T\mathbf y=(y_1,y_2)^T che obbediscono a y′=Ay\mathbf y'=A\mathbf y con A=(−50,549,549,5−50,5),y(0)=(10).A=\begin{pmatrix}-50{,}5&49{,}5\\49{,}5&-50{,}5\end{pmatrix},\qquad\mathbf y(0)=\begin{pmatrix}1\\0\end{pmatrix}.

  1. Calcolare autovalori e autovettori e la soluzione esatta. Qual è lo stiffness ratio?
  2. Determinare il massimo passo hh per la stabilità di Eulero esplicito.
  3. Con h=0,025h=0{,}025 eseguire tre passi di EE e osservare il comportamento.
  4. Con h=0,1h=0{,}1 eseguire tre passi di Eulero implicito, risolvendo il sistema, e confrontare con la soluzione esatta.

Soluzione

1. Autovalori. A=−50 I+(−0,549,549,5−0,5)A=-50\,I+\begin{pmatrix}-0{,}5&49{,}5\\49{,}5&-0{,}5\end{pmatrix}; la matrice (abba)\begin{pmatrix}a&b\\b&a\end{pmatrix} ha autovalori a±ba\pm b con autovettori (1,±1)(1,\pm1): a=−50,5a=-50{,}5, b=49,5b=49{,}5 dà λ1=−50,5+49,5=−1\lambda_1=-50{,}5+49{,}5=-1 con v1=(1,1)\mathbf v_1=(1,1) e λ2=−50,5−49,5=−100\lambda_2=-50{,}5-49{,}5=-100 con v2=(1,−1)\mathbf v_2=(1,-1). Si scrive (1,0)=12(1,1)+12(1,−1)(1,0)=\frac12(1,1)+\frac12(1,-1), quindi y(t)=12e−t(11)+12e−100t(1−1).\mathbf y(t)=\tfrac12e^{-t}\begin{pmatrix}1\\1\end{pmatrix}+\tfrac12e^{-100t}\begin{pmatrix}1\\-1\end{pmatrix}. Le due temperature si avvicinano a un valore comune (componente lenta, scala di tempo 11 s) e la loro differenza svanisce in ≈0,05\approx0{,}05 s (componente veloce, scala 0,010{,}01 s). Stiffness ratio: 1001=100\frac{100}{1}=100.

2. Stabilità di EE. Servono ∣1+hλi∣<1|1+h\lambda_i|<1 per entrambi: ∣1−h∣<1  ⟺  h<2|1-h|<1\iff h<2 e ∣1−100h∣<1  ⟺  h<0,02|1-100h|<1\iff h<0{,}02. Il vincolo più stringente è h<0,02\mathbf{h<0{,}02}, imposto dall'autovalore −100-100, anche se per t>0,05t>0{,}05 quella componente è già sparita e ci vorrebbero passi di 0,10{,}1 s per seguire e−te^{-t}.

3. EE con h=0,025h=0{,}025. I+hA=(1−1,26251,23751,23751−1,2625)=(−0,26251,23751,2375−0,2625)I+hA=\begin{pmatrix}1-1{,}2625&1{,}2375\\1{,}2375&1-1{,}2625\end{pmatrix}=\begin{pmatrix}-0{,}2625&1{,}2375\\1{,}2375&-0{,}2625\end{pmatrix}, con autovalori 1−0,025=0,9751-0{,}025=0{,}975 e 1−2,5=−1,51-2{,}5=-1{,}5.

  • u1=(−0,2625; 1,2375)\mathbf u_1=(-0{,}2625;\ 1{,}2375);
  • u2=(−0,2625⋅(−0,2625)+1,2375⋅1,2375; 1,2375⋅(−0,2625)+(−0,2625)⋅1,2375)=(1,6003; −0,6497)\mathbf u_2=\big(-0{,}2625\cdot(-0{,}2625)+1{,}2375\cdot1{,}2375;\ 1{,}2375\cdot(-0{,}2625)+(-0{,}2625)\cdot1{,}2375\big)=(1{,}6003;\ -0{,}6497);
  • u3=(−1,2241; 2,1509)\mathbf u_3=(-1{,}2241;\ 2{,}1509). La componente veloce, che nell'esatta è 12e−100t\frac12e^{-100t} e al primo passo (t=0,025t=0{,}025) vale già solo 12e−2,5=0,041\frac12e^{-2{,}5}=0{,}041, nello schema viene moltiplicata per −1,5-1{,}5 a ogni passo e cresce: u1−u22\frac{u_1-u_2}2 vale 0,5; −0,75; 1,125; −1,68750{,}5;\ -0{,}75;\ 1{,}125;\ -1{,}6875 per n=0,1,2,3n=0,1,2,3. La soluzione numerica oscilla con ampiezza crescente: a t=1t=1 (4040 passi) vale ≈5,5⋅106\approx5{,}5\cdot10^6. Per confronto con h=0,01h=0{,}01 (I+hA=(0,4950,4950,4950,495)I+hA=\begin{pmatrix}0{,}495&0{,}495\\0{,}495&0{,}495\end{pmatrix}, autovalori 0,990{,}99 e 00): u1=(0,495;0,495)\mathbf u_1=(0{,}495;0{,}495), u2=(0,49005;0,49005)\mathbf u_2=(0{,}49005;0{,}49005): la componente veloce è annullata subito, la lenta decade come 0,99n≈e−0,01n0{,}99^n\approx e^{-0{,}01n} ✓.

4. EI con h=0,1h=0{,}1. (I−hA)un+1=un(I-hA)\mathbf u_{n+1}=\mathbf u_n, con I−hA=(6,05−4,95−4,956,05)I-hA=\begin{pmatrix}6{,}05&-4{,}95\\-4{,}95&6{,}05\end{pmatrix}, det⁡=6,052−4,952=(6,05−4,95)(6,05+4,95)=1,1⋅11=12,1\det=6{,}05^2-4{,}95^2=(6{,}05-4{,}95)(6{,}05+4{,}95)=1{,}1\cdot11=12{,}1 e inversa 112,1(6,054,954,956,05)\frac1{12{,}1}\begin{pmatrix}6{,}05&4{,}95\\4{,}95&6{,}05\end{pmatrix} (autovalori della matrice di passo: 11+0,1=0,909\frac1{1+0{,}1}=0{,}909 e 11+10=0,0909\frac1{1+10}=0{,}0909).

  • u1=112,1(6,05; 4,95)=(0,5; 0,409091)\mathbf u_1=\frac1{12{,}1}(6{,}05;\ 4{,}95)=(0{,}5;\ 0{,}409091);
  • u2=112,1(6,05⋅0,5+4,95⋅0,409091; 4,95⋅0,5+6,05⋅0,409091)=112,1(5,05; 4,95)=(0,417355; 0,409091)\mathbf u_2=\frac1{12{,}1}(6{,}05\cdot0{,}5+4{,}95\cdot0{,}409091;\ 4{,}95\cdot0{,}5+6{,}05\cdot0{,}409091)=\frac1{12{,}1}(5{,}05;\ 4{,}95)=(0{,}417355;\ 0{,}409091);
  • u3=(0,376033; 0,375282)\mathbf u_3=(0{,}376033;\ 0{,}375282). Valori esatti: t=0,1t=0{,}1: (0,452441; 0,452396)(0{,}452441;\ 0{,}452396); t=0,3t=0{,}3: 12e−0,3=0,370409\frac12e^{-0{,}3}=0{,}370409 per entrambe le componenti (la componente veloce è 12e−30≈0\frac12e^{-30}\approx0). Errore di EI a t=0,3t=0{,}3 su y1y_1: 0,370409−0,376033=−5,6⋅10−30{,}370409-0{,}376033=-5{,}6\cdot10^{-3} (circa 1,5%1{,}5\%): è l'errore di ordine 1 della componente lenta, non un problema di stabilità. Nessuna esplosione: la componente veloce, moltiplicata per 0,09090{,}0909 a ogni passo, scompare dopo due passi, come nella soluzione vera.

Risultati: λ=−1, −100\lambda=-1,\ -100 (stiffness ratio 100100); EE stabile solo per h<0,02h<0{,}02 (h=0,025h=0{,}025: fattore −1,5-1{,}5 sulla componente veloce, diverge); EI con h=0,1h=0{,}1: u3=(0,376; 0,375)\mathbf u_3=(0{,}376;\ 0{,}375), errore 1,5%1{,}5\%.

Errori tipici. Scegliere hh in base a λ=−1\lambda=-1 (la componente che interessa) invece che a λ=−100\lambda=-100; calcolare l'inversa di I−hAI-hA a ogni passo invece di fattorizzarla una volta; dimenticare l'autovalore più grande in modulo.

Versione ripasso

(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 →)

A=(−50,549,549,5−50,5)A=\begin{pmatrix}-50{,}5&49{,}5\\49{,}5&-50{,}5\end{pmatrix}, y(0)=(1,0)\mathbf y(0)=(1,0).

  • λ=a±b=−1 (v=(1,1))\lambda=a\pm b=-1\ (\mathbf v=(1,1)), −100 (v=(1,−1))-100\ (\mathbf v=(1,-1)); y=12e−t(1,1)+12e−100t(1,−1)\mathbf y=\frac12e^{-t}(1,1)+\frac12e^{-100t}(1,-1); stiffness ratio 100100.

  • EE: ∣1−h∣<1|1-h|<1 e ∣1−100h∣<1⇒h<0,02|1-100h|<1\Rightarrow h<0{,}02 (vincolo dall'autovalore −100-100).

  • h=0,025h=0{,}025: I+hA=(−0,26251,23751,2375−0,2625)I+hA=\begin{pmatrix}-0{,}2625&1{,}2375\\1{,}2375&-0{,}2625\end{pmatrix}, autovalori 0,9750{,}975 e −1,5-1{,}5; u1=(−0,2625;1,2375)\mathbf u_1=(-0{,}2625;1{,}2375), u2=(1,6003;−0,6497)\mathbf u_2=(1{,}6003;-0{,}6497), u3=(−1,2241;2,1509)\mathbf u_3=(-1{,}2241;2{,}1509): diverge. h=0,01h=0{,}01: matrice con tutti gli elementi 0,4950{,}495, autovalori 0,990{,}99 e 00.

  • EI h=0,1h=0{,}1: I−hA=(6,05−4,95−4,956,05)I-hA=\begin{pmatrix}6{,}05&-4{,}95\\-4{,}95&6{,}05\end{pmatrix}, det⁡=12,1\det=12{,}1; u1=(0,5;0,409091)\mathbf u_1=(0{,}5;0{,}409091), u2=(0,417355;0,409091)\mathbf u_2=(0{,}417355;0{,}409091), u3=(0,376033;0,375282)\mathbf u_3=(0{,}376033;0{,}375282); esatto t=0,3t=0{,}3: 0,3704090{,}370409, errore 1,5%1{,}5\% (ordine 1).

  • Decomposizione: (1,0)=12(1,1)+12(1,−1)(1,0)=\frac12(1,1)+\frac12(1,-1); (abba)\begin{pmatrix}a&b\\b&a\end{pmatrix} ha autovalori a±ba\pm b. Componente veloce di EE: u1−u22=0,5; −0,75; 1,125; −1,6875\frac{u_1-u_2}2=0{,}5;\,-0{,}75;\,1{,}125;\,-1{,}6875 (fattore −1,5-1{,}5).

  • det⁡(I−hA)=1,1⋅11=12,1\det(I-hA)=1{,}1\cdot11=12{,}1; EI: autovalori del passo 0,9090{,}909 e 0,09090{,}0909; esatto t=0,3t=0{,}3: 12e−0,3=0,370409\frac12e^{-0{,}3}=0{,}370409.

Errori tipici: hh scelto su λ=−1\lambda=-1; inversa a ogni passo; autovalore massimo dimenticato.

Teoria collegata