Lançamento Oblíquo Estocástico

O lançamento oblíquo é um clásico problema da mecanica, onde um projétil e lançado com velocidade inicial $v$ em uma direção que faz um angulo $\theta$ com a horizontal. Nesse casso iremos considerar o movimento com arrasto em que a força de resistencia do ar é oposta ao movimento e proporcional a velocidade

$F_{Ar}=-kmv$

O objetivo é descrever o movimento do projétil em duas dimensões, levando em consideração tanto a força da gravidade quanto a resistência do ar então adicionar um termo de ruido no sistema e analizar seu comportamento.

Equações de Movimento

Começamos decompondo as velocidades nas direções x e y

$v_{x}=v\cos(\theta )$

$v_{y}=v\sin(\theta )$

Para as equações de movimento usamos a segunda lei de Newton nas componentes horizontais e verticais

$m{\frac {d^{2}x}{dt^{2}}}=-kmv_{x}$

$m{\frac {d^{2}y}{dt^{2}}}=-mg-kmv_{y}$

onde:

• $m$ é a massa do projétil

• $g$ é a aceleração da gravidade

• $k$ é o coeficiente de arrasto

Introduzindo $v_{x}={\frac {dx}{dt}}$ e $v_{y}={\frac {dy}{dt}}$ podemos reescrever as equações como diferenciais de primeira ordem.

${\frac {dv_{x}}{dt}}=-kv_{x}$

${\frac {dv_{y}}{dt}}=-g-kv_{y}$

Essas E.D. possuem solução analitica bastando integrar nos limites adequados. A integração numérica por Euler é dirta a partir das relações a cima e nos da uma trajetoria para comparação quando adicionar o termo de ruido:

$v_{x,j+1}=v_{x,j}-kv_{x,j}dt$

$v_{y,j+1}=v_{y,j}-gdt-kv_{y,j}dt$

Utilizando as seguintes condições iniciais $x_{0}=0$ e $y_{0}=0$ podemos atualizar a posição

$v_{x}={\frac {dx}{dt}}$

$x_{j+1}=x_{j}+v_{x,j}dt$

$v_{y}={\frac {dy}{dt}}$

$y_{j+1}=y_{j}+v_{y,j}dt$

Para modelar o ruído produzido pela densidade do ar, consideramos dois regimes. No primeiro, assumimos que a densidade é constante ao longo da trajetória, ou seja, $\rho =\rho _{0}$ e, portanto, trataremos o problema como uma EDE com ruído aditivo, onde $B_{(X(t),t)}=\rho _{0}=\beta$. Na segunda situação, consideramos que a densidade varia com a altitude do projétil, seguindo a relação $\rho =\rho _{0}e^{-{\frac {y}{H}}}$ e, nesse caso, abordaremos como uma EDE com ruído multiplicativo, onde $B_{(X(t),t)}=\rho _{(X(t))}=\beta e^{-{\frac {y}{H}}}$.

Equação diferencial estocástica

A equação diferencial estocástica a ser resolvida é:

${\frac {dX}{dt}}=A(X(t))+B(X(t))\xi (t),$

onde

$X={\begin{pmatrix}x\\y\\v_{x}\\v_{y}\end{pmatrix}},\quad A={\begin{pmatrix}v_{x}\\v_{y}\\-kv_{x}\\-g-kv_{y}\end{pmatrix}},\quad B={\begin{pmatrix}0&0&0&0\\0&0&0&0\\0&0&\beta e^{-{\frac {y}{H}}}&0\\0&0&0&\beta e^{-{\frac {y}{H}}}\end{pmatrix}},\quad \xi ={\begin{pmatrix}0\\0\\\xi _{x}(t)\\\xi _{y}(t)\end{pmatrix}}$

Caso o termo $B(X,t)$ fosse uma constante, exemplo em que a densidade do ar não muda significativamente, o ruído é dito aditivo. Então, a integração da equação é direta, resultando em:

$dX(t)=A(X(t))dt+BdW(t),$

onde $\xi (t)dt=dW(t)$ é o incremento de Wiener, um processo estocástico com distribuição gaussiana de largura ${\sqrt {dt}}$. Escrevendo a equação em diferenças finitas e substituindo os valores de $A(X(t))$ e de $B$ têm-se:

$\Delta x_{i}=v_{x_{i}}\Delta t$

$\Delta v_{x_{i}}=-kv_{x_{i}}\Delta t+\beta \Delta W_{x_{i}}$

$\Delta y_{i}=v_{y_{i}}\Delta t$

$\Delta v_{y_{i}}=\left(-g-kv_{x_{i}}\right)\Delta t+\beta \Delta W_{x_{i}}$

Entretanto, o termo $B(X,t)$ depende do processo estocástico e o ruído é dito multiplicativo, assim, devemos utilizar um cálculo estocástco especial, por exemplo o cálculo de Itô e de Stratonovich. No cálculo de Stratonovich, as regras de integração, diferenciação, mudança de variáveis etc, se mantém as mesmas do cálculo usual, e a equação é escrita como:

$dX(t)=A(X(t))dt+B(X(t))dW(t),$

com o detalhe que o argumento de $B(X(t))$ é, na realidade, calculado na média de X entre os limites do intervalo de integração. Ou seja, ao escrever a equação em diferenças finitas, obtém-se:

$\Delta X(t)=A(X(t))\Delta t+B\left(X(t)+{\frac {1}{2}}\Delta X(t)\right)\Delta W(t).$

Agora, é necessário explicitar o $\Delta X(t)$ na equação acima, porém a dependência de $B(X(t))$ não é linear em relação à altura, o que torna inviável esse procedimento.

Com isso, a atenção volta-se para o cálculo de Itô, dado por:

$dX(t)=\left[A(X(t))+{\frac {1}{2}}{\frac {\partial B}{\partial X}}B(X(t))\right]dt+B(X(t))\bullet dW(t)$

No caso em que $B$ não depende de $X(t)$, a derivada que aparece na equação é zero e é recuperado o caso de ruído aditivo já mencionado. Aqui, vale ressaltar que $B(X(t))$ é uma matriz e está sendo derivada com relação ao vetor $X(t)$, então o cálculo para as componentes de $X_{i}(t)$ são:

$dX_{i}(t)=A_{i}^{(I)}(X(t))dt+\sum _{\alpha =1}^{m}B_{i,\alpha }(X(t))\bullet dW_{\alpha }(t),$

onde

$A_{i}^{(I)}(X(t))=A_{i}(X(t))+{\frac {1}{2}}\sum _{j=1}^{n}\sum _{\alpha =1}^{m}{\frac {\partial B_{i,\alpha }}{\partial X_{j}}}B_{j,\alpha }(X(t)).$

Após a substituição dos valores de $A_{i}(X(t))$ e $B_{i}(X(t))$, obtêm-se as seguintes expressões:

$\Delta x_{i}=v_{x_{i}}\Delta t$

$\Delta v_{x_{i}}=\left[-kv_{x_{i}}-{\frac {\beta ^{2}}{2H}}e^{-{\frac {2}{H}}y_{i}}\right]\Delta t+\beta e^{-{\frac {y(t)}{H}}}\Delta W_{x_{i}}$

$\Delta y_{i}=v_{y_{i}}\Delta t$

$\Delta v_{y_{i}}=\left[-g-kv_{x_{i}}-{\frac {\beta ^{2}}{2H}}e^{-{\frac {2}{H}}y_{i}}\right]\Delta t+\beta e^{-{\frac {y_{i}}{H}}}\Delta W_{y_{i}}$

É importante ressaltar que os incrementos de Wiener que aparecem nas equações para velocidade são independentes entre si, isto é, deverá ser gerado duas variáveis aleatórias com distribuição gaussiana de largura ${\sqrt {\Delta t}}$. Além disso, um detalhe interessante é que a equação para $\Delta v_{x}$ é semelhante à equação para $\Delta v_{y}$ em que $g=0$.

Resultados

Nas figuras abaixo podemos ver as trajetórias de certo objeto para ruído aditivo e ruído multiplicativo considerando o chão em zero.

Esperávamos que as média para ruído aditivo e multiplicativos fossem iguais, dado que o ruído não interfere no cálculo da média pois a probabilidade de puxar em qualquer direção é a mesma, porém como pode ser visto no gráfico abaixo há divergência para o ruído multiplicativo.

Tanto a divergência da médias como o aumento da dispersão da trajetórias podem ser explicados pelo modelo que estamos usando na EDE-ito $(\beta e^{-{\frac {y_{i}}{H}}})$, nesse modelo conforme temos y variando para valores negativos a exponencial explode consequentemente temos aumento na intensidade do ruído, esse aumento pode ser observado no gráfico abaixo.

Temos que para poucos passos o ganho do ruído em qual direção causara grande divergência, ou seja, para termos vermos a convergência para média do ruído multiplicativo precisamos incrementar o número de trajetórias, a figura a seguir foi feita para um aumento de dez vezes no número de passos comparado ao da figura anterior.

Nos gráficos abaixo temos as trajetórias para o ruído multiplicativo, as médias para os dois métodos e a solução por euler levando em conta o solo, para as trajetória é visível a dispersão em torno da solução de euler que acontece por causa do ruído, porém como comentado anteriormente essa dispersão não influência na trajetória média, da mesma forma que as duas primeiras figuras conforme a diminuição das velocidades pelo termo dissipativo temos maior influência do ruído no trajeto que pode até mesmo ocasionar em um retorno do objeto na direção da posição inicial porém a média continua somente para a direita.

Agora para um maior valor de ruído podemos observar o gráfico a seguir, onde notamos que a média, na região de contato com o solo, fica acima da solução de euler, isso ocorre pela condição de contorno que estamos utilizando, quando a posição da partícula chega em zero invertemos o sinal da velocidade mas ao fazermos isso alteramos também o ruído que posteriormente faria a partícula tender para baixo porém agora, com a troca de sinal, estará contribuído para cima assim temos uma deslocamento da média para cima. Neste caso temos que o deslocamento para o aditivo é maior do que o do multiplicativo logo o ruído do aditivo é mais intenso.