\chapter{Studio delle Vibrazioni}
 
In questo capitolo ci occuperemo di capire come si comporta un sistema elastico, ovvero
un sistema di masse collegate da mollle. In questo modo potremo modellizzare, ad esempio,
strumenti musicali (come le corde di una chitarra o un tamburo) o un generale sistema
dinamico costituito da masse e molle.  %TODO: Inserire qualche altro esempio pratico
 
\section{Sistemi discreti}
\subsection{L'oscillatore armonico}
Il nostro scopo sarà, principalmente, risolvere sistemi complicati la cui soluzione esatta non
si riesce a calcolare (o non si riesce a farlo agevolmente) direttamente con un foglio ed una penna,
ma per capire di cosa ci occuperemo
è utile cominciare con un esempio noto. \\
\begin{wrapfigure}{l}{70mm}
\begin{tikzpicture}[scale=0.8]
 % Il terreno
 \draw[->] (0,0) -- (7,0) node[anchor=south west] {$x$};
 
 % La parete
 \draw (1,0) -- (1,2.8);
 \foreach \i in {0,...,12}
   \draw (1,0.2 * \i) -- (0.8, 0.2*\i + 0.2);
 
 % La molla
 \begin{scope}[thick]
  \foreach \i in {1,1.5,...,3.5}
    \draw (\i,0.5) -- (\i+0.25,1);
  \foreach \i in {1,1.5,...,3.5}
    \draw (\i+0.25,1) -- (\i+0.5,0.5);
 
 % E la massa
 \draw (4,0) rectangle (5.5,1.5);
 \fill[gray!28] (4,0) rectangle (5.5,1.5);
\end{scope}
 
\end{tikzpicture} \caption{L'oscillatore armonico} \label{fig:vib:oscarm}\end{wrapfigure}
Consideriamo il sistema costituito da una massa $m$ attaccata tramite una molla di costante
elastica $K$ e lunghezza a riposo nulla ad una parete (tralasciando la gravità), ovvero il sistema
in Figura~\ref{fig:vib:oscarm}.\\
Questo è un problema di cui conosciamo la soluzione esatta anche se introduciamo una dissipazione proporzionale
alla velocità $-\gamma v$. Abbiamo infatti che la forza applicata sul blocchetto è
\[
 F = -Kx - \gamma \dot{x}
\]
da cui , ricorando l'uguaglianza $F = ma = m\ddot{x}$ si ottiene l'equazione differenziale
\begin{equation}
 m \ddot{x} = -K x -\gamma \dot{x}
\end{equation}
Se aggiungiamo le condizioni iniziali $x(0) = x_0$ e $\dot{x}(0) = v_0$ la soluzione è anche unica.
Se consideriamo infatti una soluzione del tipo $x(t) = e^{\lambda t}$ abbiamo che
\[
 \dot{x}(t) = \lambda e^{\lambda t} \qquad \ddot{x}(t) = \lambda^2e^{\lambda t}
\]
e quindi $x$ è soluzione se e solo se
\[
 m \lambda^2 e^{\lambda t} = -K e^{\lambda t} - \gamma \lambda e^{\lambda t}
\]
e considerando che l'esponenziale non è mai nullo questo si verifica solo se
\[
 m \lambda^2 + \gamma \lambda + K = 0
\]
Escludendo il caso in cui il determinante è $0$ abbiamo due scelte per $\lambda$, che possono
essere reali o complesse. Nel caso reale i due esponenti ci danno una base dello spazio delle soluzioni.
Nel caso complesso possiamo osservare che, essendo le soluzioni coniugate, abbiamo che
\[
 \lambda = \alpha + \beta i \qquad \con{\lambda} = \alpha - \beta i
\]
e quindi
\[
 e^{\lambda t} = e^{\alpha}( \cos{\beta t} + i \sin{\beta t}) \qquad e^{\con{\lambda} t} =
 e^{\alpha}(\cos{\beta t} - i \sin{\beta t})
\]
dalle quali si ottiene che $e^{\alpha t}\cos{\beta t}$ e $e^{\alpha t}\sin{\beta t}$ sono ancora
una base per lo spazio delle soluzioni. \\
Nel caso con determinante nullo si può mostrare che prendendo $e^{\lambda t}$ e $te^{\lambda t}$ queste
sono ancora una base. In ogni caso, siamo in grado di risolvere esattamente il nostro problema imponendo
le condizioni iniziali.
 
Ci piacerebbe ora generalizzare questo procedimento in situazioni più generali e più complicate. In particolare,
vorremmo costruire un modello matematico di qualche situazione di reale interesse, come la seguente
 
\subsection{La corda di una chitarra}
Possiamo modellizzare la corda di una chitarra come in Figura~\ref{fig:vib:cordachitarra}, ovvero come $N$
masse collegate fra loro da $N-1$ molle e fissate, sempre tramite molle, a dei supporti laterali.
\begin{figure}[ht!]
 \begin{center}
  \begin{tikzpicture}
   % Supporto sx
   \draw (0,0) rectangle (1,2);
   \fill[gray!30] (0,0) rectangle (1,2);
 
   %% Serie di molle
   \draw[snake=snake] (1,1) -- (2,0.5);
   \draw (2.2,0.5) circle (0.2) node {$1$};
   \draw[snake=snake] (2.4,0.5) -- (3.2,1.1);
   \draw (3.4,1.1) circle (0.2) node {$2$};
   \draw[snake=snake] (3.6,1.1) -- (4.2,0.5);
   \draw (4.4,0.5) circle (0.2);
   \draw[snake=snake] (4.6,0.5) -- (5.8,1.5);
   \draw (6.0,1.5) circle (0.2);
   \draw[snake=snake] (6.2,1.5) -- (7.4,0.6);
   \draw (7.6,0.6) circle (0.2) node {$N$};
   \draw[snake=snake] (7.8,0.6) --(9,1);
 
   %% Bloccketto finale
   \draw (9,0) rectangle (10,2);
   \fill[gray!30] (9,0) rectangle (10,2);
  \end{tikzpicture}
  \caption{Modello di una corda di chitarra}
  \label{fig:vib:cordachitarra}
 \end{center}
\end{figure}
Questo modello non è più rappresentato da una equazione differenziale, ma da un sistema di equazioni.
In tutta generalità possiamo assumere che le particelle $p_i$ abbiano massa $m_i$ e che la molla che collega
$p_{i-1}$ con $p_{i}$ abbia costante elastica $k_i$. Indichiamo inoltre con $\gamma_i$ la costante di attrito
della $i$-esima particella. Possiamo scrivere l'equazione del moto di questa proiettandola sugli assi:
\begin{equation}
 \left\{ \begin{array}{ll}
 \ddot{x_i} m_i &= -k_i (x_i - x_{i-1}) - k_{i+1}(x_i - x_{i+1}) - \gamma_i \dot{x_i} \\
 \ddot{y_i} m_i &= -k_i (y_i - y_{i-1}) - k_{i+1}(y_i - y_{i+1}) - \gamma_i \dot{y_i}
 \end{array} \right.
\end{equation}
Queste possono essere risolte indipendemente per determinare il movimento orizzontale e verticale delle
particelle. Per affrontare la risoluzione è conveniente scrivere l'equazione in forma matriciale.
 
Sia $y = \trasp{(y_1, \ldots, y_n)}$, $M$ la matrice diagonale tale che $m_{ii} = m_i$ e $R$ la
matrice diagonale tale che $r_{ii} = \gamma_i$. Costruiamo infine la matrice $K$ costruita in questo
modo
\[
 K = \left[ \begin{array}{cccccc}
         \ddots & \ddots & & & & \\
          \ddots & \ddots & \ddots & & & \\
          & -k_i & k_i+k_{i+1} & k_{i+1} & \\
          & & \ddots & \ddots & \ddots & \\
        \end{array} \right]
\]
ovvero in modo che le equazioni di prima si possano scrivere nel seguente modo
\begin{equation}
 My'' + Ry' + Ky = 0
\end{equation}
Usiamo ora un piccolo trucco per riportare questa equazione differenziale del secondo ordine ad una
del primo. Definiamo $w = (y,y')$ il vettore ottenuto giustapponendo il vettore $y$ e la sua derivata.
Riscrivendo il sistema\footnote{Stiamo tacitamente assumendo che $M$ sia invertibile ma questo è piuttosto
ragionevole in quanto il determinante di questa è il prodotto delle masse e quindi è $0$ solo se almeno una
delle masse è $0$, situazione che non siamo interessati a rappresentare con il nostro modello} si ottiene
\begin{equation}
 w' = \left[ \begin{array}{cc|cc}
         \sblocko{0}{2} & \sblocke{I}{2} \\
          & & & \\ \hline
         \sblocko{M^{-1}K}{2} & \sblocke{M^{-1}R}{2} \\
          & & & \\
        \end{array} \right] w = Aw
\end{equation}
Supponiamo ora per semplicità che $A$ sia diagonalizzabile, ovvero $A = SDS^{-1}$. Si ha in questo caso
che $w = SDS^{-1}w$ e quindi $S^{-1}w = DS^{-1}w$ dove $D$ è diagonale. Ponendo $z = S^{-1}w$ si ottiene
$z_i = z_i d_{ii}$ e quindi ci si può ricondurre al caso precedente con la soluzione per componenti
$z_i(t) = e^{d_{ii} t}$.
\begin{os}
 Osserviamo che, in realtà, conosciamo la soluzione di questa equazione differenziale, e la conosciamo
 esplicitamente in ogni caso. Infatti è una verifica piuttosto immediata (e, più probabilmente, un fatto noto)
 che $\frac{d}{dt} \Exp{tA}w_0 = A\Exp(tA)w_0$ e quindi è soluzione del nostro problema di Cauchy. Sorge spontanea
la domanda: ``Perché abbiamo sviluppato tutto il ragionamento precedente?''. La riposta è che valutare questa
soluzione per un certo $t$ ha un costo computazionale molto alto. \`E quindi inaccettabile il costo che si ottiene
per la valutazione in tutti i punti che interessano a noi per ottenere un suono bem campionato (circa 44100 al secondo).
\end{os}
 
\section{Sistemi continui} \label{sec:vibsistcont}
\subsection{Generalizzazione del caso delle $N$ masse}
Supponiamo di avere una corda tesa fra due supporti fissi, che viene, in un'istante $t = 0$ messa in
una certa posizione con una data velocità iniziale. Siamo interessati a studiare il successivo moto
della corda e le sue vibrazioni, come abbiamo fatto nel caso di $N$ particelle discrete. \\
Consideriamo la seguente formalizzazione matematica di questo problema. Possiamo descrivere il moto della
corda parametrizzando i suoi punti tramite la coordinata orizzontale $x$ e introducendo una funzione
$u(x,t)$ che definisce la posizone del punto individuato dalla coordinata $x$ al tempo $t$, come si vede nella
Figura~\ref{fig:corda}. \\
\begin{wrapfigure}{l}{70mm}
 \begin{tikzpicture}
  % Blocchetto 1
  \draw (0,0) rectangle (1,2);
  \fill[gray!30] (0,0) rectangle (1,2);
 
  % Corda
  \draw[very thick] (1,1) .. controls (2.5,2) and (3.5,2) .. (5,1);
 
  % Secondo blocchetto
  \fill[gray!30] (5,0) rectangle (6,2);
  \draw (5,0) rectangle (6,2);
 
  % E qualche chiarimento :)
  \draw[->] (-0.5,-0.5) -- (6.5,-0.5) node[anchor=south] {$x$};
  \draw (1,-0.4) -- (1,-0.6) node[anchor=north] {$0$};
  \draw (5,-0.4) -- (5,-0.6) node[anchor=north] {$l$};
 
  \draw[dashed] (3,-0.5) -- (3,1.75) node[anchor=south] {$u(x,t)$};
 \end{tikzpicture}
 \caption{Rappresentazione della corda}
 \label{fig:corda}
\end{wrapfigure}
Un'analisi fisica del fenomeno\footnote{che non è nostro interesse mostrare ora} porta alla seguente
equazione differenziale:
\[
 \frac{\partial^2 u}{\partial t^2} - \sigma \frac{\partial^2 u}{\partial x^2} = 0
\]
dove $\sigma > 0$ è una costante determinata dalla gravità, dalla tensione e dalla densità
della corda. Prima di risolverla imponiamo le condizioni al contorno
\[
 u(0,t) = u(l,t)
\]
che ci permettono di controllare che gli estremi della corda siano sempre fissati ai supporti, e
le condizioni iniziali
\[
 u_0(x) = u(x,0)
\]
che dà la posizione iniziale e
\[
 v_0(x) = \frac{\partial u}{\partial t}(x,0)
\]
che esprime la velocità iniziale di ogni particella. Questa formalizzazione si può facilmente estendere
supponendo di avere una membrana introducendo $u(x,y,t)$.
 
Per ora siamo particolarmente interessati a studiare se esistono delle soluzioni in cui tutti i punti
si muovano con la stessa frequenza $\omega$, anche se eventulamente con diversa ampiezza. Si tratta
di trovare una soluzione della forma
\[
 u(x,t) = v(x) \cdot \cos{(\omega t)}
\]
Osserviamo che data una $u$ di quest forma le derivate parziali diventano
\[
 \frac{\partial^2 u}{\partial t^2} = -v(x) \omega^2 \cos(wt) \qquad \frac{\partial^2 u}{\partial x^2} = v''(x) \cos(wt)
\]
e quindi l'equazione differenziale risulta essere
\[
  \cos(wt) \left[ -\omega^2 v(x) - \sigma v''(x) \right] = 0
\]
che può essere identicamente nulla solo se
\[
 \omega^2 v(x) + \sigma v''(x) \equiv 0 \iff v''(x) = -\frac{\omega^2 v(x)}{\sigma}
\]
 
\subsection{Discretizzazione del modello}
 
Per risolvere effettivamente il problema supponiamo di suddividere l'intervallo $[0,l]$ in $n+1$
punti equidistanti $x_0, \ldots, x_{n+1}$ tali che $x_0 = 0$ e $x_{n+1} = l$. Si avrà dunque
che $x_i - x_{i-1} = h = \frac{l}{n+1}$. \\
Per comodità poniamo $v(x_i) = v_i$ la posizione dell'$i$-esimo
punto.
Assumiamo ora che la soluzione dell'equazione differenziale sia di classe $C^4$ su $[0,l]$, e facciamo
lo sviluppo in serie centrato in un generico punto $x_i$.
\[
 v(x_i + h) = v_{i+1} =  v_i + hv'(x_i) + \frac{h^2 v''(x_i)}{2} + \frac{h^3 v^{(3)}(x_i)}{6} + \frac{h^4 v^{(4)}(\xi_i)}{24}
\]
dove $\xi_i \in (x_i , x_{i+1})$.
Analogamente si ottiene
\[
 v(x_i - h) = v_{i-1} = v_i - hv'(x_i) + \frac{h^2 v''(x_i)}{2} - \frac{h^3 v^{(3)}(x_i)}{6} + \frac{h^4 v^{(4)}(\eta_i)}{24}
\]
dove $\eta_i \in (x_{i-1}, x_i)$. Sommando queste due si ha
\[
 v_{i+1} + v_{i-1} = 2v_i + h^2 v''(x_i) + \frac{h^4}{24}[v^{(4)}(\xi_i) + v^{(4)}(\eta_i)]
\]
e possiamo ancora una volta ottenere un'espressione ``quasi'' esplicita della $v''(x_i)$:
\[
 v''(x_i) = \frac{1}{h^2} ( v_{i+1} + v_{i-1} - 2v_i ) - \frac{h^2}{24}[v^{(4)}(\xi_i) + v^{(4)}(\eta_i)]
\]
Nel nostro modello $h$ è destinato ad essere piccolo, perché $h = \frac{1}{n+1}$ e intuitivamente $n$ deve essere un
numero grande se vogliamo approssimare bene il problema continuo. Abbiamo quindi che $h^2 \to 0$ quando $n \to \infty$
e vorremmo poter dire che possiamo trascurare il termine con le derivate quarte.
Se così fosse avremmo infatti ricondotto il nostro problema ad un problema di autovalori.
 
Possiamo riscrivere infatti le relazione in questo modo
\[
 Av - \frac{h^2}{24} \tau = \lambda v
\]
dove $v = \trasp{( v_1 \ldots v_n )}$ e
\[
 A = \left[ \begin{array}{cccccc}
             -2 & 1 &  & & & \\
              1 & -2 & 1 & & & \\
              & \ddots & \ddots & \ddots & & \\
              & & 1 & -2 & 1&   \\
              & & & 1 & -2 \\
            \end{array} \right]
\]
e dove $\tau$ è un opportuno vettore che contiene i termini del tipo $\tau_i = v^{(4)}(\xi_i) + v^{(4)}(\eta_i)$.
Osserviamo che possiamo dare una maggiorazione alla norma di $\tau$ perché essendo la derivata quarta di $v$ continua
su un compatto è limitata, e quindi $\forall i \: |\tau_i | \leq M$ per un opportuno $M$.
 
Ricordiamo il Teorema~\ref{te:BauerFike}; questo dice che presa una norma assoluta ed $A$ diagonalizzabile ed $F$
una matrice di perturbazione per ogni autovalore $\lambda$ di $A + F$ esiste un autovalore $\mu$ di $A$
tale che $||\lambda - \mu|| \leq \cond{T}||F||$ dove $T$ è il cambio di base che diagonalizza $A$.
Se poniamo $F = -\frac{h^2}{24}\tau \trasp v$ dove $v$ è un autovettore di $v$ è tale che
$Av - \frac{h^2}{24}\tau = \lambda v$ allora si ha
\[
 (A - \frac{h^2}{24} \tau \trasp v) v = \lambda v
\]
e quindi esiste $\mu$ autovalore di $A$ tale che\footnote{Osserviamo che, essendo $A$ reale simmetrica
il cambio di base che la diagonalizza è ortogonale e quindi ha condizionamento $1$.} $||\lambda - \mu|| \leq ||F||
= \frac{h^2}{24} \cdot \frac{||\tau\trasp v||}{\trasp vv} \leq \frac{h^2}{24\trasp vv} \sqrt{\rs{\trasp \tau v\trasp v\tau}}
= ||\tau|| \frac{h^2}{24 \trasp vv}$; questo si ottiene osservando che la norma di una matrice diade $v\trasp v$
è data dalla norma du $v$.
 
Vogliamo ora mostrare che $\frac{||\tau||}{||v||}$ è limitato superiormente da una costante, in modo
da provare la convergenza della nostra approssimazione.
 
Consideriamo la seguente serie di eguaglianze
\[
 \sqrt{h} || \tau ||_2 = ( \sum_{i=1}^{n} v^{(4)}(\xi_i)^2 + v{(4)}(\eta_i)^2)^{\frac 1 2} \sqrt{h}
 = ( \sum_{i=1}^{n} h v^{(4)}(\xi_i)^2 + h v{(4)}(\eta_i)^2)^{\frac 1 2}
\]
Osserviamo che l'ultima espressione può essere intesa, per $h \to 0$, come un'approssimazione dell'integrale
di Riemann di $v^2(x)$, ed essendo l'integrale di una funzione continua su un compatto è limitato
superiormente da una costante. Abbiamo quindi provato la tesi.
 
Introduciamo ora un metodo iterativo che possiamo vedere come una estensione del metodo delle potenze
e che ci permetterà di studiare nel dettaglio questa situazione.
 
\subsection{Il metodo delle iterazioni ortogonali di sottospazi}
Il metodo delle potenze come è stato analizzato nella Sezione~\ref{sec:metodopotenze} è basato sull'idea
di prendere un vettore e iterare su di lui l'applicazione individuata dalla matrice $A$ di cui si vogliono
calcolare gli autovalori. \\
Osserviamo che questa idea è equivalente (anzi, forse è meglio rappresentata) all'iterare un sottospazio
vettoriale di dimensione $1$. Supponiamo ora di voler calcolare, invece che 1, più autovalori delle matrice,
precisamente $k$ (dove $k < n$ ed $n$ è la dimensione delle matrice) ognuno di molteplicità 1\footnote{%
Questo solo per semplicità di esposizione. Non è complicato generalizzare l'idea}. Una naturale estensione del metodo
sarebbe iterare un sottospazio vettoriale di dimensione $k$. Questo dovrebbe convergere ad un sottospazio di dimensione
$k$ uguale alla somma degli autospazi dei primi $k$ autovalori.\\
Si può osservare subito un problema computazionale. Sebbene questo discorso sia teoricamente corretto, è facile prevedere
che il nostro spazio di dimensione $k$ convergerà verso l'autospazio relativo
all'autovalore dominante (che è di dimensione 1). Più precisamente, supponiamo di rappresentare lo spazio
con una matrice $x$ $n \times k$ le cui colonne siano una base dello spazio. Allora l'iterazione sarebbe
\[
\left\{\begin{array}{lcl}
 x_0 &=& x \\
 x_{k+1} &=& Ax_k
\end{array} \right.
\]
Per lo stesso ragionamento fatto nel caso del metodo delle potenze originale tutte le colonne di $A$ convergeranno
a multipli dello stesso vettore (quello dell'autovalore dominante).
 
Un modo di risolvere questo problema potrebbe essere assicurarsi che durante l'iterazione tutte le colonne restino
ortogonali fra loro. Consideriamo la seguente iterazione
\[
 \left\{ \begin{array}{lcl}
          x_0 &=& x \\
          T_{k} &=& A x_k \\
          Q_k R_k &=& T_k \: \text{(scomposizione $QR$ di $T_k$)} \\
          x_{k+1} &=& Q_k
         \end{array} \right.
\]
Ricordando che la scomposizione QR di una matrice $n \times k$ con $k < n$ è tale che $Q$ è rettangolare
ed è $n \times k$ e $R$ è quadrata $k \times k$. In particolare si osserva che nella nostra iterazione
le colonne di $Q_k$ sono una
base ortonormale per lo spazio generato dalle colonne di $T_k$.
Osserviamo che calcolare esplicitamente $x_k$ non è realmente necessario e possiamo riformulare
l'iterazione. Poniamo
\[
 B_k = \trasp{Q_k} A Q_k \qquad \text{e} \qquad B_k = U_k D \trasp{U_k}
\]
dove $U_K D \trasp{U_k}$ è la decomposizione spettrale di $B_k$. Si osserva che $T_{k+1} = A Q_k =
Q_k B_k = $.
% TODO: Capire l'iterazione implicita.