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êndulo 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\overrightarrow{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:

$⟨\xi (t)⟩=0,⟨\xi (t_2)\xi (t_1)⟩=\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{array}{rl}\dot{\theta }& =\omega \\ \dot{\omega }& =-2b\omega -sen(\theta )+\alpha \xi (t)\end{array}$

que pode ser escrito na forma diferencial

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

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

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

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

$\begin{array}{rl}{\theta }_{j+1}& ={\theta }_j+{\omega }_j\mathrm{\Delta }t\\ {\omega }_{j+1}& ={\omega }_j+(-2b{\omega }_j-sen({\theta }_j))\mathrm{\Delta }t+\alpha {R_G}_j\sqrt{\mathrm{\Delta }t}\end{array}$

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\mathrm{\Delta }t$

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

   $\begin{array}{rl}{\overline{\theta }}_j& =\frac{{\theta }_{j+1}^{(2)}+{\theta }_j}{2}\\ {\omega }_{j+1}^{(2)}& ={\omega }_j+f({\overline{\theta }}_j,{\omega }_j,{R_G}_j)\end{array}$

Em que $f$ é a expressão do método de Euler visto logo acima.

• Recalcular theta utilizando um omega intermediário

   $\begin{array}{rl}{\overline{\omega }}_j& =\frac{{\omega }_{j+1}^{(2)}+{\omega }_j}{2}\\ {\theta }_{j+1}& ={\theta }_j+{\overline{\omega }}_j\mathrm{\Delta }t\end{array}$

• Recalcular omega com um theta intermediário atualizado

   $\begin{array}{rl}{\overline{\theta }}_j^{(2)}& =\frac{{\theta }_{j+1}+{\theta }_j}{2}\\ {\omega }_{j+1}& ={\omega }_j+f({\overline{\theta }}_j^{(2)},{\overline{\omega }}_j,{R_G}_j)\end{array}$

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 $\mathrm{\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 $⟨E{⟩}_t\times t$ para obter o coeficiente angular, que corresponde a potência média gerada pelo ruído ($\overline{P}$).

Com as simulações executadas, foi realizado o gráfico $\overline{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 $\overline{P}=a{e}^{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:

Potência em função do ruído ($\alpha$). O painel da esquerda possui eixos em escala linear e o da direita em escala logarítmica.

Portanto, $\overline{P}$ aumenta, aproximadamente, de forma quadrática com $\alpha$. Por fim, gostaríamos de mostrar que mesmo para $\alpha$ muito pequeno, ainda existe energia sendo injetada no sistema com taxa constante, e isso não é um artefato da simulação. Para tal, foram realizadas 700 simulações com $\alpha =5\cdot 10^{-4}$ e os resultados foram comparados com $\alpha =0$

Energia média para $\alpha$ muito pequeno comparado com $\alpha$ nulo.

Energia (Com amortecimento)

Até o momento, o amortecimento foi negligenciado. Vamos, então, introduzi-lo e rodar várias simulações (700, neste caso) e observar como a média temporal da energia evolui ao longo do tempo, assim como foi feita na seção anterior. Utilizando $\beta =0.1$ e $\alpha =0.485$ foi obtido o seguinte resultado

Energia mecânica média em função do tempo com amortecimento.

claramente o comportamento neste caso é diferente do observado sem amortecimento, agora a anergia aumenta até um certo valor e permanecenele. Para explorar este novo fenômeno, os seguintes passos foram feitos para cada valor de $\beta =0.01,0.1,0.2$:

• Para diversos valores de $\alpha$, executar 700 simulações até a energia estabilizar, salvando a média da energia entre as simulações.
• Para cada conjunto de dados gerados por um determinado $\alpha$, selecionar um intervalo de tempo onde a energia está estabilizada e calcular a sua média ($⟨E⟩$).

Produzindo o gráfico de $⟨E⟩\times \alpha$ obtemos

Energia estabilizada média em função de $\alpha$ para diferentes valores de amortecimento. Os painéis da esquerda possui eixos em escala linear e os da direita em escala logarítmica.

as linhas vermelhas são os melhores ajustas de leis de potência na forma $⟨E⟩=a{\alpha }^b$. Para $\beta =0.1$ os dados utilizados no ajuste foram apenas até $\alpha =0.6$ (indicado pela reta preta vertical no gráfico), pois após esse limite, a lei de potência deixa de ser um ótimo ajuste. É chamativo o fato de todos os coeficientes, independente de $\beta$, serem aproximadamente 2.

Pêndulo invertido

O próximo pêndulo a ser considerado é um pêndulo invertido, que possui um potencial harmônico em seu ângulo, ou seja, um dos termos de sua energia potential é $\frac{k}{2}{\theta }^2$, sendo que agora $\theta$ é zero quando a haste está apontando para cima, conforme ilustrado na figura a seguir

Esquema do pêndulo invertido com movimento vertical livre na base.

Ainda, a base do pêndulo é livre para movimentar-se na direção vertical, é justamente nesse local onde será adicionado uma força ruidosa. A equação de movimento neste caso é

$\ddot{\theta }=\frac{-2b}{m}\dot{\theta }+\frac{1}{l}(g+\ddot{u})sen(\theta )-\frac{k}{m}\theta$

O primeiro termo vem da resistência do ar, o segunda se origina da gravidade e do deslocamento de $u$ e o último provém da "mola" em $\theta$. Vamos supor que existe um vínculo mantendo a base fixa, mas o local onda o pêndulo é fixado pode ser movimentar de forma aleatória na direção vertical, suponde que tal movimento pode ser modelado por ruído branco gaussiano, segue que $\ddot{u}(t)=\alpha \xi (t)$. Introduzindo a variável $\omega :=\dot{\theta }$, ficamos com o seguinte sistemas de equações na forma diferencial

$\begin{array}{rl}d\theta & =\omega dt\\ d\omega & =(\frac{-2b}{m}\omega +\frac{g}{l}sen(\theta )-\frac{k}{m}\theta )dt+sen(\theta )\frac{\alpha }{l}\stackrel{\xi (t)dt}{\overbrace{dW(t)}}\end{array}$

Note que agora o ruído é multiplicativo, em contraste com o ruído aditivo dos pêndulos anteriores, para lidar com esta complicação, no momento da integração vamos utilizar um $\theta$ médio no argumento do seno que multiplica $dW$

$\begin{array}{rl}{\theta }_{j+1}& ={\theta }_j+{\omega }_j\mathrm{\Delta }t\\ {\overline{\theta }}_j& =({\theta }_{j+1}+{\theta }_j)/2\\ {\omega }_{j+1}& ={\omega }_j+(\frac{-2b}{m}{\omega }_j+\frac{g}{l}sen({\overline{\theta }}_j)-\frac{k}{m}{\overline{\theta }}_j)\mathrm{\Delta }t+sen({\overline{\theta }}_j)\frac{\alpha }{l}{R_G}_j\sqrt{\mathrm{\Delta }t}\end{array}$

Esse tipo de pêndulo é de grande interesse em algumas áreas, como na engenharia estrutural, pois uma coluna pode ser modelada como um pêndulo invertido com uma mola no ângulo, o ruído da base pode representar um terremoto.