Introdução à dinâmica molecular clássica

só um divisor

Osciladores acoplados

Um dos problemas mais clássico de mecânica, mecânica estatística e estado sólido é encontrar a curva de dispersão de um sistema de infinitos osciladores harmônicos acoplados

O Hamiltoniano que descreve o sistema está dado por \[ H = \sum_l \frac{p_l^2}{2m}+\frac{1}{2}\beta\left( x_l - x_{l+a}\right)^2 = \sum_l H_l \nonumber \] onde $i$ diz respeito ao índice da partícula, $\beta$ é a constante de força (lei de Hook), $m$ é a massa, $x_l$ é a posição da partícula e $a$ é a separação de equilíbrio entre as partículas. A partir desta equação podemos calcular a equação de movimento: \[ \frac{dp_l}{dt} = -\frac{\partial H_l}{\partial x_l} = -\beta\left( x_l - x_{l+a} \right) + \beta\left( x_{l-a} - x_l \right) \nonumber \] e \[ \frac{dx_l}{dt} = \frac{\partial H_l}{\partial p_l} = \frac{p_l}{m} \nonumber \] de onde \[ m\frac{d^2x_l}{dt^2} = -k\left( 2x_l - x_{l+a} - x_{l-a} \right) \nonumber \]

Condições periódicas

Uma forma bem conveniente de simular o infinito é assumir as seguintes condições: \[ x_{Na} = x_0 \nonumber \] onde $N$ é o número de osciladores. Observe que a temperatura zero as partículas estão nas posições $x_a,\,x_{2a},\,x_{3a},\,\ldots,\,x_{Na}$

Modos normais de vibração

Vamos propor a seguinte solução \[ x_l = \sum_k x_k e^{ikl} \nonumber \] de onde \[ m e^{ikl} \frac{d^2x_k}{dt^2} = -\beta \left[ 2x_k e^{ikl} - x_k e^{ik\left(l+a\right)} - x_k e^{ik\left(l-a\right)} \right] \nonumber \] de onde \[ m \frac{d^2x_k}{dt^2} =-2k\left( 1 - \frac{e^{ia} - e^{-ia}}{2} \right)x_k \nonumber \] ou \[ m \frac{d^2x_k}{dt^2} =-2k\left(1-\cos ka\right)x_k \nonumber \] esta é a equação de um oscilador harmônico na qual a frequência angular está dada por \[ \omega_k = \sqrt{\frac{2k}{m}\left(1-\cos ka\right)} \nonumber \] ou \[ \omega_k = 2 \sqrt{\frac{\beta}{m}}\;\; \left|\sin \frac{ka}{2}\right| \label{relDisp} \] de onde \[ x_k = x_k\left(0\right) e^{-i\omega_k t} \nonumber \] de forma que \[ x_l = \sum_k x_k\left(0\right) e^{-i\omega_k t}\;\; e^{ikl} \nonumber \] mas, as condições periódicas exigem \[ e^{ik \left(Na\right) } = e^{ik \left(0\right) }=1 \nonumber \] como $e^{i2\pi\,l} = 1$ ($l$ inteiro), então \[ k = \frac{2\pi}{Na}l \label{vecOnda} \]

A equação $\ref{relDisp}$ se conhece como relação de dispersão e diz qual é a dependência do vetor de onda $q$ com a frequência angular $\omega$. E a sua forma está dada pela figura acima.

No limite de comprimentos de onda largos (ou $k\ll \pi/a$) temos que \[ \begin{align*} \omega_k =& 2 \sqrt{\frac{\beta}{m}}\;\; \left(\frac{ka}{2}\right)\\ =& \sqrt{\frac{a^2\beta}{m}} \; k \end{align*} \] quando $k = \pi /a$ temos a chamada frequência de Debye dada por \[ \omega_D = \sqrt{\frac{\beta}{m}} \pi \nonumber \] que se obtém supondo o sistema contínuo e não discreto como se fez aqui. Observe que a máxima frequência que o sistema adquire é tal que \[ \omega_0 = \frac{\pi}{2} \omega_D \nonumber \]

Densidade de estados

Da figura acima vemos que \[ \begin{align*} g(\omega)d\omega =& 2g(k)dk\\ =& 2g(k)\frac{dk}{d\omega}d\omega \end{align*} \] onde $g$ são as distribuições de $\omega$ e $q$. Da equação $\ref{vecOnda}$: \[ \Delta k = \frac{2\pi}{Na}\left[ (l+1) - l\right] = \frac{2\pi}{Na} \nonumber \] assim a quantidade de valores que $k$ assume, $g(k)$, numa unidade de intervalo de $k$ é \[ g(k) = \frac{Na}{2\pi}\;\;\;\;\;\;\;\;-\frac{\pi}{a} < q < \frac{\pi}{a} \nonumber \] já que deve se verificar que \[ g(k)\Delta k = 1 \nonumber \] de $\ref{relDisp}$ temos que \[ k = \frac{2}{a} \arcsin\left( \frac{\omega}{\omega_{0}} \right) \nonumber \] onde \[ \omega_{0} = 2 \sqrt{\frac{\beta}{m}} \nonumber \] dessa forma \[ \frac{dk}{d\omega} = \frac{2}{a\left( \omega^2_0 - \omega^2 \right)^{1/2}} \nonumber \] substituindo \[ g(\omega)d\omega = \frac{N}{\pi}\frac{1}{\left( \omega^2_0 - \omega^2 \right)^{1/2}}d\omega \] note que \[ g_0 = \lim_{\omega \rightarrow 0} g(\omega) = \frac{2N}{\pi\omega_0} \nonumber \]

O número de modos de vibração em que uma molécula linear pode vibrar são $3N-5$, o 5 se deve às coordenadas do centro de massa $(x,y,z)$ + duas formas de rotar. Se o movimento é restrito no plano o número de modos de vibração são $N-1$, o 1 se deve à coordenada do centro de massa.

"Fonons" clássicos em Dinâmica Molecular

Densidade de estado

Nos parágrafos prévios mostramos qual é a forma como as frequências de uma cadeia de osciladores clássico se distribui. Nesta seção vamos calcular esse mesmo comportamento a partir das velocidades das partícula da cadeia.

No artigo com DOI: 10.1103/PhysRev.188.1407 Dickey e Paskin mostram como calcular a densidade de estados de um sistema de partículas de uma forma muito simples: Sabemos que para um oscilador de amplitude $A$, fase $\phi$ e frequência $\omega$ é valido que \[ \begin{align*} x =& A \cos \left( \omega t + \phi \right)\\ v = & A \sin \left( \omega t + \phi \right) \end{align*} \] e que no caso de um oscilador clássico sua energia está dada por \[ \frac{1}{2}mA^2\omega^2 = k_BT \nonumber \] Note que se calculamos a autocorrelação de velocidades definida como \[ \gamma (t) = \frac{\left<\sum_i v_i\left(t\right)\cdot v_i\left(0\right)\right>}{\left<\sum v_i^2\left(0\right)\right>} \] obteremos \[ \begin{align*} \gamma (t) =& \frac{\left< A^2 \omega^2 \sin \left( \omega t + \phi \right) \sin \phi \right>}{\left< A^2 \omega^2 \sin^2 \phi \right>}\\ =& \frac{\left<\sin \left( \omega t + \phi \right) \sin \phi \right>}{\left< \sin^2 \phi \right>}\\ \end{align*} \] como \[ \left< \sin^2 \phi \right> = \frac{1}{2} \nonumber \] e \[ \begin{align*} \left< \sin \left( \omega t + \phi \right) \sin \phi \right> =& \frac{1}{2}\left< \cos \omega t - \cos \left( \omega t + 2\phi \right)\right>\\ =& \frac{1}{2}\cos \omega t \end{align*} \] então \[ \gamma (t) = N^{-1}\sum_i \cos \omega t \nonumber \] que, passando para o contínuo fica \[ \gamma (t) = N^{-1}\int \cos \omega t\, dt \nonumber \] Note que essa equação acima é adequada para quando temos $1$ oscilador, no caso de mais osciladores devemos esperar \[ \gamma (t) = \int f\left(\omega \right)\cos \omega t\; dt \] onde $f\left(\omega \right)$ é a probabilidade de um dado oscilador ter frequência $\omega$. Note que que $\gamma (t)$ é a transformada inversa de fourier de $f\left(\omega \right)$.

Curva de dispersão

Sabemos da Física IV que quando um feixo de luz com comprimento de onda na região dos Raios X é jogada sobre uma superfície cristalina obtemos como luz refletida a formação de um espetro de difração onde a posição dos máximos desse espetro é descrita pela conhecida lei de Bragg: \[ 2d\sin\theta = n\lambda \nonumber \] Essa descrição simplista da interação da fóton-fônon pode ser descrita de uma forma mais completa considerando o comportamento quântico das partículas no processo de espalhamento, nesse sentido se $\left| \vec{k} \right>$ e $\left| \vec{k}' \right>$ são os estados próprios da onda incidente e refletida, respetivamente. Segundo a regra de ouro de Fermi, a taxa de transição entre estes estados estará dada por \[ M_{\vec{k},\vec{k}'} = \int d e^{-i\vec{k}\cdot \vec{r}}U\left( \vec{r} \right)e^{i\vec{k}\cdot \vec{r}} \nonumber \] onde $U\left( \vec{r} \right)$ é o potencial de espalhamento e $\left< \vec{r} {\left| \vec{k} \right.} \right> = e^{i\vec{k}\cdot \vec{r}}$. Pode ser demostrado que a taxa de transição está relacionada à seção de choque ou seção de espalhamento diferencial \[ \frac{d\sigma}{d\Omega} \sim \frac{2\pi}{\hbar} \left| M_{\vec{k},\vec{k}'} \right| \nonumber \] Se assumimos que o potencial de interação é igual à soma dos potenciais iguais que formam o sistema \[ U\left( \vec{r} \right) = \sum_\alpha U_\alpha \left( \vec{r} - \vec{r}_\alpha \right) \nonumber \] pode se mostrar que \[ M_{\vec{k},\vec{k}'} = \left< {\vec{k} \left|U\left( \vec{r}\right)\right| \vec{k}'}\right> = \sum_\alpha U_\alpha \left( \vec{q} \right) e^{i\vec{q}\cdot \vec{r}_\alpha} \nonumber \] onde, $\vec{q} = \vec{k} - \vec{k}'$. $U_\alpha \left( \vec{q} \right)$ é conhecido como fator de forma atômico o qual não é outra coisa do que a transformada de Fourier do potencial. Se o fator de forma atômico é o mesmo para todas as partículas, então \[ \frac{d\sigma}{d\Omega} \sim \left| U \left( \vec{q} \right) \right|I\left( \vec{q} \right) \nonumber \] onde $I\left( \vec{q} \right)$ é a função de estrutura dada por \[ I\left( \vec{q} \right) = \left< {\sum_{\alpha,\alpha'} e^{i\,\vec{q} \cdot \left( \vec{r}_\alpha - \vec{r}_{\alpha '} \right)}} \right> \nonumber \] a qual, ao ser divida por $N$ ou pelo volumem nos da o fator de estrutura $S\left( \vec{q} \right)$ \[ S\left( \vec{q} \right) = N^{-1}\left< {\sum_{\alpha,\alpha'} e^{i\,\vec{q} \cdot \left( \vec{r}_\alpha - \vec{r}_{\alpha '} \right)}} \right> \] A validade do resultados até obtidos está determinada pela regra de ouro de fermi, isto é, interações fracas entre a onda incidente e o alvo ou, equivalentemente, livre caminho meio maior do que a distância inter partícula do alvo.

A partir da lei de Bragg sabemos qual é o comprimento caraterístico que pode ser sondado com um dada onda incidente de comprimento de onda $\lambda$, assim, para a luz visível ($\epsilon \sim 1\,eV \Rightarrow \lambda \sim 0.4 - 0.7\overset{\circ}{A} $) permite enxergar estruturas em escala de mícrons; com raio X podemos penetrar até $1mm$ e obter informação do "bulk", com elétrons ($\epsilon \sim 100\,eV \Rightarrow \lambda \sim 1\overset{\circ}{A} $) temos informação da nuvem eletrônica, com nêutrons ($\epsilon \sim 0.1\,eV \sim 400K \Rightarrow \lambda \sim 1\overset{\circ}{A} $) conseguimos ver o núcleo e os spins eletrônicos.

Funções de correlação

Vamos tentar obter a partir das funções de correlação o fator de estrutura, para isso definimos o operador densidade de partículas por unidade de volume como \[ \rho\left(\vec{r}\right) = \sum_\alpha \delta\left(\vec{r} - \vec{r}_\alpha\right) \nonumber \] o qual resulta, no caso de um fluido isotrópico, igual à densidade média $\rho = N/V$. Também definimos a densidade de correção em dois pontos como \[ \begin{align*} C_{nn}\left(\vec{r}_1,\vec{r}_2\right) =& \left< {\rho\left(\vec{r_1}\right), \rho\left(\vec{r_2}\right)} \right>\\ =& \left< {\sum_{\alpha,\alpha '}\, \delta\left(\vec{r}_1 - \vec{r}_\alpha\right) \delta\left(\vec{r}_2 - \vec{r}_{\alpha '}\right) } \right> \end{align*} \] Se o sistema possui invariância translacional então é verdade que \[ C_{nn}\left(\vec{r}_1,\vec{r}_2\right) = C_{nn}\left(\vec{r}_1 - \vec{r}_2\right) \nonumber \] de donde, se calculamos a transformada de Fourier dessa grandeza teremos \[ \begin{align*} I\left( \vec{q} \right) =& \int d\vec{r}_1 d\vec{r}_2\, e^{-i\vec{q}\cdot \left(\vec{r}_1 - \vec{r}_2\right)} \left< {\rho\left(\vec{r}_1\right), \rho\left(\vec{r}_2\right)} \right>\\ =& \left< {\rho\left(\vec{q}\right), \rho\left(-\vec{q}\right)} \right> \end{align*} \] onde \[ \begin{align*} \rho\left(\vec{q}\right) =& \int d\vec{r} e^{-i\vec{q}\cdot \vec{r}} \rho\left(\vec{r}\right) \\ =& \sum_\alpha e^{-i\vec{q}\cdot \vec{r}_\alpha} \end{align*} \] de onde vemos que os experimentos de espalhamento medem diretamente as funções de correlação.

De entre as varias funções de correlação que podem ser medidas, a função de correlação de pares merece uma atenção especial ou função de van Hole. Essa correlação se define como \[ \begin{align*} g(\vec{r},t) =& \left< { \frac{1}{N} \sum_i^N \sum_{j\neq i}^N\int \delta\left[ \vec{r} - \vec{r}_j + \vec{r}(0) \right]} \right>\\ =& \left< { \frac{1}{N}\int d\vec{r}'\sum_i^N \sum_{j\neq i}^N\delta\left[ \vec{r}' + \vec{r} - \vec{r}_j(t) \right] \delta\left[ \vec{r}' -\vec{r}_j(0) \right]} \right>\\ =& \left< { \frac{1}{N} d\vec{r}' \rho\left( \vec{r}'+\vec{r},t \right) \rho\left( \vec{r}',0 \right) } \right>\\ =& \frac{1}{\rho} \left< { \rho\left( \vec{r}(t),t \right) \rho\left( \vec{r}(0),0 \right) } \right> \end{align*} \] \[ \begin{align*} \left< \rho \right> g\left(\vec{r}_1, \vec{r}_2\right) =& \frac{1}{V}\int d\vec{r}_2 \left< { \sum_{\alpha\neq \alpha '} \delta\left( \vec{r}_1 - \vec{r}_\alpha \right) \delta\left( \vec{r}_2 - \vec{r}_{\alpha '} \right)} \right>\\ =& \frac{1}{V} \left< { \int d\vec{r}_2 \sum_{\alpha\neq \alpha '} \delta\left( \vec{r} + \vec{r}_2 - \vec{r}_\alpha \right) \delta\left( \vec{r}_2 - \vec{r}_{\alpha '} \right)} \right>\\ =& \frac{1}{V} \left< { \sum_{\alpha\neq \alpha '} \delta\left( \vec{r} + \vec{r}_\alpha - \vec{r}_{\alpha '} \right) } \right> \end{align*} \] onde $\vec{r} =\vec{r}_1 - \vec{r}_2$. Como as soma em $\alpha$ varre todos os valores para cada $\alpha '$, então todos os termos da soma em $\alpha$ são iguais, então \[ g\left(\vec{r}\right) = \frac{1}{\left< \rho \right>}\left< { \sum_{\alpha\neq 0} \delta\left( \vec{r} + \vec{r}_\alpha - \vec{r}_{0} \right) } \right> \label{gr} \] assim esta função representa a probabilidade de encontrar uma partícula a uma distancia $\vec{r}$ de outra na posição $\vec{r}_0$, normalizada pela densidade média.

Para fluidos homogêneos podemos mostrar que \[ S(\vec{q}) = \left< \rho \right>\int d \vec{r} e^{-i\,\vec{q}\cdot \vec{r}} g(\vec{r}) \]

Se levamos em consideração o tempo, a equação $\ref{gr}$ se modifica para \[ g\left(\vec{r},t\right) = \frac{1}{V}\left< { \sum_{\alpha\neq \alpha '} \delta\left( \vec{r} + \vec{r}_\alpha(t) - \vec{r}_{\alpha '}(0) \right) } \right> \] a partir desta equação podemos calcular a função intermediaria de dispersão como \[ F\left( \vec{k},t \right) = \int g\left(\vec{r},t\right) e^{-i\,\vec{k}\cdot \vec{r}} d\vec{r} \] e desta, por sua vez, calculamos o fator de estrutura dinâmico: \[ S( \vec{k},\omega) = \int F\left( \vec{k},t \right) e^{-i\,\omega\, t} dt \] Do fator de estrutura dinâmico podemos calcular a curva de dispersão simplesmente plotando $S(\vec{k},\omega)$ vs $k$, para diversos $\omega$ dessas figuras escolhemos os maiores picos (assim fixam o $k$) e depois se faz o gráfico de $\omega\, vs\, k$

Análise de Fourier com Python

Serie de Fourier

Devido ao fato de que funções $\sin$ e $\cos$ formam um conjunto completo e ortogonal podemos utilizar elas como base e expressar outras funções em serie de senos e coseno, sendo assim podemos escrever \[ f(t) = \frac{a_0}{2}+\sum_{n=1}^\infty a_n \cos nt + \sum_{n=1}^\infty b_n \sin nt \nonumber \] onde \[ a_n = \frac{1}{\pi}\int_{-\pi}^{\pi} f(t)\cos nt\, dt \nonumber \] \[ b_n = \frac{1}{\pi}\int_{-\pi}^{\pi} f(t)\sin nt\, dt \nonumber \]

Transformada de Fourier

A expansão em serie de Fourier é muito adequada quando a função é periódicas num intervalo finito ou infinito, mas caso a função seja não periódica devemos utilizar a transformada de Fourier num intervalo infinito. Por definição a transformada de Fourier está dada por \[ F(y) = A\int_{-\infty}^{\infty} f(x) e^{-ixy}dx \nonumber \] e a transformada inversa \[ f(x) = \frac{1}{2\pi A}\int_{-\infty}^{\infty} F(y) e^{-ixy}dy \nonumber \] onde $A$ depende do autor, nos aqui suporemos $A=1/\sqrt{2\pi}$. O algoritmo numérico utilizado para o calculo da transformada de Fourier se conhece como "fast Fourier Transform" ou simplesmente fft.

Transformada de Fourier discreta - DFT

Caso tenhamos dados discretos, por exemplo, dados que foram capturados a cada certo $\Delta t$, de forma que $f(m\Delta t) = 0, 1, \ldots , N-1$ e tenhamos certeza que nesse intervalo está nosso fenômeno de interesse. No caso da serie de Fourier podemos aproximar ela no intervalo $0 < t < T$ pela seguinte equação \[ f(t) = \sum_{-\infty}^{\infty} c_n e^{i2\pi nt/T} \nonumber \] onde \[ c_n = \frac{1}{T}\int_0 ^T f(t) e^{i2\pi nt/T} dt \nonumber \] observe que esta função é periódica em $T$. Afim de aproximar este resultado ao da transformada de Fourier continua fazemos \[ \Delta \omega = \frac{2\pi}{T} \nonumber \] A transformada discreta de Fourier está dada por \[ g(n\Delta \omega) = \sum_{0}^{N-1} f(m\Delta t) e^{-i2\pi mn/N} \nonumber \] onde temos $N$ dados conhecidos e ao fazer a transformada a maior frequência que obteremos é $(N-1)\Delta \omega$. A transformada inversa discreta está dada por \[ f(m\Delta t) = \frac{1}{N} \sum_{0}^{N-1} g(n\Delta \omega) e^{i2\pi mn/N} \label{tfdiscreta} \]

A transformada rápida de Fourier - FFT

Este é um método que permite reduzir em muito o tempo de calculo da transformada discreta de Fourier, enquanto a transformada discreta é da ordem de $N^2$ a transformada rápida é da ordem $N\log_2N$. O algoritmo descompõe os dados de forma a constituir grupo de dados pares e impares, por isso esse algoritmo trabalha bem utilizando dados tais que $N= 2^k$, ou dito simplificadamente, dados que são potencias de $2$.

No caso do Numpy, a função que calcula a transformada rápida de Fourier é a função np.ftt.ftt e se os dados são todos reais a função np.ftt.rfft melhora a performance. Lembre que esses algoritmos são otimizados para trabalhar com $2^k$ dados, isto é, um número par de elementos. Se $A=np.ftt.ftt(a, n)$ então $A[1:n/2]$ tem os termos com frequência positiva e $A[n/2+1:]$ tem os termos com frequência negativa. Se os dados da variável independente de entradas estão espaçados em $\Delta t$, então a máxima frequência que se obtém a partir da transformada é $f_{max} = 1.0/2\Delta t$ e a mínima frequência é $\Delta f = 1.0/2N\Delta t$ o qual é também o intervalo em frequência do sistema, ou seja, as frequência amostradas estão dadas por \[ f_i = \frac{i}{N\Delta t}\;\;\;\;\;\;\;\;\;\;i=-\frac{N}{2},\ldots ,\frac{N}{2} \] No caso do numpy, o espectro de amplitude está dado por np.abs(A) enquanto que o espetro em potencias está dado por np.abs(A)**2 e o espetro de fase é np.angle(A), se a é um sinal temporal em A = np.ftt.rfft(a).

Exemplos

No exemplo a seguir vamos calcular a transformada rápida de Fourier para obter as frequência caraterísticas da função \[ y = \sin\left(2\pi f_1 t \right) + \sin\left(2\pi f_2 t + \phi\right) \nonumber \] com $f_1 = 50$, $f_2 = 70$ e $\phi = \pi / 4$

      #/usr/bin/env python
      #-*- coding: utf-8 -*-
      
      import numpy as np
      import matplotlib.pyplot as plt
      
      pi = np.pi
      
      tmax    = 0.5
      Npontos = 500
      dt      = 1.0 / Npontos
      df      = 1.0 / tmax
      
      t   = np.arange(0, tmax, dt)
      
      y = np.sin (2* pi *50* t )+ np.sin (2* pi *70* t + pi /4)
      
      a = np.fft.fft( y )
      
      #freqs = df * np.arange(0, (n_t -1)/2. ,dtype = 'float')
      freqs  = np.fft.fftfreq(len(f), dt)
      
      plt.subplot (2, 1, 1)
      plt.plot (t, y, label = ' data de entrada ')
      plt.xlabel ('time [ s ] ')
      plt.ylabel ('y ')
      
      plt.subplot (2 ,1 ,2)
      plt.plot(freqs , abs(a), \
         label = 'abs ( Transformada de Fourier )'
      plt.xlabel ( u'frequência ')
      plt.ylabel ( 'abs ( FFT ( y )) ')
      
      plt.savefig ( 'fig_fft.png')
      plt.show ()
      

O resultado da transformada está plotado na figura acima, observe que aparecem 4 frequências, isso acontece porque como foi dito anteriormente o intervalo de interesse é de $[0,n-1/2]$, além desse intervalo se repetem as mesmas frequência mas agora com sinal negativo.

No seguinte exemplo vamos utilizar a FFT do python para aproximar a função \[ y = \left( 1 + x - x^3 \right) e^{-x^2} \nonumber \]

      #/usr/bin/env python
      #-*- coding: utf-8 -*-

      import matplotlib.pyplot as plt
      import numpy as np

      N = 2**11
      x = np.linspace(-10,10,N+1)
      a = x[2] - x[1]
      L = a*N
      f_x = (1 + x - x*3)*np.exp(-x**2)
      a_k = np.fft.rfft(f_x)
      k = 2*np.pi/L*np.arange(len(a_k))

      x_des = np.fft.fftshift(x)

      #f_nova_x = np.zeros(len(x), dtype='complex')
      #for i in range(len(x)):
        #for l in range(len(k)):
          #f_nova_x[i]  = f_nova_x[i] + a_k[l] * np.exp(1j*k[l]*x[i])

      f_nova_x = np.sum( ( np.outer(np.ones(len(x)), a_k) * \
             np.exp( np.outer(x, 1j*k) ) ), axis=1)
      
      f_nova_x = f_nova_x / len(k)


      plt.figure( figsize=(14,5) )
      plt.subplot(131)
      plt.plot(x, f_x)
      plt.xlabel('x')
      plt.ylabel('f_x')

      plt.subplot(132)
      plt.plot(x, f_nova_x.real)
      plt.xlabel('x')
      plt.ylabel('f_nova_x')

      plt.subplot(133)
      plt.plot(x_des, f_nova_x.real)
      plt.xlabel('x_deslocado')
      plt.ylabel('f_nova_x')

      plt.tight_layout()

      plt.savefig('fig07.png')
      plt.show() 
      

Para aproximar uma função utilizando serie de Fourier discreta utilizamos a equação $\ref{tfdiscreta}$. No script observe que foi utilizado o produto externo de matrizes para gerar os dados, tudo foi feito em uma linha só - linha 22 -, mas essa linha é equivalente ao dos loop entre as linhas 17 e 20 que estão comentadas. Observe que no gráfico do meio plotamos o valor aproximado de y (f_new_x) vs x e isso deu uma curva esquisita que lembra a curva original, mas já dize, os valores de frequência são tais que da metade para a frente dos dodos correspondem às frequência negativas então devemos mover da metade dos dados até o fim para o inicio do vetor, isso se faz com a função np.ftt.fftshift - linha 15 - que ordena o eixo $x$ para que este em acordo com as frequência calculadas.

Tarefas

  1. Crie um programa que aproxime a função de onda quadrada mediante $n$ termos da sua serie de Fourier
  2. Crie um programa que calcule a transformada de Fourier de uma função dada por \[ y = \sum_i^n A_i \sin\left(k_i\,x\right) \nonumber \] onde $A_i$ é um número aleatório entre $[0,1)$ e \[ k_i = \frac{2\pi}{\xi _i} \nonumber \] em que $\xi_i$ é um número aleatório entre $(0, 200]$.
  3. Dada a função do problema anterior calcule a sua autocorrelação: \[ \gamma (x) = \frac{\left<\sum_i y_i\left(x\right)\cdot y_i\left(0\right)\right>}{\left<\sum y_i^2\left(0\right)\right>} \nonumber \]
  4. Faça o gráfico da função \[ sinc(x) = \frac{\sin(x)}{x} \nonumber \] e calcule a sua transformada de Fourier e faça o gráfico.
  5. Dada a equação diferencial de Van der Pol calcule sua
  6. Usando o algoritmo de Runge-Kutta crie um programa que calcule a densidade de estados de um conjunto de $n$ osciladores acoplados nos quais $x_0=x_{n-1}$ (fixo nos extremos).