Salta al contenuto
Note per Studenti Esercizio - oscillatore armonico con Eulero esplicito, implicito e Crank-Nicolson, energia

Esercizio - oscillatore armonico con Eulero esplicito, implicito e Crank-Nicolson, energia

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 →, Moto armonico sempliceMoto rettilineo con legge oraria x = A·sin(ωt + φ): oscilla tra −A e A con periodo T = 2π/ω. È caratterizzato da a = −ω²x, l'equazione differenziale dell'oscillatore armonico; ampiezza e fase si ricavano dalle condizioni iniziali.Moto armonico semplice →.

Testo. Una massa m=1m=1 kg collegata a una molla di costante k=16k=16 N/m, senza attrito, viene lasciata da ferma con spostamento x0=0,1x_0=0{,}1 m. Con v=x′v=x' il moto è y′=Ay\mathbf y'=A\mathbf y, y=(x,v)T\mathbf y=(x,v)^T, A=(01−160)A=\begin{pmatrix}0&1\\-16&0\end{pmatrix} (pulsazione ω=4\omega=4 rad/s, periodo T=2πω=1,571T=\frac{2\pi}\omega=1{,}571 s, energia conservata E=12v2+8x2=0,08E=\frac12v^2+8x^2=0{,}08 J). Si usa h=0,05h=0{,}05 s.

  1. Ricavare il fattore di amplificazione dell'energia a ogni passo per EE, EI e CN.
  2. Eseguire il primo passo con ciascun metodo e verificare l'energia.
  3. Calcolare l'energia dopo un periodo (n=31n=31 passi) e dopo n=40n=40.
  4. Quanti periodi si possono simulare con EE prima che l'energia raddoppi? Quale metodo conviene per un'integrazione su tempi lunghi?

Soluzione

1. Fattori. Gli autovalori di AA sono λ=±4i\lambda=\pm4i, z=hλ=±0,2 iz=h\lambda=\pm0{,}2\,i. Il modulo del fattore di amplificazione per passo e per l'energia (proporzionale al quadrato dell'ampiezza):

  • EE: ∣R∣2=∣1+0,2i∣2=1+0,04=1,04|R|^2=|1+0{,}2i|^2=1+0{,}04=1{,}04;
  • EI: ∣R∣2=1∣1−0,2i∣2=11,04=0,9615|R|^2=\frac1{|1-0{,}2i|^2}=\frac1{1{,}04}=0{,}9615;
  • CN: ∣R∣2=∣1+0,1i∣2∣1−0,1i∣2=1|R|^2=\frac{|1+0{,}1i|^2}{|1-0{,}1i|^2}=1. In generale 1±ω2h21\pm\omega^2h^2 con ω2h2=(4⋅0,05)2=0,04\omega^2h^2=(4\cdot0{,}05)^2=0{,}04.

2. Primo passo da y0=(0,1; 0)\mathbf y_0=(0{,}1;\ 0).

  • EE: y1=y0+hAy0=(0,1+0,05⋅0; 0+0,05⋅(−16⋅0,1))=(0,1; −0,08)\mathbf y_1=\mathbf y_0+hA\mathbf y_0=\big(0{,}1+0{,}05\cdot0;\ 0+0{,}05\cdot(-16\cdot0{,}1)\big)=(0{,}1;\ -0{,}08). Energia 12⋅0,082+8⋅0,12=0,0032+0,08=0,0832=1,04⋅0,08\frac12\cdot0{,}08^2+8\cdot0{,}1^2=0{,}0032+0{,}08=0{,}0832=1{,}04\cdot0{,}08 ✓.
  • EI: (I−hA)y1=y0(I-hA)\mathbf y_1=\mathbf y_0, I−hA=(1−0,050,81)I-hA=\begin{pmatrix}1&-0{,}05\\0{,}8&1\end{pmatrix}, det⁡=1+0,04=1,04\det=1+0{,}04=1{,}04: y1=11,04(10,05−0,81)(0,10)=(0,096154; −0,076923)\mathbf y_1=\frac1{1{,}04}\begin{pmatrix}1&0{,}05\\-0{,}8&1\end{pmatrix}\begin{pmatrix}0{,}1\\0\end{pmatrix}=(0{,}096154;\ -0{,}076923). Energia 0,002959+0,073964=0,076923=0,08/1,040{,}002959+0{,}073964=0{,}076923=0{,}08/1{,}04 ✓.
  • CN: (I−h2A)y1=(I+h2A)y0(I-\frac h2A)\mathbf y_1=(I+\frac h2A)\mathbf y_0 con h2=0,025\frac h2=0{,}025: (1−0,0250,41)y1=(0,1−0,04)\begin{pmatrix}1&-0{,}025\\0{,}4&1\end{pmatrix}\mathbf y_1=\begin{pmatrix}0{,}1\\-0{,}04\end{pmatrix}. Dalla seconda riga v1=−0,04−0,4 x1v_1=-0{,}04-0{,}4\,x_1; nella prima x1+0,001+0,01x1=0,1⇒x1=0,0991,01=0,098020x_1+0{,}001+0{,}01x_1=0{,}1\Rightarrow x_1=\frac{0{,}099}{1{,}01}=0{,}098020, v1=−0,079208v_1=-0{,}079208. Energia 12(0,079208)2+8⋅0,0980202=0,003137+0,076863=0,080000\frac12(0{,}079208)^2+8\cdot0{,}098020^2=0{,}003137+0{,}076863=0{,}080000 ✓ (conservata).

3. Energia dopo nn passi (En/E0E_n/E_0):

nn tt (s) EE (1,04n1{,}04^n) EI (1,04−n1{,}04^{-n}) CN
1 0,05 1,040 0,962 1
31 1,55 (≈T\approx T) 3,373 0,296 1
40 2,0 4,801 0,208 1
Dopo un periodo l'energia con EE è più che triplicata, e con EI è ridotta a meno di un terzo. Le posizioni: per n=31n=31 EE dà x=0,181x=0{,}181 m, maggiore dell'ampiezza iniziale 0,10{,}1 m (impossibile per un sistema conservativo), EI dà x=0,054x=0{,}054 m, CN dà x=0,0995x=0{,}0995 m contro l'esatto 0,1cos⁡(4⋅1,55)=0,09970{,}1\cos(4\cdot1{,}55)=0{,}0997 m (errore di fase piccolo).

4. Teoria. Con EE l'energia raddoppia quando 1,04n=2⇒n=ln⁡2ln⁡1,04=17,71{,}04^n=2\Rightarrow n=\frac{\ln2}{\ln1{,}04}=17{,}7 passi, cioè t=0,88t=0{,}88 s, circa 0,560{,}56 periodi. Con hh più piccolo l'energia cresce più lentamente: per raddoppiare in 10 periodi servono 1+ω2h2≈21/(10T/h)1+\omega^2h^2\approx2^{1/(10T/h)}, cioè ω2h2≈ln⁡2 h10 T\omega^2h^2\approx\frac{\ln2\,h}{10\,T}, quindi h≈ln⁡210⋅Tω2≈2,8⋅10−3h\approx\frac{\ln2}{10\cdot T\omega^2}\approx2{,}8\cdot10^{-3} s: più di 500 passi per periodo. Per tempi lunghi si usa CN (o altri metodi che conservano l'energia): errore solo di fase. La pulsazione numerica di CN è 2harctan⁡ωh2=3,9867\frac2h\arctan\frac{\omega h}2=3{,}9867 rad/s invece di 44 (periodo 1,5761{,}576 s invece di 1,5711{,}571 s): lo scarto relativo è (ωh)212=0,33%\frac{(\omega h)^2}{12}=0{,}33\%.

Risultati: fattori energia per passo 1,041{,}04 (EE), 0,96150{,}9615 (EI), 11 (CN); dopo un periodo E/E0=3,37; 0,296; 1E/E_0=3{,}37;\ 0{,}296;\ 1; EE raddoppia l'energia in ≈18\approx18 passi.

Errori tipici. Confondere il fattore sull'ampiezza (1,04\sqrt{1{,}04}) con quello sull'energia (1,041{,}04); applicare la regola h<2∣λ∣h<\frac2{|\lambda|} (per λ=±4i\lambda=\pm4i darebbe h<0,5h<0{,}5 ma EE è instabile per ogni hh); non risolvere il sistema lineare a ogni passo di EI e CN.

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

m=1m=1, k=16k=16, ω=4\omega=4, x0=0,1x_0=0{,}1, E0=0,08E_0=0{,}08 J, h=0,05h=0{,}05, A=(01−160)A=\begin{pmatrix}0&1\\-16&0\end{pmatrix}, z=±0,2iz=\pm0{,}2i.

  • Energia per passo: EE 1+ω2h2=1,041+\omega^2h^2=1{,}04; EI 11,04\frac1{1{,}04}; CN 11.
  • Primo passo: EE (0,1;−0,08)(0{,}1;-0{,}08), E=0,0832E=0{,}0832; EI 11,04(10,05−0,81)y0=(0,096154;−0,076923)\frac1{1{,}04}\begin{pmatrix}1&0{,}05\\-0{,}8&1\end{pmatrix}\mathbf y_0=(0{,}096154;-0{,}076923), E=0,076923E=0{,}076923; CN: x1=0,0991,01=0,098020x_1=\frac{0{,}099}{1{,}01}=0{,}098020, v1=−0,079208v_1=-0{,}079208, E=0,08000E=0{,}08000.
  • n=31n=31 (t=1,55≈Tt=1{,}55\approx T): E/E0=3,37; 0,296; 1E/E_0=3{,}37;\,0{,}296;\,1; x=0,181x=0{,}181 (EE, oltre l'ampiezza 0,10{,}1), 0,0540{,}054 (EI), 0,09950{,}0995 (CN; esatto 0,09970{,}0997). n=40n=40: 4,80; 0,208; 14{,}80;\,0{,}208;\,1.
  • EE raddoppia l'energia in ln⁡2ln⁡1,04=17,7\frac{\ln2}{\ln1{,}04}=17{,}7 passi (0,560{,}56 periodi). CN: pulsazione numerica 2harctan⁡ωh2=3,9867\frac2h\arctan\frac{\omega h}2=3{,}9867, errore relativo (ωh)212=0,33%\frac{(\omega h)^2}{12}=0{,}33\%.

Errori tipici: fattore sull'ampiezza confuso con quello sull'energia; h<2∣λ∣h<\frac2{|\lambda|} applicato a λ\lambda immaginario; sistema lineare non risolto.

Teoria collegata