Pêndulos Estocásticos

Grupo : Gustavo H. Guesser, Joshua L. Kipper, Marcos Pasa.

Pêndulo Simples

Equação de movimento

Primeiramente vamos inserir ruído em um pêndulo simples, que é constituído de uma barra de comprimento $l$ , sem massa e rígida que contém uma massa $m$ pontual em sua ponta, conforme ilustrado na figura a seguir.

Esquema de um pêndulos simples em um campo gravitacional constante.

Considerando que o pêndulo está sob o efeito da gravidade e se encontra submerso em um fluido viscoso (como o ar), tal que a força de resistência que atua na massa é $-2b{\vec {v}}$ , a equação de movimento é dada por:

${\ddot {\theta }}(t)=-{\frac {2b}{m}}{\dot {\theta }}-{\frac {g}{l}}sen(\theta )$

Vamos supor que existe uma força ruidosa atuando em $m$ ( $F_{r}(t)$ ), que pode ser modelada por um ruído branco gaussiano $\xi (t)$ da seguinte forma

$F_{r}(t)=m\alpha \xi (t)$

em que $\alpha$ é a intensidade do ruído. $\xi (t)$ é caracterizado pelas seguintes propriedades:

$\langle \xi (t)\rangle =0,~~\langle \xi (t_{2})\xi (t_{1})\rangle =\delta (t_{2}-t_{1})$

Adicionando essa nova força nas equações de movimento, ficamos com

${\ddot {\theta }}(t)=-{\frac {2b}{m}}{\dot {\theta }}-{\frac {g}{l}}sen(\theta )+{\frac {\alpha }{l}}\xi (t)$

A partir de agora, por questão de simplicidade, vamos supor que $g=l=1$ , então

${\ddot {\theta }}(t)=-2b{\dot {\theta }}-sen(\theta )+\alpha \xi (t)$

Método de integração

Vamos montar um métodos para integrar o sistema no tempo. Primeiramente vamos dividir a equação em duas equações diferencias de primeira ordem, introduzindo a variável $\omega :={\dot {\theta }}$ , então ficamos com o seguinte sistema

${\begin{aligned}{\dot {\theta }}&=\omega \\{\dot {\omega }}&=-2b\omega -sen(\theta )+\alpha \xi (t)\end{aligned}}$

que pode ser escrito na forma diferencial

${\begin{aligned}d\theta &=\omega dt\\d\omega &=(-2b\omega -sen(\theta ))dt+\alpha \xi (t)dt\end{aligned}}$

mas $\xi (t)dt$ é o incremento do processo de Wiener ( $W(t)=\int _{0}^{t}\xi (t')dt'$ ), então

${\begin{aligned}d\theta &=\omega dt\\d\omega &=(-2b\omega -sen(\theta ))dt+\alpha dW(t)\end{aligned}}$

Discretizando o tempo e lembrando que a densidade de probabilidade de transição de $W(t)$ para $W(t+\Delta t)$ tem desvio padrão igual a ${\sqrt {\Delta t}}$

${\begin{aligned}\theta _{j+1}&=\theta _{j}+\omega _{j}\Delta t\\\omega _{j+1}&=\omega _{j}+(-2b\omega _{j}-sen(\theta _{j}))\Delta t+\alpha {R_{G}}_{j}{\sqrt {\Delta t}}\end{aligned}}$

em que $R_{G}$ é uma amostra de uma distribuição gaussiana com média 0 e variância 1, e o método de Euler foi utilizado para a parte determinística da equação.

Nas próximas seções será analisado a energia do sistema, e como o método de Euler não é muito bom para preservar a energia de sistemas conservativos, será utilizado o método preditor corretor (com adição do método de Heun para $\theta$ ) para a parte determinística da equação, que consiste nos seguintes passos:

Calcular um theta intermediário:

    $\theta _{j+1}^{(2)}=\theta _{j}+\omega _{j}\Delta t$

Com $\theta _{j+1}^{(2)}$ calcular um theta médio e utilizá-lo para obter um omega intermediário:

    ${\begin{aligned}{\bar {\theta }}_{j}&={\frac {\theta _{j+1}^{(2)}+\theta _{j}}{2}}\\\omega _{j+1}^{(2)}&=\omega _{j}+f({\bar {\theta }}_{j},\omega _{j},{R_{G}}_{j})\end{aligned}}$

Em que

f

é a expressão do método de Euler visto logo acima.

Recalcular theta utilizando um omega intermediário

    ${\begin{aligned}{\bar {\omega }}_{j}&={\frac {\omega _{j+1}^{(2)}+\omega _{j}}{2}}\\\theta _{j+1}&=\theta _{j}+{\bar {\omega }}_{j}\Delta t\end{aligned}}$

Recalcular omega com um theta intermediário atualizado

    ${\begin{aligned}{\bar {\theta }}_{j}^{(2)}&={\frac {\theta _{j+1}+\theta _{j}}{2}}\\\omega _{j+1}&=\omega _{j}+f({\bar {\theta }}_{j}^{(2)},{\bar {\omega }}_{j},{R_{G}}_{j})\end{aligned}}$

OBS: No cálculo de

\omega _{j+1}^{(2)}

e

\omega _{j+1}

foi utilizado o mesmo

{R_{G}}_{j}

.

Energia (Sem amortecimento)

Logo após terminar a implementação do método numérico, rapidamente notamos que a adição do ruído aumento a energia mecânica do pêndulo ( $E$ ), vamos explorar esse fenômeno, sempre utilizando $\beta =0$ . Para ilustrar esse efeito, segue uma animação do pêndulo partindo do repouso na configuração de equilíbrio estável ( $\theta (0)=0,~\omega (0)=0$ ) com $\alpha =0.1$

Pêndulo partindo do repouso com ruído.

Para realizar uma exploração quantitativa, o seguinte procedimento foi feito para vários valores de $\alpha$ :

Utilizando $\Delta t=0.01$ , integrar o sistema até $t_{f}=100$ , calculando e armazenando a energia em cada passo temporal. Repetir essa integração 700 vezes para fazer médias temporais. Como exemplo, segue os dados obtidos da energia em função do tempo (com a média temporal feita) para um determinado $\alpha$ utilizado

Energia média em função tempo.

O gráfico nos indica que o ruído gera uma potência média sobre o pêndulo.

Realizar um ajuste linear nos dados $\langle E\rangle _{t}\times t$ para obter o coeficiente angular, que corresponde a potência média gerada pelo ruído ( ${\bar {P}}$ ).

Com as simulações executadas, foi realizado o gráfico ${\bar {P}}\times \alpha$ . Notamos que os dados se alinham em linha reta com os eixos em escala logarítmica, ou seja, os mesmos seguem uma lei de potência ${\bar {P}}=ae^{b\alpha }$ ), então foi realizado outro ajuste linear para encontra o expoente $b$ , a figura a seguir ilustra os dados e os resultados do ajuste: