Aula 1: Do Oscilador Harmônico às Matrizes Circulantes

Minicurso: Osciladores Acoplados e a Equação de Klein–Gordon

Sandro Vitenti

Roteiro do minicurso

Aula 1

1. O oscilador harmônico: espaço de fase e integração numérica
2. A cadeia de osciladores
3. Matrizes circulantes e a transformada de Fourier discreta

Aula 2

4. O limite contínuo e a equação de Klein–Gordon
5. Variáveis de criação e aniquilação clássicas

Notas, código e notebooks: página do minicurso

O oscilador harmônico

Equações de Hamilton e espaço de fase

m\ddot x = -\kappa x, \qquad \omega \equiv \sqrt{\kappa/m}, \qquad H(x, p) \equiv \frac{p^2}{2m} + \frac{\kappa x^2}{2}, \qquad \dot x = \frac{\partial H}{\partial p} = \frac{p}{m}, \qquad \dot p = -\frac{\partial H}{\partial x} = -\kappa x

Once Loop Reflect

\dot x \equiv \mathrm dx/\mathrm dt. Notas: Equações de Hamilton

O gradiente e a matriz \Omega

\vec\chi \equiv \begin{pmatrix} x \\ p \end{pmatrix}, \qquad \vec\partial H \equiv \begin{pmatrix} \partial H/\partial x \\ \partial H/\partial p \end{pmatrix} = \begin{pmatrix} \kappa x \\ p/m \end{pmatrix} \dot{\vec\chi} = \begin{pmatrix} \partial H/\partial p \\ -\partial H/\partial x \end{pmatrix} = \Omega\,\vec\partial H, \qquad \Omega \equiv \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}

• Hamilton troca as componentes do gradiente: \Omega gira de 90^\circ, a velocidade é tangente a H = \text{constante}.
• \dfrac{\mathrm dH}{\mathrm dt} = \left(\vec\partial H\right)^{\mathsf T}\dot{\vec\chi} = \left(\vec\partial H\right)^{\mathsf T}\Omega\,\vec\partial H = 0.

Figura 1: Curvas H = 1/4, 1 e 9/4 (cinza) do oscilador com m = 1 e \kappa = 2. Em oito pontos da curva H = 1: o gradiente \vec\partial H (azul), perpendicular à curva, e a velocidade \dot{\vec\chi} = \Omega\,\vec\partial H (vermelho), o gradiente girado de 90^\circ no sentido horário, tangente à curva. Os dois eixos estão na mesma escala, para que os ângulos retos apareçam como tais.

Notas: Forma matricial

Forma matricial

H = \frac{1}{2} \vec\chi^{\mathsf T}\, \mathsf H\, \vec\chi, \qquad \mathsf H \equiv \begin{pmatrix} \kappa & 0 \\ 0 & 1/m \end{pmatrix}, \qquad \vec\partial H = \mathsf H\,\vec\chi \dot{\vec\chi} = \Omega\, \mathsf H\, \vec\chi \equiv A \vec\chi, \qquad A = \begin{pmatrix} 0 & 1/m \\ -\kappa & 0 \end{pmatrix}

• \Omega é a mesma para todo sistema; o sistema entra apenas por \mathsf H.
• Sobre as soluções, \mathrm d/\mathrm dt é representada por A: \ddot{\vec\chi} = A\dot{\vec\chi} = A^2\vec\chi, \qquad \frac{\mathrm d^n\vec\chi}{\mathrm dt^n} = A^n\vec\chi

Notas: Forma matricial

A exponencial de uma matriz

Série de Taylor, trocando cada derivada por uma potência de A: \vec\chi(t) = \sum_{n=0}^\infty \frac{t^n}{n!}\,\frac{\mathrm d^n\vec\chi}{\mathrm dt^n}(0) = \left(I + tA + \frac{t^2A^2}{2} + \frac{t^3A^3}{6} + \dots\right)\vec\chi(0) \equiv e^{tA}\,\vec\chi(0)

• Com um número a no lugar de A: \dot x = ax, x(t) = e^{at}x(0).
• Derivando a série termo a termo, \dfrac{\mathrm d}{\mathrm dt}e^{tA} = A\,e^{tA} e e^{0A} = I.

Notas: Solução exata

Potências de A

A^2 = \begin{pmatrix} 0 & 1/m \\ -\kappa & 0 \end{pmatrix} \begin{pmatrix} 0 & 1/m \\ -\kappa & 0 \end{pmatrix} = \begin{pmatrix} -\kappa/m & 0 \\ 0 & -\kappa/m \end{pmatrix} = -\omega^2 I A^3 = -\omega^2 A, \qquad A^4 = \omega^4 I, \qquad A^5 = \omega^4 A, \qquad A^6 = -\omega^6 I, \qquad\dots A^{2n} = (-\omega^2)^n I, \qquad A^{2n+1} = (-\omega^2)^n A Separando pares e ímpares: e^{tA} = \underbrace{\left(1 - \frac{\omega^2t^2}{2} + \frac{\omega^4t^4}{24} - \dots\right)}_{\cos\omega t} I + \underbrace{\left(t - \frac{\omega^2t^3}{6} + \frac{\omega^4t^5}{120} - \dots\right)}_{\sin(\omega t)/\omega} A

Notas: Potências de A

Solução exata

e^{tA} = \cos(\omega t)\, I + \frac{\sin(\omega t)}{\omega}\, A = \begin{pmatrix} \cos\omega t & \dfrac{\sin\omega t}{m\omega} \\[2ex] -m\omega\sin\omega t & \cos\omega t \end{pmatrix} x(t) = x(0)\cos\omega t + \frac{p(0)}{m\omega}\sin\omega t, \qquad p(t) = p(0)\cos\omega t - m\omega\, x(0)\sin\omega t

• As trajetórias são as elipses H = \text{constante}.

Notas: Soma da série

O produto simplético é uma área

\Omega(\vec\chi_1, \vec\chi_2) \equiv \vec\chi_1^{\mathsf T}\Omega\,\vec\chi_2 = x_1p_2 - p_1x_2 = \det\begin{pmatrix} x_1 & x_2 \\ p_1 & p_2 \end{pmatrix}

Figura 2: O paralelogramo de lados \vec\chi_1(t) (azul) e \vec\chi_2(t) (laranja) em três instantes, para o oscilador com m = 1, \kappa = 4 (\omega = 2) e \vec\chi_1(0) = (1;\ 0{,}4), \vec\chi_2(0) = (0{,}2;\ 1{,}2). Cada vetor percorre a sua elipse H = \text{constante} (cinza). A forma do paralelogramo muda, e a sua área, \Omega(\vec\chi_1, \vec\chi_2) = 1{,}1200, é a mesma nos três painéis.

• \frac{\mathrm d}{\mathrm dt}\Omega(\vec\chi_1, \vec\chi_2) = \vec\chi_1^{\mathsf T}\left(A^{\mathsf T}\Omega + \Omega A\right)\vec\chi_2 = 0, pois A^{\mathsf T}\Omega + \Omega A = 0 para qualquer \mathsf H simétrica (Liouville).

Notas: O produto simplético

Integração numérica: um passo exato é uma rotação

\xi \equiv \sqrt{m}\,\omega\, x, \qquad \eta \equiv \frac{p}{\sqrt m}, \qquad \varepsilon \equiv \omega\,\Delta t H = \frac{\xi^2 + \eta^2}{2}, \qquad \frac{\mathrm d}{\mathrm dt}\begin{pmatrix} \xi \\ \eta \end{pmatrix} = \omega\,\Omega\begin{pmatrix} \xi \\ \eta \end{pmatrix} M_{\rm ex}(\varepsilon) \equiv \cos\varepsilon\, I + \sin\varepsilon\, \Omega = \begin{pmatrix} \cos\varepsilon & \sin\varepsilon \\ -\sin\varepsilon & \cos\varepsilon \end{pmatrix}

Em z \equiv \xi + i\eta: z_{n+1} = e^{-i\varepsilon} z_n.

Figura 3: Um passo exato no plano (\xi, \eta): o ponto z_n vai para z_{n+1} = e^{-i\varepsilon} z_n, uma rotação de ângulo \varepsilon no sentido horário sobre a circunferência H = |z|^2/2 = \text{constante}. Na figura, \varepsilon = 0{,}6.

Notas: Variáveis adimensionais

Os métodos como aproximações da rotação

Cada método substitui M_{\rm ex}(\varepsilon) por uma matriz M(\varepsilon) polinomial em \varepsilon: \vec\chi_{n+1} = M(\varepsilon)\,\vec\chi_n.

1. Ordem. M - M_{\rm ex} = O(\varepsilon^{k+1}): erro acumulado O(\Delta t^k) depois de n = t/\Delta t passos.
2. Produto simplético. Queremos M^{\mathsf T}\Omega M = \Omega, como a evolução exata.
3. Autovalores. \vec\chi_n = M^n\vec\chi_0 cresce como |\lambda|^n. Exato: \lambda = e^{\pm i\varepsilon}, |\lambda| = 1.

Para qualquer matriz 2\times 2: M^{\mathsf T}\Omega M = \det(M)\,\Omega, \qquad \Omega(M\vec\chi_1, M\vec\chi_2) = \det(M)\,\Omega(\vec\chi_1, \vec\chi_2)

• \det M mede a mudança das áreas; em duas dimensões, simplética \Leftrightarrow \det M = 1.

Notas: Os métodos como aproximações da rotação

Autovalores: traço e determinante

\det(M - \lambda I) = \lambda^2 - (\operatorname{tr} M)\,\lambda + \det M, \qquad \lambda_1 + \lambda_2 = \operatorname{tr} M, \qquad \lambda_1\lambda_2 = \det M

• (\operatorname{tr} M)^2 < 4\det M: par complexo conjugado, |\lambda|^2 = \det M. Cresce se \det M > 1.
• \det M = 1, |\operatorname{tr} M| < 2: \lambda = e^{\pm i\theta}, 2\cos\theta = \operatorname{tr} M. Limitada; \theta \neq \varepsilon é um erro de fase.
• \det M = 1, |\operatorname{tr} M| > 2: autovalores reais, um deles com |\lambda| > 1. Cresce exponencialmente.
• Exato: \operatorname{tr} M_{\rm ex} = 2\cos\varepsilon, \theta = \varepsilon.

Notas: Os métodos como aproximações da rotação

Euler explícito

Taylor truncado: \vec\chi(t_{n+1}) = \vec\chi(t_n) + \Delta t\,A\vec\chi(t_n) + O(\Delta t^2) \vec\chi_{n+1} = (I + \Delta t\, A)\, \vec\chi_n, \qquad M_{\rm E}(\varepsilon) = I + \varepsilon\Omega = \begin{pmatrix} 1 & \varepsilon \\ -\varepsilon & 1 \end{pmatrix} = \sqrt{1+\varepsilon^2}\;M_{\rm ex}(\arctan\varepsilon)

• \operatorname{tr} M_{\rm E} = 2, \det M_{\rm E} = 1 + \varepsilon^2: \lambda = 1 \pm i\varepsilon, |\lambda| = \sqrt{1+\varepsilon^2}. Primeira ordem.
• H_n = (1+\varepsilon^2)^n H_0: espiral que se afasta da origem.
• Até o instante t: n = \omega t/\varepsilon passos e H_n/H_0 \approx e^{\varepsilon\,\omega t}. Passo menor: menos erro, mais passos.
• 20 períodos com erro de 1\%: \varepsilon < 7{,}9\times 10^{-5}, 1{,}6\times 10^6 passos.

Notas: Euler explícito

Euler simplético: dois cisalhamentos

Figura 4: Um passo do Euler simplético com \varepsilon = 0{,}8 aplicado ao disco H \le 1/2 (sombreado) e ao ponto (\xi, \eta) = (1, 0) (vermelho). Esquerda: estado inicial. Centro: após S_\eta, cada ponto desce em \varepsilon\xi. Direita: após S_\xi, cada ponto se desloca horizontalmente em \varepsilon\eta. A área sombreada, calculada pela fórmula do cadarço, é 3{,}1415 nos três painéis. O ponto vazio no painel da direita é a rotação exata de (1, 0).

Notas: Euler simplético

Euler simplético

S_\eta(\varepsilon) \equiv \begin{pmatrix} 1 & 0 \\ -\varepsilon & 1 \end{pmatrix}, \qquad S_\xi(\varepsilon) \equiv \begin{pmatrix} 1 & \varepsilon \\ 0 & 1 \end{pmatrix}, \qquad \det S_\eta = \det S_\xi = 1 M_{\rm SE}(\varepsilon) = S_\xi(\varepsilon)\, S_\eta(\varepsilon) = \begin{pmatrix} 1 - \varepsilon^2 & \varepsilon \\ -\varepsilon & 1 \end{pmatrix}, \qquad \det M_{\rm SE} = 1, \qquad \operatorname{tr} M_{\rm SE} = 2 - \varepsilon^2

• S_\xi e S_\eta geram o grupo simplético \mathrm{Sp}(2,\mathbb R) = \{M : M^{\mathsf T}\Omega M = \Omega\}: toda M com M_{12} \neq 0 é S_\eta(\cdot)\,S_\xi(M_{12})\,S_\eta(\cdot).
• Primeira ordem, mas o erro deforma em vez de dilatar.
• Conserva exatamente \tilde H_{\rm SE} \equiv (\xi^2 + \eta^2 - \varepsilon\xi\eta)/2: a trajetória é uma elipse inclinada.
• H oscila com amplitude relativa \approx\varepsilon/2, sem crescer.

Notas: Euler simplético

Verlet: três cisalhamentos

Meio passo em p, passo em x, meio passo em p: M_{\rm V}(\varepsilon) = S_\eta(\varepsilon/2)\, S_\xi(\varepsilon)\, S_\eta(\varepsilon/2) = \begin{pmatrix} 1 - \varepsilon^2/2 & \varepsilon \\ -\varepsilon + \varepsilon^3/4 & 1 - \varepsilon^2/2 \end{pmatrix} M_{\rm V} - M_{\rm ex} = \begin{pmatrix} -\varepsilon^4/24 & \varepsilon^3/6 \\ \varepsilon^3/12 & -\varepsilon^4/24 \end{pmatrix} + O(\varepsilon^5)

• Segunda ordem: M_{\rm V}(-\varepsilon) = M_{\rm V}(\varepsilon)^{-1} elimina o erro de ordem \varepsilon^2.
• \operatorname{tr} M_{\rm V} = 2 - \varepsilon^2 = 2\cos\theta_{\rm V}: avança \theta_{\rm V} = 2\arcsin(\varepsilon/2) = \varepsilon + \varepsilon^3/24 + \dots por passo.
• Estável somente se \omega\Delta t < 2. Conserva \tilde H_{\rm V} \equiv \eta^2/2 + (1 - \varepsilon^2/4)\,\xi^2/2.

Notas: Verlet

Runge–Kutta de quarta ordem

Quatro velocidades, média com os pesos 1, 2, 2, 1 de Simpson: k_1 = f(\vec\chi_n), \qquad k_2 = f\!\left(\vec\chi_n + \tfrac{\Delta t}{2}k_1\right), \qquad k_3 = f\!\left(\vec\chi_n + \tfrac{\Delta t}{2}k_2\right), \qquad k_4 = f\!\left(\vec\chi_n + \Delta t\,k_3\right) \vec\chi_{n+1} = \vec\chi_n + \frac{\Delta t}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right), \qquad M_{\rm RK4} = \left(1 - \frac{\varepsilon^2}{2} + \frac{\varepsilon^4}{24}\right) I + \left(\varepsilon - \frac{\varepsilon^3}{6}\right)\Omega

• Para f(\vec\chi) = A\vec\chi a média dá I + \Delta t\,A + \dots + \Delta t^4A^4/24, os cinco primeiros termos da série de e^{\Delta t\,A}.
• M_{\rm RK4} é a rotação com \cos\varepsilon e \sin\varepsilon truncados na ordem \varepsilon^4: quarta ordem.
• \dfrac{H_{n+1}}{H_n} = \det M_{\rm RK4} = 1 - \dfrac{\varepsilon^6}{72} + \dfrac{\varepsilon^8}{576}: erro pequeno, sempre com o mesmo sinal.

Notas: RK4

Resumo

Método	M(\varepsilon)	Ordem	\det M	\lvert\lambda\rvert
Exato	\cos\varepsilon\, I + \sin\varepsilon\, \Omega		1	1
Euler explícito	I + \varepsilon\Omega	1	1 + \varepsilon^2	\sqrt{1 + \varepsilon^2}
Euler simplético	S_\xi(\varepsilon)\, S_\eta(\varepsilon)	1	1	1 para \varepsilon < 2
Verlet	S_\eta(\varepsilon/2)\, S_\xi(\varepsilon)\, S_\eta(\varepsilon/2)	2	1	1 para \varepsilon < 2
RK4	\sum_{k=0}^{4} (\varepsilon\Omega)^k/k!	4	1 - \varepsilon^6/72 + \varepsilon^8/576	\sqrt{\det M}

• Simpléticos: \det M = 1, |\lambda| = 1, a energia oscila e não deriva.
• Não simpléticos: |\lambda| \neq 1, a energia muda com o mesmo sinal a cada passo.

Notas: Resumo

Oito passos com \varepsilon = 0{,}8

Figura 5: Oito passos de cada método com \varepsilon = 0{,}8 a partir de (\xi, \eta) = (1, 0) (pontos cheios), comparados com os oito passos exatos sobre a circunferência H = 1/2 (pontos vazios). Euler explícito sai da circunferência: H_8/H_0 = (1 + \varepsilon^2)^8 = 52{,}3. Euler simplético e Verlet permanecem sobre as elipses tracejadas \tilde H_{\rm SE} = 1/2 e \tilde H_{\rm V} = 1/2, e H/H_0 varia entre 0{,}77 e 1{,}57 (Euler simplético) e entre 0{,}84 e 1{,}00 (Verlet). RK4: H_8/H_0 = 0{,}974.

Notas: Implementação e comparação

Implementação

def passo_rk4(f, chi, dt):
    k1 = f(chi)
    k2 = f(chi + 0.5 * dt * k1)
    k3 = f(chi + 0.5 * dt * k2)
    k4 = f(chi + dt * k3)
    return chi + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6

def passo_verlet(F, x, p, dt, m):
    p = p + 0.5 * dt * F(x)
    x = x + dt * p / m
    p = p + 0.5 * dt * F(x)
    return x, p

As mesmas funções servem para a cadeia: x, p e F passam a ser vetores.

Notas: Implementação e comparação

Vinte períodos: erro na energia

Figura 6: Erro relativo na energia do oscilador harmônico após 2513 passos com \varepsilon = 0{,}05. Euler explícito: erro final 5{,}3\times 10^{2}, igual a (1+\varepsilon^2)^n - 1. RK4: 5{,}5\times 10^{-7}, igual a n\varepsilon^6/72. Euler simplético e Verlet: erro oscilante, com máximos 2{,}56\times 10^{-2} (\approx\varepsilon/2) e 6{,}25\times 10^{-4} (=\varepsilon^2/4), que não crescem com o tempo.

Vinte períodos: espaço de fase

Figura 7: Trajetórias no espaço de fase durante 20 períodos. Euler explícito: espiral que se afasta da origem (a área não é preservada). Verlet: curva fechada, a elipse \tilde H_{\rm V} = \text{constante}.

A cadeia de osciladores

O modelo

Figura 8: Esquerda: um trecho da cadeia, com as molas de constante K entre vizinhas (azul), as molas de constante K_0 que prendem cada partícula à sua posição de equilíbrio (verde, a linha pontilhada é h = 0) e o espaçamento horizontal \Delta x. Direita: as N = 8 partículas vistas como um anel; a mola entre as partículas N - 1 = 7 e 0 (laranja) é a condição de contorno periódica.

• N partículas de massa m, alturas h^j(t), molas de comprimento natural nulo, h^N \equiv h^0.

Notas: O modelo

Forças: as componentes horizontais se cancelam

Figura 9: Forças sobre a partícula j para K = 1, K_0 = 0{,}5, \Delta x = 1 e alturas h^{j-1} = 0{,}3, h^j = 0{,}8, h^{j+1} = 0{,}5. As molas de constante K exercem K(\Delta x,\ h^{j+1} - h^j) = (1;\ -0{,}3) e K(-\Delta x,\ h^{j-1} - h^j) = (-1;\ -0{,}5) (azul); as componentes horizontais (pontilhadas) se cancelam e a soma é (0;\ -0{,}8) = (0;\ K(h^{j+1} - 2h^j + h^{j-1})) (vermelho). A mola de constante K_0 exerce (0;\ -K_0h^j) = (0;\ -0{,}4) (verde).

• Molas de comprimento natural nulo: força = K \times (vetor até a outra extremidade).
• Horizontais: K\Delta x - K\Delta x = 0. Partindo do repouso, as posições horizontais não mudam.

Notas: Forças e equações de movimento

Equações de movimento

K\left(h^{j+1} - h^j\right) + K\left(h^{j-1} - h^j\right) = K\Delta x\left(\frac{h^{j+1} - h^j}{\Delta x} - \frac{h^j - h^{j-1}}{\Delta x}\right) = K\left(h^{j+1} - 2h^j + h^{j-1}\right) m\ddot h^j = F_j = K\left(h^{j+1} - 2h^j + h^{j-1}\right) - K_0\, h^j

• Diferença de inclinações: a força é proporcional à curvatura (segunda derivada discreta).
• K = 0: N osciladores independentes de frequência \sqrt{K_0/m}.
• K \neq 0: a equação de h^j contém h^{j \pm 1}, os osciladores estão acoplados.
• Exata para qualquer inclinação; com comprimento natural \ell_0, só na ordem linear.

Notas: Forças e equações de movimento

O operador de translação

(T\vec h)^j \equiv h^{j+1}, \qquad T = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \end{pmatrix}, \qquad T\begin{pmatrix} h^0 \\ h^1 \\ h^2 \\ h^3 \end{pmatrix} = \begin{pmatrix} h^1 \\ h^2 \\ h^3 \\ h^0 \end{pmatrix} \qquad (N = 4)

• Matriz de permutação: T^{-1} = T^{\mathsf T}, (T^{-1}\vec h)^j = h^{j-1}, e T^N = I.
• O 1 no canto inferior esquerdo é a condição periódica.
• No numpy: np.roll(h, -1) é T\vec h e np.roll(h, 1) é T^{-1}\vec h.

Notas: O operador de translação

A matriz D

\vec F = K\left(T\vec h - \vec h\right) + K\left(T^{-1}\vec h - \vec h\right) - K_0\,\vec h = -\left(KD + K_0 I\right)\vec h, \qquad D \equiv 2I - T - T^{-1}, \qquad m\ddot{\vec h} = \vec F

Figura 10: As matrizes T, T^{-1} e D para N = 8, com os elementos não nulos escritos nas casas. Os elementos nos cantos superior direito e inferior esquerdo vêm da mola entre as partículas N - 1 e 0.

Notas: O operador de translação

Hamiltoniana

V = \frac{K}{2}\sum_j\left(h^{j+1} - h^j\right)^2 + \frac{K_0}{2}\sum_j \left(h^j\right)^2 = \frac12\, \vec h^{\mathsf T}\left(KD + K_0 I\right) \vec h, \qquad H = \sum_j \frac{p_j^2}{2m} + V \vec\chi \equiv \begin{pmatrix} \vec h \\ \vec p \end{pmatrix} = \begin{pmatrix} h^0 \\ \vdots \\ h^{N-1} \\ p_0 \\ \vdots \\ p_{N-1} \end{pmatrix}, \qquad \vec\partial H = \begin{pmatrix} \left(KD + K_0 I\right)\vec h \\ \vec p/m \end{pmatrix}, \qquad \Omega = \begin{pmatrix} 0 & I \\ -I & 0 \end{pmatrix}, \qquad A = \Omega\,\mathsf H = \begin{pmatrix} 0 & I/m \\ -\left(KD + K_0 I\right) & 0 \end{pmatrix}

• \dot{\vec\chi} = \Omega\,\vec\partial H = A\vec\chi: o oscilador isolado com \kappa \to KD + K_0 I, blocos N \times N.
• A^2 = -\frac1m\operatorname{diag}\left(KD + K_0 I,\ KD + K_0 I\right) não é múltiplo de I: é preciso diagonalizar KD + K_0 I.

Notas: Energia e Hamiltoniana

O modo zero

\vec 1 \equiv \begin{pmatrix} 1 & 1 & \cdots & 1 \end{pmatrix}^{\mathsf T}, \qquad K_0 = 0:\ H \text{ invariante por } \vec h \to \vec h + \beta\,\vec 1, \qquad D\,\vec 1 = 0 \bar h \equiv \frac1N\sum_j h^j = \frac{\vec 1^{\mathsf T}\vec h}{N}, \qquad \bar p \equiv \frac1N\sum_j p_j = \frac PN, \qquad \dot{\bar h} = \frac{\bar p}{m}, \qquad \dot{\bar p} = -K_0\,\bar h

• O par (\bar h, \bar p) é um oscilador com \kappa = K_0, desacoplado do resto da cadeia.
• K_0 = 0: \bar p conservado, \bar h(t) = \bar h(0) + \bar p\,t/m. Partícula livre, frequência zero.
• K_0 > 0: \vec h = 0 é o equilíbrio e \bar h oscila com \sqrt{K_0/m}.
• Na Aula 2, K_0 dará o termo de massa: um campo sem massa não tem valor preferido.

Notas: O modo zero

Altura média

Figura 11: Altura média \bar h(t) da cadeia (N = 15, m = K = 1, P(0) = 1, \bar h(0) = 1/15). Com K_0 = 0, \bar h cresce linearmente até 4{,}07 em t = 60 e coincide com a reta \bar h(0) + \bar p\, t/m (tracejada), com diferença máxima 1{,}3\times 10^{-14}. Com K_0 = 0{,}2, \bar h oscila com frequência \sqrt{K_0/m} = 0{,}447 e coincide com \bar h(0)\cos(\sqrt{K_0/m}\,t) + \bar p(0)\sin(\sqrt{K_0/m}\,t)/\sqrt{mK_0} (tracejada), com diferença máxima 8{,}1\times 10^{-5}.

Notas: Simulação

Simulação

def forca_cadeia(h, K, K0):
    return K * (np.roll(h, -1) - 2 * h + np.roll(h, 1)) - K0 * h

• O mesmo passo_verlet do oscilador, com h e p vetores.
• N = 15, m = K = 1; h^0(0) = 1, partículas 1 e 2 com velocidade 0{,}5: P = 1, \bar p = 1/15.
• K_0 = 0 e K_0 = 0{,}2; \Delta t = 0{,}05, 1200 passos.

Notas: Simulação

Animação

Once Loop Reflect

Linha tracejada: \bar h(t). Com K_0 = 0 a cadeia inteira sobe; com K_0 = 0{,}2 ela oscila em torno de h = 0.

Notas: Simulação

Desvios em relação à média

Figura 12: Desvio h^j(t) - \bar h(t) de cada partícula em relação à altura média, para K_0 = 0 (esquerda) e K_0 = 0{,}2 (direita), durante 0 \le t \le 60. A perturbação inicial nas partículas 0, 1 e 2 se propaga nos dois sentidos ao redor do anel; a partícula j = 7, no lado oposto, se desloca mais de 0{,}05 pela primeira vez em t = 4{,}6 nos dois casos. Os desvios ficam entre -0{,}70 e 0{,}93.

Notas: Simulação

Matrizes circulantes e a transformada de Fourier discreta

Matrizes circulantes

C = \begin{pmatrix} b_0 & b_1 & b_2 & b_3 \\ b_3 & b_0 & b_1 & b_2 \\ b_2 & b_3 & b_0 & b_1 \\ b_1 & b_2 & b_3 & b_0 \end{pmatrix}, \qquad C^j{}_k = b_{k-j}, \qquad C = \sum_{r=0}^{N-1} b_r\, T^r

• Cada linha é a anterior deslocada para a direita. Exemplo: KD + K_0 I, com b_0 = 2K + K_0 e b_1 = b_{N-1} = -K.
• Circulantes comutam entre si: são polinômios em T.
• C é circulante se, e somente se, CT = TC: a regra é a mesma em todos os pontos do anel (invariância por translação).

Notas: Matrizes circulantes

Autovetores de T: a base de Fourier

T\vec v = \lambda\vec v \;\Rightarrow\; \vec v = T^N\vec v = \lambda^N\vec v \;\Rightarrow\; \lambda^N = 1, \qquad v^{j+1} = \lambda v^j \;\Rightarrow\; v^j = \lambda^j v^0, \qquad |v^0| = \frac{1}{\sqrt N}, \quad v^0 > 0 v_a^{\,j} \equiv \frac{e^{2\pi i a j/N}}{\sqrt N}, \qquad T \vec v_a = e^{2\pi i a/N} \vec v_a, \qquad \vec v_a^{\,\dagger} \vec v_b = \delta_{ab}

• A fase de v^0 é arbitrária; escolhemos v^0 real e positivo.

Notas: A base de Fourier

A base de Fourier

Figura 13: Esquerda: os N = 15 autovalores e^{2\pi i a/N} de T no círculo unitário, rotulados por a no intervalo simétrico. Centro: parte real de \sqrt N\, v_a^{\,j} = e^{2\pi i a j/N} para a = 1 e a = 2 (pontos) e as funções \cos(2\pi a x/N) (linhas). Direita: \cos(2\pi a x/N) para a = 7 e a = 7 - 15 = -8 passam pelos mesmos pontos da rede: nos pontos j inteiros, \vec v_7 = \vec v_{-8}.

Notas: A base de Fourier

A transformada de Fourier discreta

\vec h = \sum_a \tilde h^a\, \vec v_a, \qquad \tilde h^a \equiv \vec v_a^{\,\dagger} \vec h = \frac{1}{\sqrt N}\sum_{j=0}^{N-1} e^{-2\pi i aj/N}\, h^j, \qquad h^j = \frac{1}{\sqrt N}\sum_a e^{2\pi i aj/N}\,\tilde h^a

• a é definido módulo N: \vec v_{a+N} = \vec v_a. Usamos o intervalo simétrico, o de menor |a|.
• \vec h real: \tilde h^{-a} = \left(\tilde h^a\right)^*.
• np.fft.fft(h, norm="ortho") dá \tilde h^a com O(N\log N) operações.

Notas: A base de Fourier

Toda circulante é diagonal na base de Fourier

C\, \vec v_a = \sum_r b_r\, T^r \vec v_a = \lambda_a(C)\, \vec v_a, \qquad \lambda_a(C) = \sum_{r=0}^{N-1} b_r\, e^{2\pi i a r/N}, \qquad U^\dagger C\,U = \operatorname{diag}\big(\lambda_a(C)\big) \lambda_a(D) = 2 - e^{2\pi i a/N} - e^{-2\pi i a/N} = 2 - 2\cos\frac{2\pi a}{N} = 4\sin^2\frac{\pi a}{N}

• Os autovetores são os mesmos para toda circulante; os autovalores são a transformada de Fourier da primeira linha.
• \lambda_a(D) = 0 somente para a = 0: \vec v_0 = \vec 1/\sqrt N, o modo zero.

Notas: Toda circulante é diagonal

A relação de dispersão

\omega_a = \sqrt{\frac{K_0 + 4K\sin^2(\pi a/N)}{m}}

Figura 14: Relação de dispersão da cadeia para N = 15, m = K = 1, com K_0 = 0 e K_0 = 0{,}2. As linhas são \omega_a com a contínuo, os pontos cheios são os 15 valores inteiros de a, e os círculos vazios são as raízes dos autovalores de (KD + K_0 I)/m calculados numericamente, que coincidem com os pontos cheios com diferença máxima 2{,}3\times 10^{-15}. Os máximos são \omega_{\max} = 1{,}989 e 2{,}039, abaixo de \sqrt{(K_0 + 4K)/m} = 2 e 2{,}049 porque N é ímpar.

Notas: A relação de dispersão

Modos normais

m\,\ddot{\tilde h}^a = -m\,\omega_a^2\,\tilde h^a, \qquad \tilde h^a(t) = \tilde h^a(0)\cos\omega_a t + \frac{\tilde p^a(0)}{m\,\omega_a}\sin\omega_a t H = \sum_a \left(\frac{|\tilde p^a|^2}{2m} + \frac{m\,\omega_a^2\,|\tilde h^a|^2}{2}\right)

• N osciladores acoplados \to N osciladores independentes, um por modo.
• a = 0, K_0 = 0: \tilde h^0 = \sqrt N\,\bar h e \tilde p^0 = \sqrt N\,\bar p dão \bar h(0) + \bar p\,t/m, o modo zero.
• Para h real, N números reais: \tilde h^0 e os pares (a, -a).

Notas: Modos normais

A base real de Hartley

\zeta_a^{\,j} \equiv \frac{1}{\sqrt N}\left(\cos\frac{2\pi a j}{N} + \sin\frac{2\pi a j}{N}\right) = \frac{(1 - i)\,v_a^{\,j} + (1 + i)\,v_{-a}^{\,j}}{2}, \qquad \vec\zeta_a^{\,\mathsf T}\vec\zeta_b = \delta_{ab}

Once Loop Reflect

Os quatro modos com o mesmo relógio: cada um oscila com a sua \omega_a, sempre com a mesma forma.

• \vec\zeta_a\cos\omega_a t: onda estacionária. \vec v_a e^{-i\omega_a t}: onda viajante, para j crescente se a > 0.

Notas: A base real de Hartley

Solução exata pela FFT

ht0 = np.fft.fft(h0, norm="ortho")
pt0 = np.fft.fft(p0, norm="ortho")
ht = ht0 * np.cos(w * t) + pt0 / m * t * np.sinc(w * t / np.pi)
h = np.fft.ifft(ht, axis=1, norm="ortho").real

• FFT, evolução exata de cada modo, FFT inversa: O(N\log N) para qualquer t.
• t * np.sinc(w * t / np.pi) é \sin(\omega_a t)/\omega_a, igual a t para \omega_a = 0.
• Verlet com \Delta t dividido por 2: erro dividido por 4{,}0 (segunda ordem).

Notas: Solução exata pela FFT

Exata \times Verlet

Figura 15: Solução exata pela FFT comparada com o Verlet (\Delta t = 0{,}05) para as condições iniciais da simulação da seção anterior. Esquerda: altura h^7(t) da partícula do lado oposto do anel, com K_0 = 0{,}2; a solução exata (linha cheia) e o Verlet (tracejada) coincidem na escala do gráfico. Direita: \max_j |h^j_{\rm Verlet} - h^j_{\rm exata}|, que cresce linearmente com t e chega a 1{,}15\times 10^{-2} (K_0 = 0) e 1{,}22\times 10^{-2} (K_0 = 0{,}2) em t = 60.

Notas: Solução exata pela FFT

Estabilidade do Verlet na cadeia

\Delta t < \frac{2}{\omega_{\max}}, \qquad \omega_{\max} = \max_a\omega_a \le \sqrt{\frac{K_0 + 4K}{m}}

Figura 16: Energia da cadeia (N = 15, m = K = 1, K_0 = 0{,}2, condições iniciais da simulação da seção anterior) em 300 passos de Verlet, 2\% abaixo e 2\% acima do limite \Delta t = 2/\omega_{\max} = 0{,}981. Abaixo do limite, E/E_0 oscila entre 0{,}51 e 1{,}04. Acima, apenas os modos a = \pm 7 têm \omega_a\Delta t > 2, e a energia cresce por um fator 2{,}2226 por passo, igual a \lambda^2 para o autovalor \lambda de maior módulo da matriz do Verlet com \varepsilon = \omega_7\Delta t; após 300 passos, E/E_0 = 7\times 10^{102}.

Notas: Estabilidade do Verlet na cadeia