Exibir código-fonte para Equação de Fokker-Planck

'''Grupo: Álison Soares, Rodrigo Avancini Lara e Samuel Huff Dieterich'''

A equação de Fokker-Planck foi aplicada em primeiro modo em problemas relacionados ao movimento Browniano, como veremos à seguir. Nesse caso, lidando com flutuações originadas de vários pequenos distúrbios, as partículas de interesse se chocavam com as moléculas do meio, provocando uma trajetória imprevisível. Por conta dessas flutuações, é impossível determinar a posição exata dessas partículas. Porém, é possível determinar a probabilidade de encontrá-las em determinada região.

Esta equação pode ser obtida a partir da equação de Langevin e fornece a probabilidade de encontrar determinada partícula em uma posição <math>x</math> em certo instante <math>t</math>. Além disso, no presente trabalho temos como objetivo estudar as soluções analítica e numérica da equação para um dado exemplo, através do método FTCS explícito.

== Introdução ==

=== Movimento browniano ===

O movimento browniano foi descoberto pelo botanista Robert Brown, em 1827. Durante seu estudo sobre vida microscópica, ele percebeu pequenas partículas de pólen de plantas se movendo de maneira aleatória no líquido que ele estava estudando e, notando que se tratava de partículas de sujeira, e não seres vivos, chegou a conclusão que era um fenômeno físico, e não biológico, que causava este movimento. <ref name=Feyman_Lecture>The Brownian Movement, The Feynman Lectures on Physics. URL: https://www.feynmanlectures.caltech.edu/I_41.html.</ref>

Posteriormente, foi provado que este fenômeno se dava pelos efeitos do movimento molecular. Em um meio com uma temperatura qualquer <math> >0K </math>, há vibração e movimento molecular. Embora haja conservacão de energia, quando a partícula interage com as moléculas do meio, a energia cinética desta partícula se altera (assim como a das moléculas). Contudo, a soma destas energias é a energia interna do fluido, como descreve o Teorema da Equiparação. <ref name=>The Equipartition Theorem, University of Oxford. URL: http://vallance.chem.ox.ac.uk/pdfs/Equipartition.pdf.</ref>

Em 1905, Albert Einstein propôs uma teoria para descrever tal movimento matematicamente. Primeiramente, ele se propôs a descrever o quão longe uma partícula browniana se desloca em um determinado intervalo de tempo. Como a partícula está sujeita a <math>10^{14}</math> colisões por segundo (provindas das moléculas do meio), a mecânica clássica é incapaz de resolver este sistema <ref name=Feyman_Lecture>The Brownian Movement, The Feynman Lectures on Physics. URL: https://www.feynmanlectures.caltech.edu/I_41.html.</ref><ref>Stachel, J., et al.; The Collected Papers of Albert Einstein, 1989, Princeton University Press. URL: http://users.physik.fu-berlin.de/~kleinert/files/eins_brownian.pdf.</ref>.  Para resolver este problema, ele abordou o problema pela ótica da mecânica estatística. 

Inicialmente, ele considerou o incremento da posição da partícula num espaço unidimensional (<math>x</math>), num determinado intervalo de tempo (<math>\tau</math>). Considerando que existe uma probabilidade aleatória da partícula se mover dentro do intervalo <math>\tau</math>, ele definiu uma função para a densidade de probabilidade (<math>\varphi(\Delta x)</math>). Sabendo que o número de partículas é constante dentro do meio, ele expandiu a densidade (<math>\rho</math>) deste em uma série de Taylor:


:<math>
\rho(x,t) + \tau \left(\frac{\partial \rho(x,t)}{\partial t}\right) + \dots = \rho(x, t+\tau)
=\int_{-\infty}^{+\infty} \rho(x+\Delta x, t) \cdot \varphi(\Delta x) d\Delta x 
</math>


que é, por definição, <math>\varphi</math>. Continuando, 


:<math>
=\rho(x,t) \cdot \int_{-\infty}^{+\infty}  \varphi(\Delta x) d\Delta x + \frac{\partial \rho}{\partial x} \cdot \int_{-\infty}^{+\infty} \Delta x \cdot \varphi(\Delta x) d \Delta x + \frac{\partial^2\rho}{\partial x^2} \int_{-\infty}^{+\infty} \frac{\Delta x}{2} \cdot \varphi(\Delta x) d\Delta x + \dots.
</math>


Pela definição da probabilidade,


:<math>
\int_{-\infty}^{+\infty}  \varphi(\Delta x) d\Delta x=1
</math>


e as integrais dos termos pares da série são nulos devido à simetria do espaço. 


Temos então,

:<math>
\begin{alignat}{2}
\rho(x, t+\tau)
& = \rho(x,t) \cdot 1 + 0 + \frac{\partial^2 \rho}{\partial x^2} \int_{-\infty}^{+\infty} \frac{\Delta x}{2} \cdot \varphi(\Delta x) d\Delta x  + \dots
\\
& = \rho(x,t) + \frac{\partial^2 \rho}{\partial x^2} \int_{-\infty}^{+\infty} \frac{\Delta x}{2} \cdot \varphi(\Delta x) d\Delta x  + \dots.
\end{alignat}
</math>


Esta equação nos leva à igualdade


:<math>
\frac{\partial\rho}{\partial t} =  \frac{\partial^2\rho}{\partial x^2} \cdot \int_{-\infty}^{+\infty} \frac{\Delta x}{2} \cdot \varphi(\Delta x) d\Delta x + \mathcal{O}.
</math>


Podemos interpretar a integral como o coeficiente de difusão <math>D</math>:


:<math>
D = \int_{-\infty}^{+\infty} \frac{\Delta x}{2} \cdot \varphi(\Delta x) d\Delta x.
</math>


O que nos dá a equação da difusão


:<math>
\frac{\partial\rho}{\partial t} = D \frac{\partial^2\rho}{\partial x^2}
</math>

=== Equação de Langevin ===

Em 1908, 3 anos após os estudos de Albert Einstein em processos aleatórios e movimento aleatório, Paul Langevin (1872-1946) apresentou um novo método para o movimento browniano que - segundo Langevin - era "infinitamente mais simples" que a solução proposta por Einstein. <ref name=Langevin_original>P. Langevin, "Sur la théorie de mouvement Brownien" C.R. Acad. Sci. Paris , 146 (1908) pp. 530–533, https://www.physik.uni-augsburg.de/theo1/hanggi/History/Langevin1908.pdf.</ref><ref name=Gardiner>Gardiner, C.W. (1985). Handbook of stochastic methods - for physics, chemistry and the natural sciences, Second Edition. ''Springer series in synergetics''.</ref>
Para interpretar o movimento browniano, Einstein derivou e resolveu uma equação diferencial parcial descrevendo a evolução temporal da densidade de probabilidade para a partícula. Já Langevin aplicou a segunda lei de Newton na forma diferencial para essa partícula.

Para uma partícula browniana de massa <math>m</math> em um líquido com viscosidade <math>\eta</math>, existem duas forças que agem sobre o seu movimento<ref name=Gardiner>Gardiner, C.W. (1985). Handbook of stochastic methods - for physics, chemistry and the natural sciences, Second Edition. ''Springer series in synergetics''.</ref>:

1. Arrasto pela viscosidade.
2. Força de flutuação.

Considerando que a partícula é relativamente grande em comparação com as distâncias médias entre as moléculas do líquido e esta se movimento nesse meio com velocidade <math>v</math>, experimenta uma resistência pela viscosidade. Essa força é descrita pela Lei de Stokes que, para uma partícula esférica com diâmetro <math>a</math>, corresponde a <math>-6\pi\eta a \ dx/dt</math>.

A segunda força foi proposta por Langevin para descrever o efeito de constantes impactos das moléculas do líquido sobre as partículas de estudo. Assim, como essa força possui uma origem aleatório, esta deveria ser positiva ou negativa de maneira equiprovável e cuja magnitude fosse suficiente para manter a agitação da partícula. Caso contrário, a viscosidade iria parar o movimento dessa partícula.

Com isso, a equação que descreve o movimento a partir da posição da partícula - em 1D na direção <math>x</math> - é dado pela Lei de Newton como:


:<math>
m\frac{d^2x}{dt^2} = -\gamma \frac{dx}{dt} + R(t)
</math>

onde <math>m</math> é a massa da partícula, <math>\gamma>0</math> é o coeficiente de fricção devido a viscosidade do líquido e <math>R</math> é postulado a força de Langevin representando as flutuações de pressão devido ao movimento térmico das moléculas que compõem o líquido. <ref>Langevin equation. Encyclopedia of Mathematics. URL: http://encyclopediaofmath.org/index.php?title=Langevin_equation&oldid=47575.</ref>


Essa equação é conhecida como a equação de Langevin e foi o primeiro exemplo de equação diferencial estocástica, isto é, uma equação diferencial com um termo, <math>R(t)</math> nesse caso, cuja solução em algum sentido também é uma função aleatória <ref name=Gardiner>Gardiner, C.W. (1985). Handbook of stochastic methods - for physics, chemistry and the natural sciences, Second Edition. ''Springer series in synergetics''.</ref>. Essa função foi desenvolvida para possuir as seguintes propriedades:


:<math>
\langle R(t) \rangle = 0
\qquad
\langle R(t) R(t') \rangle = D \cdot \delta(t-t')
</math>

onde <math>\langle f \rangle</math> descreve o valor médio ou esperado de uma função <math>f</math>, <math>D</math> é o coeficiente de difusão e <math>\delta</math> é a função delta.


A primeira propriedade afirma que o movimento é aleatório de forma que não existe nenhuma tendência de sentido para a partícula se locomover. Assim, é dito que <math>R(t)</math> se trata de um ruído branco gaussiano. Já a segunda propriedade mostra que a força em um dado tempo <math>t</math> é descorrelacionada de uma força para qualquer outro tempo <math>t'</math> (propriedade de Markov). <ref>Wikipédia: Langevin equation. URL: https://en.wikipedia.org/wiki/Langevin_equation.</ref>

=== Dedução da equação de Fokker-Planck ===

Para a dedução foi utilizado a seguinte referência como base: ''From Langevin to Fokker-Planck equation''.<ref>From Langevin to Fokker-Planck equation. URL: http://cgl.elte.hu/~racz/Stoch-diff-eq.pdf.</ref>

Como vimos anteriormente, o movimento browniano pode ser descrito pela equação de Langevin, a qual podemos resolver sem nenhum problema. Contudo, alternativamente, podemos fazer uso da equação de Fokker-Planck e considerar uma densidade de probabilidade em relação a posição e o tempo, levando em conta diferentes perturbações estocásticas. Para isso, iremos agora deduzir essa equação, a partir da analise de Langevin.

A descrição do movimento browniano foi anteriormente tratada em termos da equação


:<math>
m\frac{d^2x}{dt^2} = -\gamma \frac{dx}{dt} + R(t)
</math>

onde <math>m</math> é a massa da partícula, <math>\gamma>0</math> é o coeficiente de fricção devido a viscosidade do líquido e <math>R</math> é postulado a força de Langevin.


Para o nosso propósito, realizaremos uma abordagem mais direta e geral, assumindo a equação de Langevin através de um potencial <math>U(x)</math>, nos deixando com a seguinte expressão:


:<math>
\frac{dx}{dt} = -\mu \frac{dU}{dx} + \eta(t),
</math>

onde <math>\mu=\frac{1}{6\pi\eta a}</math> e <math>\eta(t) = \mu R(t)</math>.


Contudo, para entender as equações diferenciais estocásticas, precisamos discretizar as equações para posteriormente aplicar os métodos numéricos. Temos que o sistema de discretização mais simples possível nos leva ao valor de <math>x</math> no tempo <math>t + \varepsilon</math> a partir da interação


:<math>
x(t + \varepsilon) = x(t) - \mu\frac{dU}{dx}\varepsilon + \eta_\varepsilon(t).
</math>


A tarefa é compreender como o agora ruído gaussiano <math>\eta_{\varepsilon}(t)</math> afeta esse sistema discretizado. Para isso, escreveremos a evolução temporal da distribuição de probabilidade <math>\rho(x,t)</math> definida como a probabilidade da partícula estar em <math>x</math> no tempo <math>t</math>. Então, deve-se olhar para o limite <math>t \rightarrow \infty</math> de <math>\rho(x,t)</math>.
Vamos realizar a tarefa com a equação discretizada de Langevin comentada acima. Por razões de simplificação, discutiremos essa abordagem com a forma 


:<math>
x(t + \varepsilon) = x(t) + v(x) \varepsilon + \eta_{\varepsilon}(t)
</math>

onde <math>v(x) \equiv −\mu \frac{dU}{dx}</math>.


A equação da probabilidade <math>\rho(x,t)</math> pode ser derivada usando a equação de Chapman-Kolmogorov, a qual também iremos escrever em uma forma discretizada no tempo, da forma


:<math>
\rho(x,t + \varepsilon) = \int_{-\infty}^{+\infty} W(x,y;\varepsilon) \rho(y,t) dy
</math>

onde <math>W(x,y; \varepsilon)</math> é a probabilidade de que a partícula se mova para <math>x</math> no instante <math>t + \varepsilon</math>, desde que tenha começado em <math>y</math> no instante <math>t</math>. 


Lembrando que: para que a condição acima mova a partícula de <math>y</math> a <math>x</math> no tempo <math>\varepsilon</math>, o termo de ruído deve ter apenas o valor correspondente que satisfaz a equação


:<math>
x = y + v(y) \varepsilon + \eta_{\varepsilon}(t).
</math>


A probabilidade de <math>\eta_{\varepsilon}(t)</math> acima é dada a partir da distribuição de probabilidade do incremento estocástico <math>\eta{\varepsilon}
</math>

:<math>
\rho_{G} (\eta \varepsilon) = \frac{1}{\sqrt{4\pi D\varepsilon}} \exp\left({\frac{-\eta_\varepsilon^2}{4D\varepsilon}}\right),
</math>


Assim <math>W(x,y; \varepsilon)</math> é obtido como


:<math>
\begin{alignat}{2}
W(x,y; \varepsilon) & = \rho_{G} (\eta_{\varepsilon} = x − y − v(y)\varepsilon) \\
& = \frac{1}{\sqrt{4\pi D\varepsilon}} \exp\left({\frac{−[x−y−v(y)\varepsilon]^2}{4Dε}}\right).
\end{alignat}
</math>


Substituindo o <math>W(x,y; \varepsilon)</math> acima na equação de Chapman-Kolmogorov discretizada e expandindo o lado esquerdo para <math>\varepsilon</math>, temos


:<math>
\rho(x,t) + \partial_{t}\rho(x,t) \varepsilon = \int_{-\infty}^{+\infty} \frac{\exp\left(\frac{-[x-y-v(y) \varepsilon]^2}{4D\varepsilon}\right)}{\sqrt{4\pi D \varepsilon}} \rho(y,t) dy
</math>


Como último passo agora, vamos expandir a integral do lado direito para ordenar <math>\varepsilon</math>. Notamos que quando <math>\varepsilon \rightarrow 0</math>, a Gaussiana vai para uma função delta e a integral resulta em <math>\rho(x,t)</math> que irá se cancelar com <math>\rho(x,t)</math> no lado esquerdo. Ao final, quando reorganizamos os termos, temos a seguinte equação denominada de Equação de Fokker-Planck


:<math>
\frac{\partial \rho(x,t)}{\partial t} = - \frac{\partial v(x) \rho(x,t)}{\partial x} + D \frac{\partial^2 \rho(x,t)}{\partial x^2}
</math>

ou, se substituirmos <math>v(x)</math> pelo potencial <math>v(x) = −\mu \partial U(x)/\partial x</math>, temos que


:<math>
\frac{\partial \rho(x,t)}{\partial t} = \frac{\partial}{\partial x} \frac{\partial U}{\partial x} \rho(x,t) + D \frac{\partial^2 \rho(x,t)}{\partial x^2}
</math>

onde, em ambos os casos, <math>D</math> é o coeficiente de difusão no movimento browniano, com valor <math>D = \mu k_B T = \frac{k_B T}{6 \pi \eta a}</math>, cujas variáveis foram discutidas anteriormente na equação de Langevin.


Esta distribuição pode ainda depender de um conjunto de <math>N</math> macrovariáveis <math>x_i</math>, de tal maneira que o movimento browniano em questão pode ser representado por uma equação de Fokker-Planck da forma <ref>Wikipédia: Fokker-Planck Equation. URL: https://en.wikipedia.org/wiki/Fokker-Planck_equation.</ref>:


:<math>
\frac{\partial \rho(x,t)}{\partial t} = - \left[\sum_{i=1}^{N} \frac{\partial}{\partial x_{i}} D_{i}^{1} (x_{1},...,x_{N}) + \sum_{i=1}^{N} \sum_{j=1}^{N} \frac{\partial^2}{\partial x_{i} \partial x_{j}} D_{ij}^{2} (x_{i},...,x_{N})\right]\rho,
</math>

onde <math>D^1</math> é o termo que é dado por um vetor e <math>D^2</math> é o termo de difusão, dado por uma matriz.

== Exemplos: acelerador de partículas ==

O exemplo e explicações trazidas nessa seção foram retirados do artigo: ''Numerical solution of two dimensional Fokker-Planck equations'' <ref>M.P. Zorzano, H. Mais, L. Vazquez, Numerical solution of two dimensional Fokker—Planck equations, Applied Mathematics and Computation, 1999, https://doi.org/10.1016/S0096-3003(97)10161-8.</ref>.

Uma aplicação amplamente utilizada para a equação de Fokker-Planck como uma Equação Diferencial Estocástica (''Stochastic Differential Equations'', SDE) é para a modelagem da dinâmica de partículas únicas em aceleradores sob influência de ruído. As fontes de ruído são, por exemplo, campos aleatórios, movimento aleatório do solo ou flutuações quânticas da radiação.

A partir do cálculo de probabilidade fornecido pela equação de Fokker-Planck, algumas propriedades físicas de interesse são: qual é o comportamento de longo prazo da dinâmica; qual a probabilidade da partícula colidir com a câmara de vácuo e então ser perdida; quais são as flutuações médias da partícula em torno da trajetória do acelerador; e qual é a evolução temporal da densidade de probabilidade?

Para responder essas perguntas, apenas alguns casos mais específicos possuem soluções analíticas, especialmente em problemas com maior dimensionalidade. Por esse motivo, é interessante e prático estudar a equação de Fokker-Planck e os problemas relacionados a partir de métodos numéricos.

Para esse fim, serão abordados dois exemplos: o oscilador harmônico amortecido, o qual possuí solução analítica e poderá ser usado como comparação para verificar a acurácia do método; e o problema de Duffing, oscilador não linear, com amortecimento e ruído multiplicativo e aditivo, que não apresenta solução exata.

=== Equações ===

As equações diferenciais típicas encontradas na física do acelerador de partícula são


:<math>
\begin{cases}
\dot{x} = v,
\\
\dot{v} = -a_1(x) - a_2(x,v) + \sqrt{2\sigma}R(t),
\end{cases}
</math>

onde <math>R(t)</math> é o ruído branco gaussiano e <math>\sigma</math> é uma constante.


Para o caso do oscilador harmônico amortecido e com ruído, temos


:<math>
\begin{cases}
\dot{x} = v,
\\
\dot{v} = -Kx - \gamma v + \sqrt{2\sigma}R(t),
\end{cases}
</math>

onde <math>K</math> e <math>\gamma</math> também são constantes.


Já para o oscilador estocástico de Duffing, não-linear, com amortecimento, ruído multiplicativo e aditivo temos


:<math>
\begin{cases}
\dot{x} = v,
\\
\dot{v} = -\omega^2 \left[ \left( \alpha + \sqrt{2D_{11}}R_1(t) \right) x + \epsilon x^3 \right] - 2\tau \omega \dot{x} + \sqrt{2D_{22}}R_2(t),
\end{cases}
</math>

onde <math>\omega</math>, <math>\alpha</math>, <math>\epsilon</math>, <math>\tau</math>, <math>D_{11}</math> e <math>D_{22}</math> são constantes dadas.


Escrevendo esses exemplos em termos da equação de Fokker-Planck obtemos


:<math>
\frac{\partial \rho}{\partial t} = -\frac{\partial v}{\partial x}\rho + \frac{\partial}{\partial v} \left(\gamma v + Kx + \sigma \frac{\partial}{\partial v}\right)\rho
</math>

para o primeiro caso e


:<math>
\frac{\partial \rho}{\partial t} = -\frac{\partial v}{\partial x}\rho + \frac{\partial}{\partial v} [(a_1(x) + a_2(x,v))\rho] + \frac{1}{2}\frac{\partial^2}{\partial v^2}(2D\rho)
</math>

:<math>
a_1(x) =\omega^2x(\alpha+\epsilon x^2), \quad a_2(x,v) = 2\tau\omega v, \quad D=D_{11}\omega^4x^2+D_{22}
</math>

para o segundo exemplo.

=== Método FTCS explícito ===

Para os fins de estudo desses casos, o método numérico mais simples para descrever a função <math>\rho</math> é a partir do FTCS (''Foward Time Central Space'') explícito <ref>Pichler L., Masud A., Bergman L.A. (2013) Numerical Solution of the Fokker–Planck Equation by Finite Difference and Finite Element Methods—A Comparative Study. In: Papadrakakis M., Stefanou G., Papadopoulos V. (eds) Computational Methods in Stochastic Dynamics. Computational Methods in Applied Sciences, vol 26. Springer, Dordrecht. https://doi.org/10.1007/978-94-007-5134-7_5.</ref>. Como ficará mais claro a partir da discretização das equações, o método FTCS calcula o resultado de uma função no tempo seguinte em um dado ponto a partir dos valores dos "vizinhos" desse ponto. Para a técnica explícita, podemos calcular valor da função no tempo <math>n+1</math> a partir apenas de valores já conhecidos no tempo atual <math>n</math>. Assim, a aplicação do método se torna bem simples, contudo, perdemos precisão/estabilidade em comparação com métodos mais sofisticados.

O primeiro passo para o método é discretizar os termos da função de interesse <math>\rho(x,v,t)</math> nas derivadas parciais necessárias: <math>\partial/\partial x</math>, <math>\partial/\partial v</math>, <math>\partial^2/\partial v^2</math>. Para isso, usamos a expansão da série de Taylor e, para melhor precisão, a derivada centrada (excluindo a derivada no tempo). Ademais, usaremos a notação dos índices para indicar os valores de <math>x</math>, <math>v</math> e <math>t</math> por <math>i</math>, <math>j</math> e <math>n</math>, respectivamente.


:<math>
\frac{\partial \rho}{\partial t}
\rightarrow
\frac{\rho(x,v,t+\Delta t) - \rho(x,v,t)}{\Delta t}
+ \mathcal{O}(\Delta t^2)
\rightarrow
\frac{\rho_{i,j}^{n+1}-\rho_{i,j}^n}{\Delta t}
</math>


:<math>
\frac{\partial \rho}{\partial x}
\rightarrow
\frac{1}{2}
\left[ 
    \frac{\rho(x+\Delta x,v,t) - \rho(x,v,t)}{\Delta x} 
    +
    \frac{\rho(x-\Delta x,v,t) - \rho(x,v,t)}{\Delta x}
\right]
+ \mathcal{O}(\Delta x^2)
\approx
\frac{\rho(x+\Delta x,v,t)-\rho(x-\Delta x,v,t)}{2 \Delta x}
\rightarrow
\frac{\rho_{i+1,j}^n + \rho_{i-1,j}^n}{2\Delta x}
</math>


:<math>
\frac{\partial \rho}{\partial v}
\rightarrow
\frac{1}{2}
\left[ 
    \frac{\rho(x,v+\Delta v,t) - \rho(x,v,t)}{\Delta v} 
    +
    \frac{\rho(x,v,t) - \rho(x,v-\Delta v,t)}{\Delta v}
\right]
+ \mathcal{O}(\Delta v^2)
\approx
\frac{\rho(x,v+\Delta v,t)-\rho(x,v-\Delta v,t)}{2 \Delta v}
\rightarrow
\frac{\rho_{i,j+1}^n + \rho_{i,j-1}^n}{2\Delta v}
</math>


:<math>
\frac{\partial^2 \rho}{\partial v^2}
\rightarrow
\frac{\rho(x,v+\Delta v, t)-2\rho(x,v,t)+\rho(x,v-\Delta v,t)}{(\Delta v)^2}
+\mathcal{O}(\Delta v^3)
\rightarrow
\frac{\rho_{i,j+1}^n-2\rho_{i,j}^n+\rho_{i,j-1}^n}{(\Delta v)^2}
</math>

=== Aplicação ===

Com as derivadas parciais discretizadas, podemos agora reescrever as funções de Fokker-Planck para cada caso de estudo em busca do valor de <math>\rho_{i,j}^{n+1}</math>, isto é, o valor de cada ponto no espaço <math>xv</math> para o tempo <math>t+\Delta t</math>.

==== Caso 1: oscilador harmônico amortecido ====

Primeiramente, reescrevemos a função original expandindo os termos e aplicando a regra da cadeia onde for necessário.


:<math>
\begin{alignat}{2}
\frac{\partial \rho}{\partial t}
& =
-\frac{\partial v}{\partial x}\rho + \frac{\partial}{\partial v} \left(\gamma v + Kx + \sigma \frac{\partial}{\partial v}\right)\rho
\\[1em]
& =
- \frac{\partial}{\partial x}(v \rho)
+ \frac{\partial}{\partial v}(\gamma v + Kx)\rho
+ \sigma \frac{\partial^2}{\partial v^2}\rho
\\[1em]
& =
- v \frac{\partial}{\partial x} \rho
+ \frac{\partial}{\partial v} (\gamma v \rho) 
+ \frac{\partial}{\partial v} (Kx \rho)
+ \sigma \frac{\partial^2}{\partial v^2}\rho
\\[1em]
& =
- v \frac{\partial}{\partial x} \rho
+ \gamma \rho
+ \gamma v \frac{\partial}{\partial v} \rho  
+ Kx \frac{\partial}{\partial v} \rho
+ \sigma \frac{\partial^2}{\partial v^2}\rho
\\[1em]
& =
- v \frac{\partial}{\partial x} \rho
+ \gamma \rho
+ (\gamma v+Kx) \frac{\partial}{\partial v} \rho  
+ \sigma \frac{\partial^2}{\partial v^2}\rho
\end{alignat}
</math>


Agora, reescrevemos a mesma função substituindo as derivadas parciais e demais variáveis necessárias (<math>x</math>) pelas versões discretizadas.


:<math>
\frac{\rho_{i,j}^{n+1}-\rho_{i,j}^n}{\Delta t}
=
- v \left( \frac{\rho_{i+1,j}^n-\rho_{i-1,j}^n}{2\Delta x} \right)
+ \gamma \rho_{i,j}^n
+ (\gamma v +Kx^n) \left( \frac{\rho_{i,j+1}^n-\rho_{i,j-1}^n}{2\Delta v} \right)
+ \sigma \left( \frac{\rho_{i,j+1}^n-2\rho_{i,j}^n+\rho_{i,j-1}^n}{(\Delta v)^2} \right)
</math>


:<math>
\rho_{i,j}^{n+1}
=
\rho_{i,j}^n
+ \Delta t
\left[
- v \left( \frac{\rho_{i+1,j}^n-\rho_{i-1,j}^n}{2\Delta x} \right)
+ \gamma \rho_{i,j}^n
+ (\gamma v +Kx^n) \left( \frac{\rho_{i,j+1}^n-\rho_{i,j-1}^n}{2\Delta v} \right)
+ \sigma \left( \frac{\rho_{i,j+1}^n-2\rho_{i,j}^n+\rho_{i,j-1}^n}{(\Delta v)^2} \right)
\right]
</math>


Para exercício, usaremos um <math>v</math> constante de forma que apenas será necessário atualizar a variável <math>x</math> para cada passo de tempo.


:<math>
x^{n+1} = x^n+v\Delta t
</math>


Por fim, é necessário definir uma condição inicial e uma condição de contorno. Uma vez que a função descreve uma densidade de probabilidade, a literatura sugere, até para comparar com a solução analítica quando disponível, uma delta de Kronecker <math>\delta(x-x_0,v-v_0)</math>. Assim, em uma dada coordenada <math>(x,v)</math> no tempo <math>t_0</math> teremos uma probabilidade 1 de encontrar a partícula. Como condição de contorno, dado que a função em um dado ponto precisa dos pontos vizinhos, indicamos que a "borda" da matriz será sempre igual a zero.


Dado essas equações, escrevemos um programa em C para gerar os dados e o gnuplot para criar o gráfico/animação. Os códigos desenvolvidos estarão disponíveis na próxima seção "Código desenvolvido".


<center>
<li style="display: inline-block;">[[Arquivo:caso1.gif|800px|thumb|middle|Animação da evolução da probabilidade de encontrar uma partícula com dada posição <math>x</math> e <math>v</math> a partir da equação de um oscilador harmônico amortecido. Valores da animação: <math>\Delta x = 1\times 10^{-3}</math>, <math>\Delta v = 1\times 10^{-3}</math>, <math>\Delta t = 1\times 10^{-6}</math>, <math>x_0=1</math>, <math>v=2</math>, <math>k=0.5</math>, <math>\gamma=0.5</math>, <math>\sigma=0.5</math>.]]</li>
</center>


A partir da referência estudada e da física do problema, é evidente que a animação apresentada está aquém do esperado. Pelos testes pelo grupo responsável por esse documento, o programa desenvolvido se mostrou muito instável e impreciso: alterações mínimas de alguns parâmetros alteravam consideravelmente o resultado observado. Esse comportamento "imprevisível" se tornou ainda evidente ao alterar o tamanho dos passos (''dt'', ''dx'', ''dv'') de forma que tanto aumentando quanto diminuindo uma ordem de grandeza a função "explodia" - isto é - crescia exponencialmente. Assim, não foi possível aumentar a resolução da "simulação".

Além disso, em alguns casos, a função retornava valores negativos, o que é inadequado para uma densidade de probabilidade. Ademais, a ordem de grandeza para o tempo usada na simulação desse trabalho está muito abaixo da utilizada na referência base para esse problema, mostrando que a evolução temporal pode ter algum problema. No entanto, o problema que destacamos é a falta de "suavidade" na curva gerada, tanto é que em alguns quadros é possível ver mudanças bruscas no formato da função.

Dessa forma, apesar dos esforços da equipe, não fomos capazes - até o momento dessa edição - de encontrar a causa desses problemas. No entanto, segue algumas hipóteses que podem ser avaliadas para futuras revisões deste material:

- Verificar as equações utilizadas para descrever os problemas e a passagem para a equação de Fokker-Planck. A literatura encontrada não apresenta os passos usados para esse cálculo e nem a referência específica para cada etapa. Assim, o grupo precisou confiar que esse artigo estava correto.

- Verificar a discretização da função para o método FTCS explícito.

- Verificar a estabilidade do método, pois, apesar de não ter ser comentado no material de referência, é possível que o método explícito não seja adequado para os casos estudados.

- Verificar possíveis problemas no código, tanto o código para os dados quanto para a visualização.

- Verificar problemas de truncamento ou erro numérico ou ainda valores inadequados para as constantes utilizadas.


==== Caso 2: problema de Duffing ====

Assim como no caso anterior, reescrevemos a equação para uma forma conveniente.


<math>
\frac{\partial \rho}{\partial t} = -\frac{\partial v}{\partial x}\rho + \frac{\partial}{\partial v} [(a_1(x) + a_2(x,v))\rho] + \frac{1}{2}\frac{\partial^2}{\partial v^2}(2D\rho)
</math>

<math>
a_1(x) =\omega^2x(\alpha+\epsilon x^2), \quad a_2(x,v) = 2\tau\omega v, \quad D=D_{11}\omega^4x^2+D_{22}
</math>


<math>
\begin{alignat}{2}
\frac{\partial \rho}{\partial t} 
& =
- \frac{\partial v}{\partial x}\rho 
+ \frac{\partial}{\partial v} [\omega^2x(\alpha+\epsilon x^2) + 2\tau\omega v]\rho
+ \frac{\partial^2}{\partial v^2}(D_{11}\omega^4x^2+D_{22})\rho
\\[1em]
& =
- v\frac{\partial}{\partial x}\rho
+ [\omega^2x(\alpha+\epsilon x^2)+2\tau\omega v]\frac{\partial}{\partial v}\rho 
+ 2\tau\omega\rho
+ (D_{11}\omega^4x^2+D_{22})\frac{\partial^2}{\partial v^2}\rho
\end{alignat}
</math>


Novamente, substituímos os termos necessários pela forma discretizada.


<math>
\frac{\rho_{i,j}^{n+1}-\rho_{i,j}^n}{\Delta t}
=
- v \left( \frac{\rho_{i+1,j}^n-\rho_{i-1,j}^n}{2\Delta x} \right)
+ [\omega^2{x^n}(\alpha+\epsilon {x^n}^2)+2\tau\omega v]
\frac{\rho_{i,j+1}^n + \rho_{i,j-1}^n}{2\Delta v}
+ 2\tau\omega \rho_{i,j}^n
+ (D_{11}\omega^4{x^n}^2+D_{22})
\frac{\rho_{i,j+1}^n-2\rho_{i,j}^n+\rho_{i,j-1}^n}{(\Delta v)^2}
</math>


<math>
\rho_{i,j}^{n+1}
=
\rho_{i,j}^n
+
\Delta t
\left\{
- v \left( \frac{\rho_{i+1,j}^n-\rho_{i-1,j}^n}{2\Delta x} \right)
+ [\omega^2{x^n}(\alpha+\epsilon {x^n}^2)+2\tau\omega v]
\frac{\rho_{i,j+1}^n + \rho_{i,j-1}^n}{2\Delta v}
+ 2\tau\omega \rho_{i,j}^n
+ (D_{11}\omega^4{x^n}^2+D_{22})
\frac{\rho_{i,j+1}^n-2\rho_{i,j}^n+\rho_{i,j-1}^n}{(\Delta v)^2}
\right\}
</math>


Considerando mais uma vez que a velocidade da partícula será constante,podemos atualizar o valor de <math>x</math> pela expressão: <math>x^{n+1}=x^n+v\Delta t</math>.

Para condições iniciais e condições de contorno, usamos o mesmo princípio definido no caso anterior: condição inicial dada por <math>\delta(x-x_0,v-v_0)</math> e condição de contorno definindo as "paredes" com valor zero para qualquer tempo <math>t</math>.


Uma vez que o grupo não foi capaz de encontrar bons resultados para o caso anterior (mais simples), não foi desenvolvido nenhum código para esse exemplo até que o problema anterior fosse resolvido. Esperamos que em uma futura edição esse material possa trazer alguma visualização para o caso com ruído aditivo e multiplicativo descrito pelo problema de Duffing.

== Código desenvolvido ==

Função auxiliar utilizada do arquivo de cabeçalho ''pointers.h'':
<source lang='c'>
double ***create_3d_double_pointer(int x, int y, int z){
  double ***pointer;
  int i, j;

  assert((x > 0) && (y > 0) && (z > 0));
  pointer = (double ***)malloc(x * y * z * sizeof(double **));
  if (pointer == NULL){
    fprintf(stderr, "Error in routine create_3d_double_pointer\n");
    exit(8);
  }
  else {
    for (i = 0; i < x; i++){
      pointer[i] = (double **)malloc(y * z * sizeof(double *));
      if (pointer[i] == NULL){
        fprintf(stderr, "Error in routine create_3d_double_pointer\n");
        exit(8);
      }
    }
    for (i = 0; i < x; i++){
      for (j = 0; j < y; j++){
        pointer[i][j] = (double *)malloc(z * sizeof(double));
        if (pointer[i][j] == NULL){
          fprintf(stderr, "Error in routine create_3d_double_pointer\n");
          exit(8);
        }
      }
    }
  }
  return (pointer);
}
</source>


Código utilizado para gerar os dados para o "caso 1". 


<source lang='c'>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pointers.h"

int main(){

  // Declara as variáveis
  FILE *fp;
  char output[50];
  double ***f;
  int i, j, k, c, Lx, Lv;
  double time, tmax, dt, dx, dv, v, x, K, gamma, sigma, epsilon; 
  
  // Configura os valores iniciais
  time = 0.0;       // Tempo
  tmax = 0.0001;    // Tempo máximo
  Lx = 50;          // Quantidade de pontos em x
  Lv = 100;         // Quantidade de pontos em v
  dt = 0.000001;    // Passo em t
  dx = 0.001;       // Passo em x
  dv = 0.001;       // Passo em v
  x = 1;            // Posição inicial
  v = 2;            // Velocidade
  K = 0.5;          // Constante
  gamma = 0.5;      // Constante
  sigma = 0.5;      // Constante
  epsilon = 10E-4;  // Valor para comparação

  // Inicializa a matriz 3D
  f = create_3d_double_pointer(Lx, Lv, 2);
  // Todos os pontos com valores zero
  for (i=0; i<Lx; i++){
    for (j=0; j<Lv; j++){
      for (k=0; k<2; k++){
        f[i][j][k] = 0.0;
      }
    }
  }

  // Condição inicial
  f[Lx/2][Lv/2][0] = 1.0;

  // Arquivo para output
  sprintf(output, "caso1.dat");
  fp = fopen(output, "w");

  c = 0;
  // Loop da "simulação"
  while (time<tmax){
    for (i=1; i<(Lx-1); i++){
      for (j=1; j<(Lv-1); j++){
        // Atualiza a probabilidade
        f[i][j][1-c] = f[i][j][c]+dt*(-v*(f[i+1][j][c]-f[i-1][j][c])/(2*dx)+gamma*f[i][j][c]+(gamma*v+K*x)*(f[i][j+1][c]-f[i][j-1][c])/(2*dv)+sigma*(f[i][j+1][c]-2*f[i][j][c]+f[i][j-1][c])/(dv*dv));
        // Salva o resultado em arquivo
        if (fabs(f[i][j][1-c]) > epsilon){
          fprintf(fp, "%f %f %f\n", (double)(i-Lx/2)*dx+x, (double)(j-Lv/2)*dv+v, f[i][j][1-c]);
        }
      }
      // Adiciona linha em branco (arquivo XYZ - gnuplot)
      if (fabs(f[i][j][1-c]) > epsilon){
        fprintf(fp, "\n");
      }
    }
    // Atualiza a posição
    x = x + v*dt;
    c = 1 - c;
    // Adiciona duas linhas em branco (animação - gnuplot)
    fprintf(fp, "\n\n");
    time += dt;
  }

  fclose(fp); // Fecha arquivo
  free(f);    // Libera ponteiro

  return 0;
}
</source>


Código utilizado no gnuplot para gerar a animação.


<source lang='gnuplot'>
# Limpa configurações anteriores
unset xrange
unset yrange
unset zrange

# Gera estatísticas dos dados
stats "caso1.dat" name "A"

# Estilo para gráfico
set dgrid3d

# Ponto de vista
set view 60, 30, 1, 1

# Configura o range do plot
set xrange [A_min_x:A_max_x]
set yrange [A_min_y:A_max_y]
set zrange [0:1]

# Formato de saída
set term gif animate delay 25 background "#212121" enhanced font "Arial,10" fontscale 1.0 size 800, 550

set title "Partícula em acelerador:\noscilador harmônico amortecido" offset character 0, -3 font "Oswald,28" tc rgb "white"
set xlabel "Posição (x)" offset character -2, -2 font ",18" tc rgb "white" rotate by -9
set ylabel "Velocidade (v)" offset character -3, -3 font ",18" tc rgb "white" rotate parallel
set zlabel "Densidade de probabilidade" offset character -1.5, -1.5 font ",18" tc rgb "white" rotate parallel

set tics tc rgb "white" offset character -1, -1
set border lc rgb "white"

set output "caso1.gif"

do for [i=0:int(A_blocks-1)] {  
    splot 'caso1.dat' index i with pm3d notitle
    }
</source>

== Referências ==
<references/>