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{\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$, tal que sua componente tangencial ($F_{r}(t)$) 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étodo 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 gera um aumento na 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:

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, ${\bar {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 permanece nele. 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 ($\langle E\rangle$).

Produzindo o gráfico de $\langle E\rangle \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 $\langle E\rangle =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 segundo 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 onde o pêndulo é fixado pode se movimentar de forma aleatória na direção vertical, supondo 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 os seguintes sistemas de equações na forma diferencial

${\begin{aligned}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}}\overbrace {dW(t)} ^{\xi (t)dt}\end{aligned}}$

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{aligned}\theta _{j+1}&=\theta _{j}+\omega _{j}\Delta t\\{\bar {\theta }}_{j}&=(\theta _{j+1}+\theta _{j})/2\\\omega _{j+1}&=\omega _{j}+({\frac {-2b}{m}}\omega _{j}+{\frac {g}{l}}sen({\bar {\theta }}_{j})-{\frac {k}{m}}{\bar {\theta }}_{j})\Delta t+sen({\bar {\theta }}_{j}){\frac {\alpha }{l}}{R_{G}}_{j}{\sqrt {\Delta t}}\end{aligned}}$

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

Uma questão interessante neste modelo é sobre a estabilidade da configuração $\theta _{0}=\omega _{0}=0$ quando é adicionado ruído e amortecimento. Se $k$ é grande o suficiente, sem ruído, o equilíbrio é estável, mas é possível sair dessa configuração adicionando ruído, cujo valor limite vai depender do amortecimento. Podemos demostrar isso integrando o sistema, partindo da condição inicial $\theta _{0}=10^{-4},\omega _{0}=0$, com valores de $\alpha$ muito próximos:

os seguintes valores foram utilzados

• g = l = 1
• k = 1.1

No gráfico da esquerda, o ângulo oscilou um pouco e permaneceu em 0, já no outro, eventualmente, o ângulo explodiu. Abaixo segue uma animação dessa situação, mas com $\alpha =1$ para o pêndulo sair do equilíbrio mais rápido

Pêndulo invertido partindo do equilíbrio com ruído.

Pêndulo Duplo

O pêndulo duplo estocástico é um sistema dinâmico que combina a complexidade intrínseca do pêndulo duplo com a introdução de elementos de aleatoriedade ou incerteza, tornando o comportamento do sistema ainda mais imprevisível e caótico. O pêndulo duplo em si é um exemplo clássico de um sistema caótico, onde pequenas variações nas condições iniciais podem resultar em trajetórias drasticamente diferentes. Quando um termo estocástico é adicionado, por exemplo, na forma de uma força externa aleatória ou de flutuações nos parâmetros do sistema, a análise e a previsão do movimento se tornam desafiadoras.

Pêndulo Duplo não estocástico.

Equação de movimento

O pêndulo duplo consiste em dois pêndulos acoplados, onde o segundo pêndulo está suspenso na extremidade do primeiro. As equações de movimento para o pêndulo duplo sem termos estocásticos podem ser derivadas utilizando as equações de Lagrange, considerando as coordenadas angulares $\theta _{1}$ e $\theta _{2}$ como as variáveis generalizadas e são dadas por:

${\ddot {\theta }}_{1}={\frac {-g(2m_{1}+m_{2})\sin {\theta _{1}}-m_{2}g(\theta _{1}-2\theta _{2})-2\sin {(\theta _{1}-\theta _{2})}m_{2}({\dot {\theta _{2}}}^{2}l_{2}+{\dot {\theta _{1}}}^{2}l_{1}\cos {(\theta _{1}-\theta _{2})})}{l_{1}(2m_{1}+m_{2}-m_{2}\cos {(2\theta _{1}-2\theta _{2})})}}$

${\ddot {\theta }}_{2}={\frac {2\sin {(\theta _{1}-\theta _{2})}({\dot {\theta _{1}}}^{2}l_{1}(m_{1}+m_{2})+g(m_{1}+m_{2})\cos {\theta _{1}}+{\dot {\theta _{2}}}^{2}l_{2}m_{2}\cos {(\theta _{1}-\theta _{2})})}{l_{2}(2m_{1}+m_{2}-m_{2}\cos {(2\theta _{1}-2\theta _{2})})}}$

Mantendo $m_{1}$, $m_{2}$, $l$, e $g$ fixos, defina :

$f_{1}(\theta _{1},{\dot {\theta _{1}}},\theta _{2},{\dot {\theta _{1}}}):={\frac {-g(2m_{1}+m_{2})\sin {\theta _{1}}-m_{2}g(\theta _{1}-2\theta _{2})-2\sin {(\theta _{1}-\theta _{2})}m_{2}({\dot {\theta _{2}}}^{2}l_{2}+{\dot {\theta _{1}}}^{2}l_{1}\cos {(\theta _{1}-\theta _{2})})}{l_{1}(2m_{1}+m_{2}-m_{2}\cos {(2\theta _{1}-2\theta _{2})})}}$

$f_{2}(\theta _{1},{\dot {\theta _{1}}},\theta _{2},{\dot {\theta _{1}}}):={\frac {2\sin {(\theta _{1}-\theta _{2})}({\dot {\theta _{1}}}^{2}l_{1}(m_{1}+m_{2})+g(m_{1}+m_{2})\cos {\theta _{1}}+{\dot {\theta _{2}}}^{2}l_{2}m_{2}\cos {(\theta _{1}-\theta _{2})})}{l_{2}(2m_{1}+m_{2}-m_{2}\cos {(2\theta _{1}-2\theta _{2})})}}$

assim podemos escrever as equações dinâmicas de maneira mais compacta,

${\ddot {\theta }}_{1}=f_{1}(\theta _{1},{\dot {\theta _{1}}},\theta _{2},{\dot {\theta _{1}}})$

${\ddot {\theta }}_{2}=f_{2}(\theta _{1},{\dot {\theta _{1}}},\theta _{2},{\dot {\theta _{1}}})$

Suponha que haja um força externa ruidosa agindo sobre $m_{1}$ tal que, $F_{1}^{r}=\sigma _{1}\xi _{1}(t)$, onde $\sigma _{1}$ além de representar a intensidade do ruído gaussiano, terá consigo todas as constantes agrupadas, suponha o mesmo para $m_{1}$. Assim as equações de movimento se tornam:

${\ddot {\theta }}_{1}=f_{1}(\theta _{1},{\dot {\theta _{1}}},\theta _{2},{\dot {\theta _{1}}})+\sigma _{1}\xi _{1}(t)$

${\ddot {\theta }}_{2}=f_{2}(\theta _{1},{\dot {\theta _{1}}},\theta _{2},{\dot {\theta _{1}}})+\sigma _{2}\xi _{2}(t)$

Portanto, obtemos um sistema de equações diferencias estocásticas com ruído gaussiano aditivo.

Método de integração

Para a integração numérica desse sistema utilizaremos o método Runge-Kutta 4. Seja $\omega _{i}={\dot {\theta _{i}}}$, reduzindo a ordem das equações diferencias, o sistema se torna:

${\begin{aligned}\omega _{1}&={\dot {\theta _{1}}}\\{\dot {\omega }}_{1}&=f_{1}(\theta _{1},\omega _{1},\theta _{2},\omega _{2})+\sigma _{1}\xi _{1}(t)\\\omega _{2}&={\dot {\theta _{2}}}\\{\dot {\omega }}_{2}&=f_{2}(\theta _{1},\omega _{1},\theta _{2},\omega _{2})+\sigma _{2}\xi _{2}(t)\end{aligned}}$

na forma diferencial:

${\begin{aligned}\omega _{1}dt&=d{\theta _{1}}\\d\omega _{1}&=f_{1}(\theta _{1},\omega _{1},\theta _{2},\omega _{2})dt+\sigma _{1}dW_{1}(t)\\\omega _{2}dt&=d{\theta _{2}}\\d\omega _{2}&=f_{2}(\theta _{1},\omega _{1},\theta _{2},\omega _{2})dt+\sigma _{2}dW_{2}(t)\end{aligned}}$

onde $dW_{i}(t)=\xi _{i}(t)dt$ é o incremento do processo de Wiener. Discretizando as equações diferencias e aplicando o método RK4, obtemos:

${\begin{aligned}k_{1}^{\theta _{1}}&=dt\cdot \omega _{1}^{n}\\k_{1}^{\omega _{1}}&=dt\cdot f_{1}(\theta _{1}^{n},\omega _{1}^{n},\theta _{2}^{n},\omega _{2}^{n})\\k_{2}^{\theta _{1}}&=dt\cdot (\omega _{1}^{n}+k_{1}^{\omega _{1}}/2)\\k_{2}^{\omega _{1}}&=dt\cdot f_{1}(\theta _{1}^{n}+k_{1}^{\theta _{1}}/2,\omega _{1}^{n}+k_{1}^{\omega _{1}}/2,\theta _{2}^{n},\omega _{2}^{n})\\k_{3}^{\theta _{1}}&=dt\cdot (\omega _{1}^{n}+k_{2}^{\omega _{1}}/2)\\k_{3}^{\omega _{1}}&=dt\cdot f_{1}(\theta _{1}^{n}+k_{2}^{\theta _{1}}/2,\omega _{1}^{n}+k_{2}^{\omega _{1}}/2,\theta _{2}^{n},\omega _{2}^{n})\\k_{4}^{\theta _{1}}&=dt\cdot (\omega _{1}^{n}+k_{3}^{\omega _{1}})\\k_{4}^{\omega _{1}}&=dt\cdot f_{1}(\theta _{1}^{n}+k_{3}^{\theta _{1}},\omega _{1}^{n}+k_{3}^{\omega _{1}},\theta _{2}^{n},\omega _{2}^{n})\\\theta _{1}^{n+1}&=\theta _{1}^{n}+(k_{1}^{\theta _{1}}+2k_{2}^{\theta _{1}}+2k_{3}^{\theta _{1}}+k_{4}^{\theta _{1}})/6\\\omega _{1}^{n+1}&=\omega _{1}^{n}+(k_{1}^{\omega _{1}}+2k_{2}^{\omega _{1}}+2k_{3}^{\omega _{1}}+k_{4}^{\omega _{1}})/6+\sigma _{1}R_{Gn}{\sqrt {dt}}\end{aligned}}$

novamente, $R_{G}$ é uma amostra de distribuição gaussiana de média 0 e variância 1. É completamente análogo para $\omega _{2}$ e $\theta _{2}$. Utilizamos o RK4 pois necessitávamos de um alta precisão no cálculos por se tratar de um sistema caótico, simulações feitas sem o ruído externo mostram oscilações na energia total da ordem de $10^{-4}$, o que se espera de um método de quarta ordem.

Retrato de fase

O retrato de fase é uma ferramenta gráfica que representa a evolução do sistema dinâmico em termos de variáveis de estado, geralmente o ângulo $\theta$ e a velocidade angular $\omega$ no caso de sistemas oscilatórios como pêndulos. Este tipo de gráfico permite uma visualização clara de como o sistema se comporta ao longo do tempo, identificando padrões, ciclos limitados, ou comportamentos caóticos. Quando o sistema é perturbado por ruídos, o retrato de fase mostrará como essas perturbações afetam a trajetória do sistema, introduzindo variações que podem transformar oscilações regulares em comportamentos mais complexos ou caóticos.

Caso 1 : Ruído em θ₁

Pêndulo duplo estocástico com ruído branco em $\theta _{1}$, $\sigma =0.5$ e $dt=0.01$.

Neste cenário, o ruído afeta apenas o primeiro pêndulo, introduzindo variações na trajetória de $\theta _{1}$ e $\omega _{1}$. O retrato de fase para $\theta _{1}$ versus $\omega _{1}$ mostra uma dispersão maior em comparação com $\theta _{2}$ versus $\omega _{2}$, onde o comportamento pode permanecer mais regular se não houver acoplamento forte entre os pêndulos.

Caso 2 : Ruído em θ₂

Pêndulo duplo estocástico com ruído branco em $\theta _{2}$, $\sigma =0.5$ e $dt=0.01$.

Com ruído aplicado apenas ao segundo pêndulo, $\theta _{2}$ e $\omega _{2}$ exibem uma dispersão significativa no retrato de fase, enquanto $\theta _{1}$ versus $\omega _{1}$ pode permanecer mais regular. A interação entre os pêndulos, no entanto, pode ainda transferir alguma irregularidade de $\theta _{2}$ para $\theta _{1}$.

Caso 3 : Ruídos Iguais em θ₁ e θ₂

Pêndulo duplo estocástico com ruído branco igual em $\theta _{1}$ e $\theta _{2}$, $\sigma =0.5$ e $dt=0.01$.

Quando o mesmo ruído é aplicado a ambos os pêndulos, as trajetórias de $\theta _{1}$ e $\theta _{2}$ serão correlacionadas, resultando em retratos de fase que se comportam de maneira similar. Pode-se observar um padrão sincrônico onde ambos os pêndulos respondem de maneira similar às perturbações, mantendo uma relação quase constante entre $\theta _{1}$ e $\theta _{2}$.

Caso 4 : Ruídos Diferentes em θ₁ e θ₂

Pêndulo duplo estocástico com ruídos brancos diferentes em $\theta _{1}$ e $\theta _{2}$, $\sigma =0.5$ e $dt=0.01$.

Aplicar ruídos diferentes a cada pêndulo pode resultar no comportamento mais complexo e caótico. As trajetórias de $\theta _{1}$ versus $\omega _{1}$ e $\theta _{2}$ versus $\omega _{2}$ podem mostrar padrões bastante distintos, refletindo a falta de correlação entre as perturbações. Esse cenário pode conduzir a uma maior imprevisibilidade no comportamento do sistema como um todo, especialmente em regiões do espaço de fase onde os pêndulos interagem fortemente.

Energia

Afim de compreender como se comporta a energia mecânica do pêndulo duplo com o acréscimo de uma força externa ruidosa, vamos analisar 4 sistemas de pêndulo duplo onde o ruído é adicionado de alguma maneira específica. Alguns pontos a serem observados são :

• Em todas as simulações fixamos $m_{1}$ = $m_{2}$ = $l_{1}$ = $l_{2}$ = $1$ ;
• As condições inicias serão sempre $\theta _{1}(t=0)$ = $\theta _{2}(t=0)$ = $\omega _{1}(t=0)$ = $\omega _{2}(t=0)$ = $0$ ;
• Nos casos onde são calculadas grandezas médias, foram realizadas 100 simulações.