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

só um divisor

Resolução de equações diferencial de II grau

Se uma equação contem as derivadas ou diferenciáveis de uma ou mais variáveis dependente com respeito a uma ou mais variáveis independentes, então se diz que é uma equação diferencial, Exemplo \begin{equation} \frac{d^2x}{dt^2} = \frac{F}{m} \nonumber \end{equation} A equação acima é uma equação diferencial de II ordem. Provavelmente seja a equação diferencial mais famosa que existe em Física. Futuramente vocês aprenderam métodos matemáticos que permitiram resolver analiticamente essa equação (quando possível, é claro) mas neste mini curso aprenderemos a tratar esse problema de forma numérica. No caso de essa equação as vezes é possível reduzir uma ordem na equação e obtendo assim duas equações de I ordem: \begin{equation} \frac{dv}{dt} = \frac{F}{m} \nonumber \end{equation} e \begin{equation} \frac{dx}{dt} = v \nonumber \end{equation} Assim vamos ver alguns método que nos permitam resolver essas equações de I ordem

Método de Euler

Consideremos a equação diferencial \begin{equation} \begin{array}{rcl} \frac{dv}{dt} &= &a(r,v)\\ \frac{dr}{dt} &= &v(t) \end{array} \nonumber \end{equation} a solução dessa equação diferencial mediante a aplicação a técnica de derivada para adiante mostrada anteriormente, isto é \begin{equation} \begin{array}{rcl} \frac{v(t+\Delta t)-v(t)}{\Delta t} & = & a(r,t)\\ \frac{r(t+\Delta t)-r(t)}{\Delta t} & = & v(t)\\ \end{array} \nonumber \end{equation} o que pode ser reagrupado como: \begin{equation} \begin{array}{rcl} v(t+\Delta t) & = & v(t) + a(r,t) \Delta t\\ r(t+\Delta t) & = & r(t) + v(t) \Delta t\\ \end{array} \nonumber \end{equation} Como o método de derivada para frente tem erro $O(h^2)$ então o método de Euler também tem o mesmo erro $O(\Delta t^2)$

Utilizando a nova notação introduzida no calculo de integrais podemos reescrever o método de Euler como \begin{equation} \begin{array}{rcl} v_{i+1} & = & v_i + a_i \Delta t\\ r_{i+1} & = & r_i + v_i \Delta t\\ \end{array} \nonumber \end{equation} De onde podemos ver propriamente o algoritmo, é necessário dar um valor inicial de aceleração, velocidade, posição e o intervalo de tempo ($a_i, v_i, r_i, \Delta t$) e o algoritmo devolve a velocidade e posição deslocado um $\Delta t$ no tempo ($v_{i+1}, r_{i+1}$, como esse novos valores podemos realimentar o algoritmo e calcular a velocidade e posição ao tempo $2\Delta t$ e assim subsecivamente ate atingir o tempo desejado $t = n\Delta t$. Nesse sentido o algoritmo de Euler é um algoritmo iterativo (como todos os que resolvem eq. diferenciáveis)

Método de Euler-Cromer

O método de Euler não funciona muito bem com problemas de tipo oscilatório (ele não conserva a energia, do ponto de vista físico) por tanto devemos modificar o método a fim de que funciones, a modificação é muito sutil, e está descrita a continuação \begin{equation} \begin{array}{rcl} v(t+\Delta t) & = & v(t) + a(r,t) \Delta t\\ r(t+\Delta t) & = & r(t) + v(t+\Delta t) \Delta t\\ \end{array} \nonumber \end{equation} ou na notação nova \begin{equation} \begin{array}{rcl} v_{i+1} & = & v_i + a_i \Delta t\\ r_{i+1} & = & r_i + v_{i+1} \Delta t\\ \end{array} \nonumber \end{equation} a única mudança que foi realizada está no calculo da posição da partícula, em lugar de usar a velocidade ao tempo $t$ foi utilizada a nova velocidade calculada previamente ao tempo $t+\Delta t$

Exemplo 1

É bem conhecido que vários núcleos atômicos são instáveis. Um exemplo típico é o isótopo nuclear $U^{235}U$ (o núcleo de urânio contem 143 nêutrons e 92 prótons, assim um total de 235 núcleons), o qual tem a probabilidade, pequena, mas não insignificante, de decair em dois núcleons de aproximadamente a metade do seu tamanho. Desafortunadamente não é possível prever quando um dado átomo vai decair, o mais que podemos dizer é qual é a probabilidade de que num determinado momento esse átomo decaia. Assim se inicialmente ($t=0$) temos uma amostra com $N_0$ átomos de $U^{235}U$ podemos dizer quanto teremos um tempo $t$ posterior. De fato, se sabe que a equação diferencial que descreve o decaimento de um conjunto de átomos é \begin{equation} \frac{dN}{dt} = -\frac{N}{\tau} \nonumber \end{equation} onde $\tau$ mede de alguma forma a taxa de decaimento que sofre o material que, para o caso do $U^{235}$, está dada por $\tau = 1\times 10^9$ anos.

Nosso objetivo é resolver numericamente a equação diferencial do decaimento radioativo, supondo que temos $100$ átomos inicialmente. utilizando um intervalo de integração de $\Delta t = 0.05s$ e um maximo de tempo de $t = 5s$

Solução

Antes de proceder ao calculo dessa equação diferencial, é bom saber que a solução analítica dessa equação é \begin{equation} N = N_0\exp\,\left(-\frac{N}{\tau}\right) \nonumber \end{equation} onde $N_0$ é a quantidade inicial do átomos e $\tau$ sua vida media.

Para resolver numericamente observamos que podemos utilizar diretamente o método de Euler ou mesmo o método de Euler-Cromer. Usando o Euler, chamemos de $N_0$ o número de átomos iniciais, suponhamos que vamos calcular o número de átomos depois de transcorrido $t_{max}$ anos. Assim o algoritmo deve ser similar a:

ENTRADA: $\Delta t$, $t_{max}$, $\tau$, $N_0$

$ \begin{array}{rcl} n & \leftarrow & t_{max} / \Delta t \\ N_1 & \leftarrow & N_0\\ t_1 & \leftarrow & 0 \end{array} \nonumber $
Para $i = 1, \ldots , (n-1)$ faça
$ \begin{array}{rcl} N_{i+1} & \leftarrow & N_i + \frac{N_i}{\tau} \Delta t \\ t_{i+1} & \leftarrow & t_i + \Delta t \\ \end{array} \nonumber $
fin do faça

SAÍDA: {$(x_1,t_1),\,(x_2,t_2),\, \dots , \, (x_n,t_n)$}
PARE

Exemplo 2

Calcule a trajetória de uma bala de canhão considerando que a resistência do ar está dada por $-Bv^2$.

Solução

Definamos como positivo para acima e direita. Este movimento representa na realidade duas equações diferenciais, uma para $y$ e outra para $x$ \begin{equation} \begin{array}{rcl} \frac{dv_y}{dt} & = & -g - Fa_{y}\\ \frac{dv_x}{dt} & = & -Fa_{x} \end{array} \nonumber \end{equation} onde $Fa$ é a força de arrastro dada por $Fa = -B \sqrt{v_x^2+v_y^2}$. Podemos decompor a força de arrastro da seguinte forma $\vec{F}a = Fa_x \hat{x} + \vec{F}a \hat{y}$

Da figura se observa que o vetor força de arrastro forma o mesmo ângulo $\theta$ que o vetor velocidade com a horizontal, dessa forma é evidente que podemos escrever \begin{equation} \begin{array}{rcl} Fa_x & = & Fa (v_x / v)\\ Fa_y & = & Fa (v_y / v)\\ \end{array} \nonumber \end{equation} ou, igualmente \begin{equation} \begin{array}{rcl} Fa_x & = & -B v_x\\ Fa_y & = & -B v_y\\ \end{array} \nonumber \end{equation} Assim, reescrevemos as equações diferencias que descrevem o movimento como \begin{equation} \frac{dv_y}{dt} = -g - B v_y \nonumber \end{equation} e \begin{equation} \frac{dv_x}{dt} = -B v_x \nonumber \end{equation} de onde podemos aplicar o algoritmo de Euler a cada uma das equações, o que resulta num algoritmo como

ENTRADA: $\Delta t$, $x_0$, $y_0$, $N_0$, $(vx)_0$, $(vy)_0$

$ \begin{array}{lcl} y_1 & \leftarrow & y_0\\ (vy)_1 & \leftarrow & (vy)_0\\ x_1 & \leftarrow & x_0\\ (vx)_1 & \leftarrow & (vx)_0\\ t_1 & \leftarrow & 0\\ \end{array} \nonumber $

faça

$ \begin{array}{rcl} (vx)_{i+1} & \leftarrow & (vx)_i - B( vx )_i \Delta t\\ x_{i+1} & \leftarrow & x_i - (vx)_i \Delta t\\ & & \\ (vy)_{i+1} & \leftarrow & (vy)_i - \left( g - B(vy)_i \right) \Delta t\\ y_{i+1} & \leftarrow & y_i - (vy)_i \Delta t\\ & & \\ t_{i+1} & \leftarrow & t_i + \Delta t \\ \end{array} \nonumber $
Se $y <= 0.0$ Sai do faça
fin do faça
$y \leftarrow 0$
chama Neville$\left(x, y , \{(y_1,x_1),\,(y_2,x_2),\, \ldots , \,(y_n,x_n)\}\right)$
$ \begin{array}{lcl} x_n & \leftarrow & x\\ y_n& \leftarrow & y\\ \end{array} \nonumber $

SAÍDA: {$(x_1,y_1),\,(x_2,y_2),\, \ldots , \, (x_n,y_n)$}
PARE

Exemplo 3

Considere um pendulo de massa $m$ dependurado do teto rigido por uma corda sem massa e inextensível. Se inicialmente o pendulo se solta livremente desde um ponto no qual a corda faz um ângulo de $\theta_0$ radianos descreva o movimento do pendulo.

Solução

A segunda lei de Newton nos diz que a força é igual à massa vezes à aceleração da partículas ao longo do arco circunferência que é a trajetória da partícula, ou seja, $Fe = md^2 s / dt^2$. O deslocamento ao longo deste arco é $s=l\theta$, onde $l$ é o comprimento da corda. Se agora assumimos que $\theta$ é tão pequeno de forma que $\sin \theta \approx \theta$, obtemos a equação (provavelmente familiar) do movimento \begin{equation} \frac{d^2\theta}{dt^2}=-\frac{g}{l}\theta \nonumber \end{equation} Esta é a equação típica de um oscilados harmônico simples. Como sabem a solução desta equação é

sendo $\Omega = \sqrt{g/l}$, e $\theta_0$ e $\phi$ constantes que dependem da velocidade e deslocamento inicial. Ou seja, as oscilações feitas por um pendulo simples são sinusoidais com o tempo e sem mantem por tempo indeterminando nesse estado.

Para resolver este problema numericamente fazemos \begin{equation} \begin{array}{lcl} \frac{d\omega}{dt} & = & -\frac{g}{l}\theta\\ \frac{d\theta}{dt} & = & \omega\\ \end{array} \nonumber \end{equation}

Movimento planetário

Para ver esse problema (extenso) clique aqui

Tarefa

  1. É conhecido que a equação diferencial que descrever o decaimento radioativo de isótopos radioativos, como o $U^{235}$ é descrita pela equação diferencial
    $\frac{dN}{dt}=-\frac{N}{\tau}$
    Onde $\tau$ é a meia-vida do isótopo, $N$ é a quantidade do isótopo ao tempo $t$. No caso específico do $U^{235}$, a meia-vida deste isótopo é de $4,5\times10^9$ anos. Cria um programa que resolva a equação diferencial anterior e compare o efeito da escolha de diversas $\Delta t$
  2. A velocidade de um corpo que cai livremente perto da superfície da Terra é descrito pela equação
    $\frac{dv}{dt}=-g$
    onde $v$ é a velocidade e $g=-9,8m/s^2$ é a aceleração devida à gravidade da Terra. Escreva um programa que empregue método de Euler para calcular a solução da equação anterior, isto é, calcule $v$ em função de $t$. Por simplicidade assuma que a velocidade inicial é zero, isto é, o objeto inicia o movimento partindo do repouso. Repita o calculo para diferentes valores de passo de tempo, compare seus resultados com a solução exata.
  3. Frequentemente se observa que o atrito experimentado por um corpo incrementa com o aumento da velocidade dos corpos sobre o qual atua esta força. Um exemplo disto é o paraquedista. A ideia do paraquedismo é produzir uma força de arrastro considerável maior de aquela que seria produzida se não se utilizasse o paraquedas. Uma primeira descrição deste tipo de dinâmica pode ser obtida a traves da equação diferencial
    $\frac{dv}{dt}=a-bv$
    onde $a$ e $b$ são constantes. Podemos pensar $a$ como tendo origem numa força que se aplica ao corpo, tal como a gravidade, enquanto que $b$ tem origem no atrito experimentado pelo corpo. Note que a força de atrito é negativa (dessa forma assumimos $b>0$), tal que se opõe ao movimento, e aumenta sua magnitude com o aumento da velocidade. Use o método de Euler para resolver essa equação diferencial para $v$ como função de $t$. Uma escolha de parâmetros convenientes é $a=10$ e $b=1$. Poderão perceber que assim que $t\rightarrow \infty$ a velocidade se aproxima de um valor constante conhecido como velocidade terminal.
  4. Problemas de crescimento populacionais frequentemente aparecem como equação diferenciais de primeira ordem. Por exemplo, a equação
    $\frac{dN}{dt}=aN-bN^2$
    deve descrever a variação do número de indivíduos, $N$, como o tempo numa dada população. Aquí o primeiro termo $aN$ corresponde ao nascimento de novos membros, enquanto que o termo $-bN^2$ corresponde aos falecimentos. O termo correspondente ao número de falecimentos escala com $N^2$ porque a medida que a população aumenta, resulta mais difícil encontrar comida. Iniciaremos a resolução da equação diferencial anterior considerando $b=0$. Compare seu resultado com a solução exata. Posteriormente utilize $a=10$ e $b=3$. De uma explicação intuitiva para seus resultados.