Física Computacional - Contribuições do usuário [pt-br]

Equação de Lotka-Volterra Competitiva Estocástica

2024-09-03T02:07:35Z

Gustavobpasquali:

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos. Por fim, descreveremos o método de Fokker-Planck, utilizando um exemplo simples para a demonstração.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Método de Fokker-Planck==
A Equação de Fokker-Planck trata da equação diferencial da distribuição de probabilidade das variáveis de um sistema de Equações Diferenciais Estocásticas. Avaliaremos sua aplicação no sistema Lotka-Volterra competitivo para duas populações.

Para um processo estocástico vetorial <math> X = [X_1, X_2, ... , X_n] </math>, a equação de Fokker-Planck tem a forma

<center>
<math>
\frac{\partial P(x,t)}{\partial t} = \sum^n_{i=1} \frac{\partial[K_i(x,t)P(x,t)]}{\partial x_i} + \frac{1}{2} \sum^n_{i,j = 1}\frac{\partial^2[D_{ij}(x,t)P(x,t)]}{\partial x_i\partial x_j}
</math>
</center>

em que

<center>
<math>
K_i(x,t) = -A_i(x,t)
</math>
</center>

e
<center>

<math>
D_{ij}(x,t) = \sum^m_{\alpha=1}B_{i,\alpha}(x,t)B_{j,\alpha}(x,t)
</math>
</center>

Assim, para o nosso caso, temos

<center>
<math>
\frac{\partial P(x_1,x_2,t)}{\partial t} = -\frac{\partial}{\partial x_1}\left[r_1x_1\left(1-\frac{x_1+\alpha_{12}x_2}{k_1}\right)P(x_1,x_2,t) \right] -\frac{\partial}{\partial x_2}\left[r_2x_2\left(1-\frac{x_2+\alpha_{21}x_1}{k_2}\right)P(x_1,x_2,t) \right] +\frac 12 \beta_1^2 \frac{\partial^2x_1^2P(x_1,x_2,t)}{\partial x_1^2} +\frac 12 \beta_2^2 \frac{\partial^2x_2^2P(x_1,x_2,t)}{\partial x_2^2}
</math>
</center>

Resolvendo as derivadas por regra do produto e separando nas derivadas parciais de <math> P(x_1,x_2,t) </math> ficamos com
<center>
<math>

\frac{\partial P(x_1,x_2,t)}{\partial t } = \delta_1 \frac{\partial^2 P}{\partial x_1^2} + \delta_2 \frac{\partial^2 P}{\partial x_2^2} + \gamma_1 \frac{\partial P}{\partial x_1} + \gamma_2 \frac{\partial P}{\partial x_2} + \omega P

</math>
</center>

com

<center>
<math>
\delta_1 = \frac{1}{2} \beta_1^2 x_1^2
</math>

<math>
\delta_2 = \frac{1}{2} \beta_2^2 x_2^2
</math>

<math>
\gamma_1 = \frac{r_1 x_1^2}{k_1} + \frac{r_1 x_1 \alpha_{12} x_2}{k_1} - r_1x_1 + 2\beta_1^2x_1
</math>

<math>
\gamma_2 = \frac{r_2 x_2^2}{k_2} + \frac{r_2 x_2 \alpha_{21} x_1}{k_2} - r_2x_2 + 2\beta_2^2x_2
</math>

<math>
\omega = \frac{2r_1x_1}{k_1} + \frac{2r_2x_2}{k_2} + \frac{r_1\alpha_{12}x_2}{k_1} + \frac{r_2\alpha_{21}x_1}{k_2} - r_1 - r_2 + \beta_1^2 + \beta_2^2
</math>

</center>

Discretizando a equação diferencial acima, usando que

<center>
<math>
\frac{\partial P(x,t)}{\partial t} = \frac{P_{i,j}^{n+1} - P_{i,j}^{n}}{dt}
</math>

<math>
\frac{\partial P(x,t)}{\partial x_1} = \frac{1}{2} \frac{P_{i+1,j}^{n} - P_{i-1,j}^{n}}{dx_1}
</math>

<math>
\frac{\partial P(x,t)}{\partial x_2} = \frac{1}{2} \frac{P_{i,j+1}^{n} - P_{i,j-1}^{n}}{dx_2}
</math>

<math>
\frac{\partial^2 P(x,t)}{\partial x_1^2} = \frac{P_{i+1,j}^{n} + P_{i-1,j}^{n} - 2P_{i,j}^{n}}{dx_1^2}
</math>

<math>
\frac{\partial^2 P(x,t)}{\partial x_2^2} = \frac{P_{i,j+1}^{n} + P_{i,j-1}^{n} - 2P_{i,j}^{n}}{dx_2^2}
</math>

</center>

e que <math> dx_1 = dx_2 </math>, chegamos, enfim, em

<center>
<math>
P_{i,j}^{n+1} = P_{i,j}^{n}\left( 1 - \frac{2\delta_1 dt}{dx^2} - \frac{2\delta_2 dt}{dx^2} + \omega dt \right) + P_{i+1,j}^{n} \left( \frac{\delta_1 dt}{dx^2} + \frac{\gamma_1 dt}{2dx} \right) + P_{i,j+1}^{n} \left( \frac{\delta_2 dt}{dx^2} + \frac{\gamma_2 dt}{2dx} \right) + P_{i-1,j}^{n} \left( \frac{\delta_1 dt}{dx^2} - \frac{\gamma_1 dt}{2dx} \right) + P_{i,j-1}^{n} \left( \frac{\delta_2 dt}{dx^2} - \frac{\gamma_2 dt}{2dx} \right)
</math>
</center>

que foi resolvida numericamente.

Abaixo, temos um gráfico considerando duas populações, com constantes todas iguais e baixo ruído. O eixo das abcissas indica a população <math>x_1</math>, o eixo das ordenadas, a população <math>x_2</math> e a cor representa a probabilidade da população possuir aquele número de indivíduos.

<center>
[[Arquivo:Tmp.gif]]
</center>

== Código ==

Os principais programas utilizados no trabalho se encontram no repositório https://github.com/gugabXD/trab2-metcompC

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T17:47:57Z

Gustavobpasquali: /* Método de Fokker-Planck */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Método de Fokker-Planck==
A Equação de Fokker-Planck trata da equação diferencial da distribuição de probabilidade das variáveis de um sistema de Equações Diferenciais Estocásticas. Avaliaremos sua aplicação no sistema Lotka-Volterra competitivo para duas populações.

Para um processo estocástico vetorial <math> X = [X_1, X_2, ... , X_n] </math>, a equação de Fokker-Planck tem a forma

<center>
<math>
\frac{\partial P(x,t)}{\partial t} = \sum^n_{i=1} \frac{\partial[K_i(x,t)P(x,t)]}{\partial x_i} + \frac{1}{2} \sum^n_{i,j = 1}\frac{\partial^2[D_{ij}(x,t)P(x,t)]}{\partial x_i\partial x_j}
</math>
</center>

Onde

<center>
<math>
K_i(x,t) = -A_i(x,t)
</math>
</center>
e
<center>
<math>
D_{ij}(x,t) = \sum^m_{\alpha=1}B_{i,\alpha}(x,t)B_{j,\alpha}(x,t)
</math>
</center>

Assim, para o nosso caso, temos:
<center>
<math>
\frac{\partial P(x_1,x_2,t)}{\partial t} = -\frac{\partial}{\partial x_1}\left[r_1x_1\left(1-\frac{x_1+\alpha_{12}x_2}{k_1}\right)P(x_1,x_2,t) \right] -\frac{\partial}{\partial x_2}\left[r_2x_2\left(1-\frac{x_2+\alpha_{21}x_1}{k_2}\right)P(x_1,x_2,t) \right] +\frac 12 \beta_1^2 \frac{\partial^2x_1^2P(x_1,x_2,t)}{\partial x_1^2} +\frac 12 \beta_2^2 \frac{\partial^2x_2^2P(x_1,x_2,t)}{\partial x_2^2}
</math>
</center>

Resolvendo as derivadas por regra do produto e separando nas derivadas parciais de <math> P(x_1,x_2,t) </math> ficamos com
<center>
<math>

\frac{\partial P(x_1,x_2,t)}{\partial t } = \delta_1 \frac{\partial^2 P}{\partial x_1^2} + \delta_2 \frac{\partial^2 P}{\partial x_2^2} + \gamma_1 \frac{\partial P}{\partial x_1} + \gamma_2 \frac{\partial P}{\partial x_2} + \omega P

</math>
</center>

com

<center>
<math>
\delta_1 = \frac{1}{2} \beta_1^2 x_1^2
</math>

<math>
\delta_2 = \frac{1}{2} \beta_2^2 x_2^2
</math>

<math>
\gamma_1 = \frac{r_1 x_1^2}{k_1} + \frac{r_1 x_1 \alpha_{12} x_2}{k_1} - r_1x_1 + 2\beta_1^2x_1
</math>

<math>
\gamma_2 = \frac{r_2 x_2^2}{k_2} + \frac{r_2 x_2 \alpha_{21} x_1}{k_2} - r_2x_2 + 2\beta_2^2x_2
</math>

<math>
\omega = \frac{2r_1x_1}{k_1} + \frac{2r_2x_2}{k_2} + \frac{r_1\alpha_{12}x_2}{k_1} + \frac{r_2\alpha_{21}x_1}{k_2} - r_1 - r_2 + \beta_1^2 + \beta_2^2
</math>

</center>

Discretizando a equação diferencial acima, usando que

<center>
<math>
\frac{\partial P(x,t)}{\partial t} = \frac{P_{i,j}^{n+1} - P_{i,j}^{n}}{dt}
</math>

<math>
\frac{\partial P(x,t)}{\partial x_1} = \frac{1}{2} \frac{P_{i+1,j}^{n} - P_{i-1,j}^{n}}{dx_1}
</math>

<math>
\frac{\partial P(x,t)}{\partial x_2} = \frac{1}{2} \frac{P_{i,j+1}^{n} - P_{i,j-1}^{n}}{dx_2}
</math>

<math>
\frac{\partial^2 P(x,t)}{\partial x_1^2} = \frac{P_{i+1,j}^{n} + P_{i-1,j}^{n} - 2P_{i,j}^{n}}{dx_1^2}
</math>

<math>
\frac{\partial^2 P(x,t)}{\partial x_2^2} = \frac{P_{i,j+1}^{n} + P_{i,j-1}^{n} - 2P_{i,j}^{n}}{dx_2^2}
</math>

</center>

e que <math> dx_1 = dx_2 </math>, chegamos, em fim, em

<center>
<math>
P_{i,j}^{n+1} = P_{i,j}^{n}\left( 1 - \frac{2\delta_1 dt}{dx^2} - \frac{2\delta_2 dt}{dx^2} + \omega dt \right) + P_{i+1,j}^{n} \left( \frac{\delta_1 dt}{dx^2} + \frac{\gamma_1 dt}{2dx} \right) + P_{i,j+1}^{n} \left( \frac{\delta_2 dt}{dx^2} + \frac{\gamma_2 dt}{2dx} \right) + P_{i-1,j}^{n} \left( \frac{\delta_1 dt}{dx^2} - \frac{\gamma_1 dt}{2dx} \right) + P_{i,j-1}^{n} \left( \frac{\delta_2 dt}{dx^2} - \frac{\gamma_2 dt}{2dx} \right)
</math>
</center>

que foi resolvida numericamente.

Abaixo, temos um gráfico considerando duas populações, com constantes todas iguais e baixo ruído. O eixo das abcissas indica a população <math>x_1</math>, o eixo das ordenadas, a população <math>x_2</math> e a cor representa a probabilidade da população possuir aquele número de indivíduos.

<center>
[[Arquivo:Tmp.gif]]
</center>

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T16:48:55Z

Gustavobpasquali: /* Método de Fokker-Planck */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Método de Fokker-Planck==
A Equação de Fokker-Planck trata da equação diferencial da distribuição de probabilidade das variáveis de um sistema de Equações Diferenciais Estocásticas. Avaliaremos sua aplicação no sistema Lotka-Volterra competitivo para duas populações.

Para um processo estocástico vetorial <math> X = [X_1, X_2, ... , X_n] </math>, a equação de Fokker-Planck tem a forma

<center>
<math>
\frac{\partial P(x,t)}{\partial t} = \sum^n_{i=1} \frac{\partial[K_i(x,t)P(x,t)]}{\partial x_i} + \frac{1}{2} \sum^n_{i,j = 1}\frac{\partial^2[D_{ij}(x,t)P(x,t)]}{\partial x_i\partial x_j}
</math>
</center>

Onde

<center>
<math>
K_i(x,t) = -A_i(x,t)
</math>
</center>
e
<center>
<math>
D_{ij}(x,t) = \sum^m_{\alpha=1}B_{i,\alpha}(x,t)B_{j,\alpha}(x,t)
</math>
</center>

Assim, para o nosso caso, temos:
<center>
<math>
\frac{\partial P(x_1,x_2,t)}{\partial t} = -\frac{\partial}{\partial x_1}\left[r_1x_1\left(1-\frac{x_1+a_{12}x_2}{k_1}\right)P(x_1,x_2,t) \right] -\frac{\partial}{\partial x_2}\left[r_2x_2\left(1-\frac{x_2+a_{21}x_1}{k_2}\right)P(x_1,x_2,t) \right] +\frac 12 \beta_1^2 \frac{\partial^2x_1^2P(x_1,x_2,t)}{\partial x_1^2} +\frac 12 \beta_2^2 \frac{\partial^2x_2^2P(x_1,x_2,t)}{\partial x_2^2}
</math>
</center>

Resolvendo as derivadas por regra do produto e separando as equações parciais de <math> P(x_1,x_2,t) </math>

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T16:44:15Z

Gustavobpasquali: /* Método de Fokker-Planck */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Método de Fokker-Planck==
A Equação de Fokker-Planck trata da equação diferencial da distribuição de probabilidade das variáveis de um sistema de Equações Diferenciais Estocásticas. Avaliaremos sua aplicação no sistema Lotka-Volterra competitivo para duas populações.

Para um processo estocástico vetorial <math> X = [X_1, X_2, ... , X_n] </math>, a equação de Fokker-Planck tem a forma

<center>
<math>
\frac{\partial P(x,t)}{\partial t} = \sum^n_{i=1} \frac{\partial[K_i(x,t)P(x,t)]}{\partial x_i} + \frac{1}{2} \sum^n_{i,j = 1}\frac{\partial^2[D_{ij}(x,t)P(x,t)]}{\partial x_i\partial x_j}
</math>
</center>

Onde

<center>
<math>
K_i(x,t) = -A_i(x,t)
</math>
</center>
e
<center>
<math>
D_{ij}(x,t) = \sum^m_{\alpha=1}B_{i,\alpha}(x,t)B_{j,\alpha}(x,t)
</math>
</center>

Assim, para o nosso caso, temos:
<center>
<math>
\frac{\partial P(x_1,x_2,t)}{\partial t} = -\frac{\partial}{\partial x_1}\left[r_1x_1\left(1-\frac{x_1+a_{12}x_2}{k_1}\right)P(x_1,x_2,t) \right] -\frac{\partial}{\partial x_2}\left[r_2x_2\left(1-\frac{x_2+a_{21}x_1}{k_2}\right)P(x_1,x_2,t) \right] +\frac 12 \beta_1^2 \frac{\partial^2x_1^2P(x_1,x_2,t)}{\partial x_1^2} +\frac 12 \beta_2^2 \frac{\partial^2x_2^2P(x_1,x_2,t)}{\partial x_2^2}
</math>
</center>

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T16:29:54Z

Gustavobpasquali: /* Método de Fokker-Planck */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Método de Fokker-Planck==
A Equação de Fokker-Planck trata da equação diferencial da distribuição de probabilidade das variáveis de um sistema de Equações Diferenciais Estocásticas. Avaliaremos sua aplicação no sistema Lotka-Volterra competitivo para duas populações.

Para um processo estocástico vetorial <math> X = [X_1, X_2, ... , X_n] </math>, a equação de Fokker-Planck tem a forma

<center>
<math>
\frac{\partial P}{\partial t} = \sum^n_{i=1} \frac{\partial[K_i(x,t)P(x,t)]}{\partial x_i} + \frac{1}{2} \sum^n_{i,j = 1}\frac{\partial^2[D_{ij}(x,t)P(x,t)]}{\partial x_i\partial x_j}
</math>
</center>

Onde

<center>
<math>
K_i(x,t) = -A_i(x,t)
</math>
</center>
e
<center>
<math>
D_{ij}(x,t) = \sum^m_{\alpha=1}B_{i,\alpha}(x,t)B_{j,\alpha}(x,t)
</math>
</center>

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T16:13:58Z

Gustavobpasquali:

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Método de Fokker-Planck==

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T16:12:53Z

Gustavobpasquali:

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Consideraremos duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

== Método de Fokker-Planck ==

=== Resultados ===

==== Caso Um ====

Com duas populações, vamos exemplificar nesta seção o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

==== Caso Dois ====

Agora, consideraremos uma pequena variação nas constantes de crescimento das populações, com um aumento de duas vezes para <math>r_1</math>. Todos os outros parâmetros são iguais ao Caso Um. Os resultados seguem nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes r1 eh 2.png]]
</center>
<center>
[[Arquivo:Duas populacoes r1 eh 2 fases.png]]
</center>

Mesmo com a população um possuindo taxa de crescimento duas vezes maior, a espécie com maior amplitude de ruído sobrevive ao final da simulação.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right) + x_3 \beta{_3} \xi{_3} (t)
</math>
</center>

=== Resultados ===

==== Caso Um ====

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

==== Caso Dois ====

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

Vemos que, mesmo considerando as relações de dominância, a população com maior ruído sobressai no final, sendo a única espécie viva ao final da simulação.

 

== Referências ==

# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T14:58:23Z

Gustavobpasquali: /* Resolução Numérica pelo Método de Itô */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Dividiremos, para tanto, o trabalho em três partes principais, considerando duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase, e generalizando para N populações. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X_1(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X_1(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X_1}B(X_1(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

Com duas populações, vamos exemplificar o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

Outras simulações foram feitas considerando diversas constantes, porém o comportamento com o ruído era praticamente o mesmo daquele sem ruído.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right)</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right)</math>
</center>

=== Resultados ===

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

== Equações para N Populações ==

<center>
<math>
\frac{dx_i}{dt} = r_ix_i\left(1 - \frac{\sum_{j=1}^{N}\alpha_{ij}x_j}{K_i}\right)
</math>
</center>

 

</source>

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T14:56:42Z

Gustavobpasquali: /* Resolução Numérica pelo Método de Itô */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Dividiremos, para tanto, o trabalho em três partes principais, considerando duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase, e generalizando para N populações. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de Langevin

<center>
<math>
A(X(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

Com duas populações, vamos exemplificar o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

Outras simulações foram feitas considerando diversas constantes, porém o comportamento com o ruído era praticamente o mesmo daquele sem ruído.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right)</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right)</math>
</center>

=== Resultados ===

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

== Equações para N Populações ==

<center>
<math>
\frac{dx_i}{dt} = r_ix_i\left(1 - \frac{\sum_{j=1}^{N}\alpha_{ij}x_j}{K_i}\right)
</math>
</center>

 

</source>

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T14:56:28Z

Gustavobpasquali: /* Resolução Numérica pelo Método de Itô */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Dividiremos, para tanto, o trabalho em três partes principais, considerando duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase, e generalizando para N populações. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para <math> x_1 </math>

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de langevin

<center>
<math>
A(X(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

Com duas populações, vamos exemplificar o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

Outras simulações foram feitas considerando diversas constantes, porém o comportamento com o ruído era praticamente o mesmo daquele sem ruído.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right)</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right)</math>
</center>

=== Resultados ===

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

== Equações para N Populações ==

<center>
<math>
\frac{dx_i}{dt} = r_ix_i\left(1 - \frac{\sum_{j=1}^{N}\alpha_{ij}x_j}{K_i}\right)
</math>
</center>

 

</source>

Equação de Lotka-Volterra Competitiva Estocástica

2024-08-27T14:55:42Z

Gustavobpasquali: /* Resolução Numérica pelo Método de Itô */

== Introdução ==

As Equações de Lotka-Volterra fornecem um modelo para a previsão de sistemas biológicos considerando diversas relações entre populações. Exploraremos no vigente trabalho a relação de competitividade, juntamente da adição de um termo estocástico multiplicativo. Dividiremos, para tanto, o trabalho em três partes principais, considerando duas e três populações, mostrando os gráficos de evolução temporal do número de indivíduos de cada espécie e os espaços de fase, e generalizando para N populações. Em todas as etapas serão mostrados os resultados considerando e desconsiderando o ruído para fins comparativos.

== Equações para Duas Populações ==

Como uma primeira análise, considerando somente duas populações distintas, vamos explorar o modelo logístico, partindo das Equações de Lotka-Volterra, de duas espécies disputando um território, que pode ser descrito pelo seguinte par de relações:
<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right)
</math>
</center>
Nesse par, <math>x_1</math> e <math>x_2</math> representam as duas populações consideradas; <math>r_1</math> e <math>r_2</math> indicam o crescimento inerente per-capita, a capacidade de uma espécie em se reproduzir; <math>K_1</math> e <math>K_2</math> retratam a capacidade de carga, o número de indivíduos limite que o meio ambiente consegue suportar considerando o nicho ecológico ao que a espécie pertence; e <math>\alpha_{12}</math> e <math>\alpha_{21}</math>, o efeito que a espécie um tem na espécie dois e vice-versa.

Para a adição de ruído ao sistema, vamos supor que o termo estocástico seja incorporado em cada equação de maneira a afetar somente a população <math> x_1 </math> ou a população <math> x_2 </math>, ou seja, estamos considerando que o ruído atinja cada população de maneira proporcional a sua quantidade de indivíduos em determinado período de tempo. Como um exemplo prático, podemos pensar no termo estocástico como a presença de um agente externo que infecte uma das populações isoladamente. Claro que o efeito resultante de existir uma doença se espalhando na espécie <math> x_1 </math> afeta o crescimento da espécie <math> x_2 </math>, entretanto esse efeito não é direto, sendo uma consequência da correlação entre o par de equações.

Além disso, antes de de fato colocarmos o fator estocástico nas equações, precisamos de explicar qual o tipo de ruído a ser inserido. Para tanto, será introduzido o ruído branco, que é definido a partir de suas propriedades estatísticas que configuram sua média e sua função de correlação temporal, podendo ser expressas matematicamente por

<center>
<math>
<\xi (t) > = 0
</math>
</center>

<center>
<math>
<\xi (t_1)\xi (t_2) > = \delta (t_2 - t_1)
</math>
</center>

A primeira relação indica que o ruído branco possui média igual a zero, enquanto que a segunda indica o fato de os ruídos considerando diferentes tempos são independentes e, sendo assim, dá o nome branco ao ruído, pois a intensidade espectral de um processo estocástico é a Transformada de Fourier da função de correlação, sendo ela constante por se tratar de uma Delta de Dirac, logo todas as frequências estão presentes com mesma intensidade, com branco fazendo uma analogia a essa característica.

Portanto, podemos reescrever as relações supracitadas considerando o ruído como

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2}\right) + x_2 \beta{_2} \xi{_2} (t)
</math>
</center>

com <math> \beta_1 </math> e <math> \beta_2 </math> sendo as amplitudes dos ruídos e <math> \xi_1 (t) </math> e <math> \xi_2 (t) </math> os ruídos brancos.

=== Resolução Numérica pelo Método de Itô ===

Para então resolvermos as equações numericamente, vamos recorrer ao cálculo diferencial estocástico. Iniciamos, para esse fim, tomando a forma de Langevin

<center>
<math>
dX = A(X(t), t)dt + B(X(t), t)\xi (t)dt
</math>
</center>

Deveríamos, como próxima etapa, integrar ambos os lados da equação da seguinte forma

<center>
<math>
\int_X^{X + \Delta X}dX' = \int_t^{t + \Delta t}A(X(t'), t')dt' + \int_t^{t + \Delta t}B(X(t'), t')\xi (t')dt'
</math>
</center>

Como não podemos resolvê-la pelos métodos convencionais, pois <math> \xi (t) </math> torna a função não diferenciável apesar de ser contínua, utilizamos a relação

<center>
<math>
B(X(t'),t') = B\left(\frac{X(t') + X(t' + \Delta t')}{2}, t'\right)
</math>
</center>

Essa troca se justifica pelo fato de estarmos trocando o argumento de B pela média entre os valores de X no início e no fim do intervalo de integração, correspondendo a uma melhor estimação. Agora, tomando a integral, ficamos com a aproximação

<center>
<math>
\Delta X(t) = A(X(t), t) \Delta t + B\left(X(t) + \frac{1}{2}\Delta X(t), t\right)\Delta W(t)
</math>
</center>

Por seguinte, expandimos B em série de Taylor, excluímos o termo de ordem maior que dt e fazemos a substituição <math> \Delta t </math> por <math> \Delta W^2 </math>

<center>
<math>
dX(t) = \underbrace{\left(A(X(t), t) + \frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t)\right)}_{A^{(I)}(X(t), t)} dt + B(X(t), t)\bullet dW(t)
</math>
</center>

Com isso, temos que a Equação diferencial de Itô é dada por

<center>
<math>
dX(t) = A^{(I)}(X(t), t) dt + B(X(t), t) \bullet dW(t)
</math>
</center>

Portanto, partindo da nossa relação para x1

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) + x_1 \beta{_1} \xi{_1} (t)
</math>
</center>

Temos, isolando <math> dx_1 </math>

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right) dt + x_1 \beta{_1} \underbrace{\xi{_1} (t) dt}_{dW}
</math>
</center>

A partir disso, sabemos os coeficientes A(X(t),t) e B(X(t),t) de langevin

<center>
<math>
A(X(t),t) = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1}\right)
</math>
</center>

<center>
<math>
B(X(t),t) = x_1 \beta{_1}
</math>
</center>

<center>
<math>
\frac{1}{2}\frac{\partial B}{\partial X}B(X(t), t) = \frac{1}{2}\beta_1^2x_1
</math>
</center>

Assim, substituindo os valores para encontrar a equação diferencial de Itô, obtemos

<center>
<math>
dx_1 = r_1x_1\left(1 - \frac{x_1 + \alpha_{12}x_2}{K_1} + \frac{1}{2}\beta_1^2x_1\right) dt + x_1 \beta{_1} dW
</math>
</center>

Analogamente, obtemos também a equação para a variável <math> x_2 </math>

<center>
<math>
dx_2 = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1}{K_2} + \frac{1}{2}\beta_2^2x_2\right) dt + x_2 \beta{_2} dW
</math>
</center>

=== Resultados ===

Com duas populações, vamos exemplificar o caso em que a competitividade se dá de maneira igual, ou seja, todas as constantes possuem mesmo valor independentemente da população. A diferença, então, se dará a partir do ruído, que será maior para uma delas.

Utilizando <math>r_1</math> = <math>r_2</math> = <math>\alpha_{12}</math> = <math>\alpha_{21}</math> = 1, <math>K_1</math> = <math>K_2</math> = 60 e população inicial de 10 indivíduos, a evolução temporal e o espaço de fases do sistema de duas populações podem ser vistos nas figuras abaixo.

<center>
[[Arquivo:Duas populacoes.png]]
</center>
<center>
[[Arquivo:Duas populacoes fases.png]]
</center>

Podemos verificar que, a partir da adição de um ruído maior na população 2, a população 1 cedeu quando o tempo atingiu valor próximo a 30 unidades, comportamento este que não é visto no gráfico sem ruído, no qual as duas populações convivem de maneira igual e atingindo cada uma metade de sua população limite. Realizaram-se mais de 20 simulações e todas tiveram o mesmo efeito, com a população 2 prosperando em consequência do desaparecimento da população 1.

Outras simulações foram feitas considerando diversas constantes, porém o comportamento com o ruído era praticamente o mesmo daquele sem ruído.

== Equações para Três Populações ==

Os métodos de resolução numérica para três populações são completamente análogos aos métodos para duas populações; portanto, a matemática envolvida será omitida a partir de agora e somente serão discutidos os resultados. O conjunto de equações utilizado para três espécies distintas é

<center>
<math>
\frac{dx_1}{dt} = r_1x_1\left(1 -\frac{x_1 + \alpha_{12}x_2 + \alpha_{13}x_3}{K_1} \right)
</math>
</center>
<center>
<math>
\frac{dx_2}{dt} = r_2x_2\left(1 - \frac{x_2 + \alpha_{21}x_1 + \alpha_{23}x_3}{K_2}\right)</math>
</center>
<center>
<math>
\frac{dx_3}{dt} = r_3x_3\left(1 - \frac{x_3 + \alpha_{31}x_1 + \alpha_{32}x_2}{K_1}\right)</math>
</center>

=== Resultados ===

Como primeira análise, vamos ao caso análogo ao anterior de duas populações. Vamos considerar espécies semelhantes, trocando entre elas somente o fator de amplitude do ruído, com o maior fator sendo da população 2. Utilizamos <math>K_1</math> = <math>K_2</math> = <math>K_3</math> = 60, população inicial de 10 indivíduos e as outras constantes como um. As figuras a seguir mostram a evolução temporal e o espaço de fases do sistema.

<center>
[[Arquivo:Tres populacoes ctes iguais.png]]
</center>
<center>
[[Arquivo:Tres populacoes ctes iguais fases.png]]
</center>

Podemos verificar que, assim como na seção anterior, a população com maior fator de ruído prosperou enquanto que as outras tiveram sua extinção em um tempo de 5 unidades. Quando desconsideramos o ruído, vemos que as três espécies convivem com mesmo comportamento populacional.

Para um segundo caso, vamos considerar os mesmos parâmetros de antes, modificando apenas os coeficientes entre as populações, dados na seguinte tabela.

<center>
{| class="wikitable"
|+ Coeficientes <math>\alpha</math>
|-
!
! <math>\alpha_{i1}</math>
! <math>\alpha_{i2}</math>
! <math>\alpha_{i3}</math>
|-
! <math>\alpha_{1j}</math>
| 1
| 2
| 1
|-
! <math>\alpha_{2j}</math>
| 1
| 1
| 2
|-
! <math>\alpha_{3j}</math>
| 2
| 1
| 1
|}
</center>

Eles indicam que a população 1 tem dominância sobre a população 2, a população 2 tem dominância sobre a população 3 e a população 3 tem dominância sobre a população 1. As imagens podem ser vistas a seguir.

<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1.png]]
</center>
<center>
[[Arquivo:Tres populacoes 1 2 1 1 1 2 2 1 1 fases.png]]
</center>

== Equações para N Populações ==

<center>
<math>
\frac{dx_i}{dt} = r_ix_i\left(1 - \frac{\sum_{j=1}^{N}\alpha_{ij}x_j}{K_i}\right)
</math>
</center>

 

</source>

Equação de Schrödinger Unidimensional

2024-04-27T20:25:01Z

Gustavobpasquali: /* Estabilidade do Método de Crank-Nicolson */

Conhecer o estado de uma partícula, na Mecânica Clássica, é saber, em um determinado momento, sua velocidade e sua posição. A partir destes valores, é possível determinar os seus futuros estados utilizando as relações fornecidas pelas Leis de Newton. Na Mecânica Quântica, porém, o estado físico não é mais caracterizado por uma quantidade discreta de valores numéricos, e sim, por uma função. A chamada Função de Onda, <math> \Psi </math>, nos informa a evolução temporal do estado quântico de uma partícula, seguindo a Equação de Schrödinger. Em uma dimensão, a Equação de Schrödinger é dada por:

<math> -\frac{\hbar}{2m}\frac{\partial^2\Psi}{\partial x^2} + V(x)\Psi = i\hbar\frac{\partial\Psi}{\partial t} </math>

== O Método De Crank-Nicolson==

A fim de resolvê-la numericamente para diferentes tipos de potenciais, usou-se o Método de Crank-Nicolson de resolução de EDPs. Este método baseia-se em combinar os métodos explícito e implícito (com igual contribuição) para aumentar sua estabilidade.

O método explícito FTCS (Forward Time Central Space) consiste em tomar derivadas temporais "para frente" e manter a derivada espacial centrada.

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t} </math>

<math> \frac{\partial\Psi(x,t)}{\partial x} = \frac{\Psi(x+\Delta x,t) - \Psi(x,t)}{\Delta x} </math>; <math> \frac{\partial^2\Psi(x,t)}{\partial^2 x} = \frac{(\Psi(x + \Delta x,t) - \Psi(x+\Delta x -\Delta x,t)) - ( \Psi(x,t) - \Psi(x - \Delta x,t) )}{\Delta x \Delta x} = \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} </math>

Tomando <math> \hbar=m=1 </math> (unidades atômicas), ao inserirmos as expressões acima na Equação de Schrödinger, temos que:

<math> \frac{i}{2} \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} -i V(x)\Psi(x,t) = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t}</math>

Simplificando a notação para <math> \Psi(x,t)=\Psi_j^n </math>, onde <math> j </math> representa a posição e <math> n </math> o tempo, e reorganizando os termos, ficamos com:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

O método implícito é semelhante, porém a derivada temporal é "para trás".

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t) - \Psi(x,t-\Delta t)}{\Delta t} </math>

O que nos leva a ter:

<math> \Psi_j^{n} = \Psi_{j}^{n-1} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

que, com uma substituição de variável (que é muda), nos dá, enfim:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - i\Delta t V_j \Psi_{j}^{n+1}</math>

Como dito, o Método de Crank-Nicolson é uma combinação dos métodos implícito e explícito, cada um tendo o mesmo peso de 1/2. Ao somarmos as duas equações (cada uma multiplicada por 0.5), ficaremos com a seguinte expressão:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n} + \Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - \frac{i\Delta t}{2} V_j (\Psi_{j}^{n} + \Psi_{j}^{n+1})</math>

que, ao reorganizarmos para isolar os <math> \Psi_j^{n+1}</math>, resulta em:

<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>

Para resolver a equação em simultaneamente para ambos os tempos e em todos os pontos, escrevemos a equação em forma matricial:
:<math>
\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{j}^{n+1}\end{pmatrix}
=
\begin{pmatrix}
1 - 2b - \frac{i\Delta t V_j}{2} & +b & 0 & \cdots & 0 & 0 & 0\\ +b & 1 - 2b - \frac{i\Delta t V_j}{2} & +b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & +b & 1 - 2b - \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{j}^{n} \end{pmatrix} ;
</math>
com
<math> b=\frac{i \Delta t}{4 \Delta x^2}</math>.

Podemos notar que a diferença entre as matrizes que multiplicam <math> \Psi_j^{n+1}</math> e <math> \Psi_j^{n}</math> é somente o sinal das componente imaginárias. Assim, a equação final pode ser escrita como:

<math> A\Psi_j^{n+1} = A^*\Psi_j^{n}</math>, onde A é a matriz <math>\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix} </math>

===Estabilidade do Método de Crank-Nicolson===
Para verificar a estabilidade, devemos primeiro olhar para a equação que define o método:
<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>
Para simplificar, chamaremos <math>\frac{i \Delta t}{4 \Delta x^2}</math> de b. Nesse caso, temos:
<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - b\Psi_{j+1}^{n+1} - b\Psi_{j-1}^{n+1} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + b\Psi_{j+1}^{n} +b\Psi_{j-1}^{n}</math>
Definimos os Modos de Fourier:

<math>\Psi_{j}^{n} \approx A^{n}e^{iqjh}</math>

Assim, obtemos:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}e^{iqjh}- bA^{n+1}e^{iq(j-1)h} - bA^{n+1}e^{iq(j+1)h} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n}e^{iqjh} + bA^{n}e^{iq(j-1)h} +bA^{n}e^{iq(j+1)h}</math>

Dividindo tudo por <math>e^{iqjh}</math>:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}- bA^{n+1}e^{-iqh} - bA^{n+1}e^{iqh} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n} + bA^{n}e^{-iqh} +bA^{n}e^{iqh}</math>

Escrito de outra forma:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(e^{-iqh}-2+e^{iqh})\right] = A^{n}\left[1-\frac{i\Delta t}{2} V_j +b(e^{-iqh}-2+e^{iqh})\right] </math>

Utilizando relações trigonométricas:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(-sin^2(\frac{qh}2))\right] = A^{n}\left[1-\frac{i\Delta t}{2} V_j +b(-sin^2(\frac{qh}2))\right]</math>

Assim, obtemos o fator de amplificação:

<math>\left|\frac{A^{n+1}}{A^n}\right| = \left|\frac{1-\frac{i\Delta t}{2} V_j -b(sin^2(\frac{qh}2))}{1+\frac{i\Delta t}{2} V_j +b(sin^2(\frac{qh}2))}\right| \leq 1</math>

Portanto, o método é incondicionalmente estável.

== Exemplos de Potenciais ==

=== Oscilador Harmônico ===

O oscilador harmônico quântico consiste de uma partícula em uma zona onde o potencial possui comportamento harmônico, ou seja, depende de <math> x^2 </math>. Em nosso exemplo, inserimos um pacote de formato gaussiano no interior de um potencial harmônico e observamos sua evolução temporal.

[[Arquivo:H_oscilator.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

# Definição de constantes
ħ = 1
xmax = 100
tmax = 150
m = 1
Δt = 0.2
Δx = 0.25
σ = 2
k = 75
ω = 0.1
x0 = xmax/2

# Definição do espaço de integração e do pacote de onda inicial
x = collect(0:Δx:xmax)
Ψ_HO = zeros(Complex{Float64}, length(x))
for (i, val) in enumerate(x)
Ψ_HO[i] = 0.5*exp(-(val - xmax/3)^2/(2*σ^2))*exp(im*k*(val-xmax/3))
end

# Potencial do oscilador harmônico
function V_HO(x)
return 0.5 * m * (ω^2) * ((x - x0)^2)
end

# Pré-cálculo do potencial do oscilador harmônico
V_HO_values = [V_HO(x_i) for x_i in x[2:end-1]]

# Matrizes de evolução temporal
A_HO = zeros(Complex{Float64}, length(x)-2, length(x)-2)
for i in 1:length(x)-2, j in 1:length(x)-2
if i == j
A_HO[i, j] = 1 + im * (Δt / (2ħ)) * (ħ^2 / (m * Δx^2) + V_HO_values[i])
elseif i == j + 1 || i == j - 1
A_HO[i, j] = -im * (ħ * Δt) / (4 * m * Δx^2)
end
end

invA_HO = inv(A_HO)
B_HO = conj(A_HO)

# Criação do GIF
anim = @animate for t in 0:Δt:tmax
Ψ_HO[2:end-1] .= invA_HO * (B_HO * Ψ_HO[2:end-1])
plot(x, abs2.(Ψ_HO), xlabel="Position", title="Time: $t", ylims=(0, 1), label="|Ψ|²")
plot!(x, V_HO.(x), label="Potential V(x)")
end

gif_path = "/gif/path/harmonic_oscilator.gif"
gif(anim, gif_path, fps=45)

</source> 

=== Barreira (Tunelamento) ===

Na mecânica quântica partículas podem existir, mesmo que com baixa probabilidade, em regiões que classicamente seriam proibidas. Este é o caso da barreira, onde uma partícula atravessa uma região onde sua energia total é negativa. Chamamos este fenômeno de tunelamento. A seguir, há um exemplo onde inserimos um pacote de formato gaussiano se deslocando em direção a uma barreira de potencial, de forma que parte da onda é refletida de volta e outra tunela pela barreira.

[[Arquivo:Tunneling.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

# Definição de constantes
ħ = 1
xmax = 500
tmax = 200
m = 1
Δt = 0.2
Δx = 0.25
σ = 1.5
k = 50
ω = 0.1

# Definição do espaço de integração e do pacote de onda inicial
x = collect(0:Δx:xmax)
Ψ = zeros(Complex{Float64}, length(x))
for (i, val) in enumerate(x)
Ψ[i] = Ψ[i] = 0.5*exp(-(val - xmax/2.5)^2/(2*σ^2))*exp(-im*k*(val))
end

# Potencial barreira
function V(arg)
if arg > xmax*0.5 && arg <= xmax*0.51
return 0.1
else
return 0
end
end

# Pré-cálculo do potencial V
V_values = [V(x_i) for x_i in x[2:end-1]]

# Matrizes de evolução temporal A
A = zeros(Complex{Float64}, length(x)-2, length(x)-2)
for i in 1:length(x)-2, j in 1:length(x)-2
if i == j
A[i, j] = 1 + im * (Δt / (2ħ)) * (ħ^2 / (m * Δx^2) + V_values[i])
elseif i == j + 1 || i == j - 1
A[i, j] = -im * (ħ * Δt) / (4 * m * Δx^2)
end
end

invA = inv(A)
B = conj(A)

# Criação do GIF
anim = @animate for t in 0:Δt:tmax
global Ψ[2:length(x)-1] = invA * (B * Ψ[2:length(x)-1])
plot(x, abs2.(Ψ), xlabel="Position", title="Time: $t", ylims=(0, 0.1), label="|Ψ|²")
plot!(x, V.(x), label="Potential V(x)")
end every 10

gif_path = "/gif/path/tunneling.gif"
gif(anim, gif_path, fps = 30)
</source> 

=== Poço Finito ===

O contrário da barreira é o chamado de poço de potencial. Ao invés da partícula atravessá-lo por completo, uma parte do pacote de onda é refletido, uma parte é transmitido e outra ainda fica preso no meio do poço. Veja o exemplo.

[[Arquivo:Pit.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_p(x)
if x<L/3
return 0.1
elseif x>2L/3
return 0.1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_p = [if i==j; 1 - 2b - V_p(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_p = @. conj(A_p)
IB_p = inv(B_p)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/4)^2/(8))*exp(-im*75*(val-L/4))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_p = @. V_p(x)

t=0
@gif for t in 0:dt:50
ψ[2:end-1] = IB_p[2:end-1,2:end-1]*A_p[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_p,c="blue",label="V(x)")
plot!(ylim=[0,0.15],title="t=$t")
end

</source> 

== Referências ==

# Griffiths, DAVID. Introduction to Quantum Mechanics 3ed. 2018.
# Scherer, CLÁUDIO. Métodos Computacionais da Física. 2010.

Equação de Schrödinger Unidimensional

2024-04-23T22:53:14Z

Gustavobpasquali: /* Estabilidade do Método de Crank-Nicolson */

Conhecer o estado de uma partícula, na Mecânica Clássica, é saber, em um determinado momento, sua velocidade e sua posição. A partir destes valores, é possível determinar os seus futuros estados utilizando as relações fornecidas pelas Leis de Newton. Na Mecânica Quântica, porém, o estado físico não é mais caracterizado por uma quantidade discreta de valores numéricos, e sim, por uma função. A chamada Função de Onda, <math> \Psi </math>, nos informa a evolução temporal do estado quântico de uma partícula, seguindo a Equação de Schrödinger. Em uma dimensão, a Equação de Schrödinger é dada por:

<math> -\frac{\hbar}{2m}\frac{\partial^2\Psi}{\partial x^2} + V(x)\Psi = i\hbar\frac{\partial\Psi}{\partial t} </math>

== O Método De Crank-Nicolson==

A fim de resolvê-la numericamente para diferentes tipos de potenciais, usou-se o Método de Crank-Nicolson de resolução de EDPs. Este método baseia-se em combinar os métodos explícito e implícito (com igual contribuição) para aumentar sua estabilidade.

O método explícito FTCS (Forward Time Central Space) consiste em tomar derivadas temporais "para frente" e manter a derivada espacial centrada.

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t} </math>

<math> \frac{\partial\Psi(x,t)}{\partial x} = \frac{\Psi(x+\Delta x,t) - \Psi(x,t)}{\Delta x} </math>; <math> \frac{\partial^2\Psi(x,t)}{\partial^2 x} = \frac{(\Psi(x + \Delta x,t) - \Psi(x+\Delta x -\Delta x,t)) - ( \Psi(x,t) - \Psi(x - \Delta x,t) )}{\Delta x \Delta x} = \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} </math>

Tomando <math> \hbar=m=1 </math> (unidades atômicas), ao inserirmos as expressões acima na Equação de Schrödinger, temos que:

<math> \frac{i}{2} \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} -i V(x)\Psi(x,t) = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t}</math>

Simplificando a notação para <math> \Psi(x,t)=\Psi_j^n </math>, onde <math> j </math> representa a posição e <math> n </math> o tempo, e reorganizando os termos, ficamos com:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

O método implícito é semelhante, porém a derivada temporal é "para trás".

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t) - \Psi(x,t-\Delta t)}{\Delta t} </math>

O que nos leva a ter:

<math> \Psi_j^{n} = \Psi_{j}^{n-1} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

que, com uma substituição de variável (que é muda), nos dá, enfim:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - i\Delta t V_j \Psi_{j}^{n+1}</math>

Como dito, o Método de Crank-Nicolson é uma combinação dos métodos implícito e explícito, cada um tendo o mesmo peso de 1/2. Ao somarmos as duas equações (cada uma multiplicada por 0.5), ficaremos com a seguinte expressão:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n} + \Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - \frac{i\Delta t}{2} V_j (\Psi_{j}^{n} + \Psi_{j}^{n+1})</math>

que, ao reorganizarmos para isolar os <math> \Psi_j^{n+1}</math>, resulta em:

<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>

Para resolver a equação em simultaneamente para ambos os tempos e em todos os pontos, escrevemos a equação em forma matricial:
:<math>
\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{j}^{n+1}\end{pmatrix}
=
\begin{pmatrix}
1 - 2b - \frac{i\Delta t V_j}{2} & +b & 0 & \cdots & 0 & 0 & 0\\ +b & 1 - 2b - \frac{i\Delta t V_j}{2} & +b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & +b & 1 - 2b - \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{j}^{n} \end{pmatrix} ;
</math>
com
<math> b=\frac{i \Delta t}{4 \Delta x^2}</math>.

Podemos notar que a diferença entre as matrizes que multiplicam <math> \Psi_j^{n+1}</math> e <math> \Psi_j^{n}</math> é somente o sinal das componente imaginárias. Assim, a equação final pode ser escrita como:

<math> A\Psi_j^{n+1} = A^*\Psi_j^{n}</math>, onde A é a matriz <math>\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix} </math>

===Estabilidade do Método de Crank-Nicolson===
Para verificar a estabilidade, devemos primeiro olhar para a equação que define o método:
<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>
Para simplificar, chamaremos <math>\frac{i \Delta t}{4 \Delta x^2}</math> de b. Nesse caso, temos:
<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - b\Psi_{j+1}^{n+1} - b\Psi_{j-1}^{n+1} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + b\Psi_{j+1}^{n} +b\Psi_{j-1}^{n}</math>
Definimos os Modos de Fourier:

<math>\Psi_{j}^{n} \approx A^{n}e^{iqjh}</math>

Assim, obtemos:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}e^{iqjh}- bA^{n+1}e^{iq(j-1)h} - bA^{n+1}e^{iq(j+1)h} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n}e^{iqjh} + bA^{n}e^{iq(j-1)h} +bA^{n}e^{iq(j+1)h}</math>

Dividindo tudo por <math>e^{iqjh}</math>:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}- bA^{n+1}e^{-iqh} - bA^{n+1}e^{iqh} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n} + bA^{n}e^{-iqh} +bA^{n}e^{iqh}</math>

Escrito de outra forma:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(e^{-iqh}-2+e^{iqh})\right] = A^{n}\left[1-\frac{i\Delta t}{2} V_j +b(e^{-iqh}-2+e^{iqh})\right] </math>

Utilizando relações trigonométricas:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(-sin^2(\frac{qh}2))\right] = A^{n}\left[1-\frac{i\Delta t}{2} V_j +b(-sin^2(\frac{qh}2))\right]</math>

Assim, obtemos o fator de amplificação:

<math>\frac{A^{n+1}}{A^n} = \frac{1-\frac{i\Delta t}{2} V_j -b(sin^2(\frac{qh}2))}{1+\frac{i\Delta t}{2} V_j +b(sin^2(\frac{qh}2))} \leq 1</math>

Portanto, o método é estavel.

== Exemplos de Potenciais ==

=== Oscilador Harmônico ===

O oscilador harmônico quântico consiste de uma partícula em uma zona onde o potencial possui comportamento harmônico, ou seja, depende de <math> x^2 </math>. Em nosso exemplo, inserimos um pacote de formato gaussiano no interior de um potencial harmônico e observamos sua evolução temporal.

[[Arquivo:H_oscilator.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

# Definição de constantes
ħ = 1
xmax = 100
tmax = 150
m = 1
Δt = 0.2
Δx = 0.25
σ = 2
k = 75
ω = 0.1
x0 = xmax/2

# Definição do espaço de integração e do pacote de onda inicial
x = collect(0:Δx:xmax)
Ψ_HO = zeros(Complex{Float64}, length(x))
for (i, val) in enumerate(x)
Ψ_HO[i] = 0.5*exp(-(val - xmax/3)^2/(2*σ^2))*exp(im*k*(val-xmax/3))
end

# Potencial do oscilador harmônico
function V_HO(x)
return 0.5 * m * (ω^2) * ((x - x0)^2)
end

# Pré-cálculo do potencial do oscilador harmônico
V_HO_values = [V_HO(x_i) for x_i in x[2:end-1]]

# Matrizes de evolução temporal
A_HO = zeros(Complex{Float64}, length(x)-2, length(x)-2)
for i in 1:length(x)-2, j in 1:length(x)-2
if i == j
A_HO[i, j] = 1 + im * (Δt / (2ħ)) * (ħ^2 / (m * Δx^2) + V_HO_values[i])
elseif i == j + 1 || i == j - 1
A_HO[i, j] = -im * (ħ * Δt) / (4 * m * Δx^2)
end
end

invA_HO = inv(A_HO)
B_HO = conj(A_HO)

# Criação do GIF
anim = @animate for t in 0:Δt:tmax
Ψ_HO[2:end-1] .= invA_HO * (B_HO * Ψ_HO[2:end-1])
plot(x, abs2.(Ψ_HO), xlabel="Position", title="Time: $t", ylims=(0, 1), label="|Ψ|²")
plot!(x, V_HO.(x), label="Potential V(x)")
end

gif_path = "/gif/path/harmonic_oscilator.gif"
gif(anim, gif_path, fps=45)

</source> 

=== Barreira (Tunelamento) ===

Na mecânica quântica partículas podem existir, mesmo que com baixa probabilidade, em regiões que classicamente seriam proibidas. Este é o caso da barreira, onde uma partícula atravessa uma região onde sua energia total é negativa. Chamamos este fenômeno de tunelamento. A seguir, há um exemplo onde inserimos um pacote de formato gaussiano se deslocando em direção a uma barreira de potencial, de forma que parte da onda é refletida de volta e outra tunela pela barreira.

[[Arquivo:Tunneling.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

# Definição de constantes
ħ = 1
xmax = 500
tmax = 200
m = 1
Δt = 0.2
Δx = 0.25
σ = 1.5
k = 50
ω = 0.1

# Definição do espaço de integração e do pacote de onda inicial
x = collect(0:Δx:xmax)
Ψ = zeros(Complex{Float64}, length(x))
for (i, val) in enumerate(x)
Ψ[i] = Ψ[i] = 0.5*exp(-(val - xmax/2.5)^2/(2*σ^2))*exp(-im*k*(val))
end

# Potencial barreira
function V(arg)
if arg > xmax*0.5 && arg <= xmax*0.51
return 0.1
else
return 0
end
end

# Pré-cálculo do potencial V
V_values = [V(x_i) for x_i in x[2:end-1]]

# Matrizes de evolução temporal A
A = zeros(Complex{Float64}, length(x)-2, length(x)-2)
for i in 1:length(x)-2, j in 1:length(x)-2
if i == j
A[i, j] = 1 + im * (Δt / (2ħ)) * (ħ^2 / (m * Δx^2) + V_values[i])
elseif i == j + 1 || i == j - 1
A[i, j] = -im * (ħ * Δt) / (4 * m * Δx^2)
end
end

invA = inv(A)
B = conj(A)

# Criação do GIF
anim = @animate for t in 0:Δt:tmax
global Ψ[2:length(x)-1] = invA * (B * Ψ[2:length(x)-1])
plot(x, abs2.(Ψ), xlabel="Position", title="Time: $t", ylims=(0, 0.1), label="|Ψ|²")
plot!(x, V.(x), label="Potential V(x)")
end every 10

gif_path = "/gif/path/tunneling.gif"
gif(anim, gif_path, fps = 30)
</source> 

=== Poço Finito ===

O contrário da barreira é o chamado de poço de potencial. Ao invés da partícula atravessá-lo por completo, uma parte do pacote de onda é refletido, uma parte é transmitido e outra ainda fica preso no meio do poço. Veja o exemplo.

[[Arquivo:Pit.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_p(x)
if x<L/3
return 0.1
elseif x>2L/3
return 0.1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_p = [if i==j; 1 - 2b - V_p(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_p = @. conj(A_p)
IB_p = inv(B_p)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/4)^2/(8))*exp(-im*75*(val-L/4))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_p = @. V_p(x)

t=0
@gif for t in 0:dt:50
ψ[2:end-1] = IB_p[2:end-1,2:end-1]*A_p[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_p,c="blue",label="V(x)")
plot!(ylim=[0,0.15],title="t=$t")
end

</source>

Equação de Schrödinger Unidimensional

2024-04-23T17:31:58Z

Gustavobpasquali: /* Estabilidade do Método de Crank-Nicolson */

Conhecer o estado de uma partícula, na Mecânica Clássica, é saber, em um determinado momento, sua velocidade e sua posição. A partir destes valores, é possível determinar os seus futuros estados utilizando as relações fornecidas pelas Leis de Newton. Na Mecânica Quântica, porém, o estado físico não é mais caracterizado por uma quantidade discreta de valores numéricos, e sim, por uma função. A chamada Função de Onda, <math> \Psi </math>, nos informa a evolução temporal do estado quântico de uma partícula, seguindo a Equação de Schrödinger. Em uma dimensão, a Equação de Schrödinger é dada por:

<math> -\frac{\hbar}{2m}\frac{\partial^2\Psi}{\partial x^2} + V(x)\Psi = i\hbar\frac{\partial\Psi}{\partial t} </math>

== O Método De Crank-Nicolson==

A fim de resolvê-la numericamente para diferentes tipos de potenciais, usou-se o Método de Crank-Nicolson de resolução de EDPs. Este método baseia-se em combinar os métodos explícito e implícito (com igual contribuição) para aumentar sua estabilidade.

O método explícito FTCS (Forward Time Central Space) consiste em tomar derivadas temporais "para frente" e manter a derivada espacial centrada.

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t} </math>

<math> \frac{\partial\Psi(x,t)}{\partial x} = \frac{\Psi(x+\Delta x,t) - \Psi(x,t)}{\Delta x} </math>; <math> \frac{\partial^2\Psi(x,t)}{\partial^2 x} = \frac{(\Psi(x + \Delta x,t) - \Psi(x+\Delta x -\Delta x,t)) - ( \Psi(x,t) - \Psi(x - \Delta x,t) )}{\Delta x \Delta x} = \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} </math>

Tomando <math> \hbar=m=1 </math> (unidades atômicas), ao inserirmos as expressões acima na Equação de Schrödinger, temos que:

<math> \frac{i}{2} \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} -i V(x)\Psi(x,t) = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t}</math>

Simplificando a notação para <math> \Psi(x,t)=\Psi_j^n </math>, onde <math> j </math> representa a posição e <math> n </math> o tempo, e reorganizando os termos, ficamos com:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

O método implícito é semelhante, porém a derivada temporal é "para trás".

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t) - \Psi(x,t-\Delta t)}{\Delta t} </math>

O que nos leva a ter:

<math> \Psi_j^{n} = \Psi_{j}^{n-1} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

que, com uma substituição de variável (que é muda), nos dá, enfim:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - i\Delta t V_j \Psi_{j}^{n+1}</math>

Como dito, o Método de Crank-Nicolson é uma combinação dos métodos implícito e explícito, cada um tendo o mesmo peso de 1/2. Ao somarmos as duas equações (cada uma multiplicada por 0.5), ficaremos com a seguinte expressão:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n} + \Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - \frac{i\Delta t}{2} V_j (\Psi_{j}^{n} + \Psi_{j}^{n+1})</math>

que, ao reorganizarmos para isolar os <math> \Psi_j^{n+1}</math>, resulta em:

<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>

Para resolver a equação em simultaneamente para ambos os tempos e em todos os pontos, escrevemos a equação em forma matricial:
:<math>
\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{j}^{n+1}\end{pmatrix}
=
\begin{pmatrix}
1 - 2b - \frac{i\Delta t V_j}{2} & +b & 0 & \cdots & 0 & 0 & 0\\ +b & 1 - 2b - \frac{i\Delta t V_j}{2} & +b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & +b & 1 - 2b - \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{j}^{n} \end{pmatrix} ;
</math>
com
<math> b=\frac{i \Delta t}{4 \Delta x^2}</math>.

Podemos notar que a diferença entre as matrizes que multiplicam <math> \Psi_j^{n+1}</math> e <math> \Psi_j^{n}</math> é somente o sinal das componente imaginárias. Assim, a equação final pode ser escrita como:

<math> A\Psi_j^{n+1} = A^*\Psi_j^{n}</math>, onde A é a matriz <math>\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix} </math>

==Estabilidade do Método de Crank-Nicolson==
Para verificar a estabilidade, devemos primeiro olhar para a equação que define o método:
<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>
Para simplificar, chamaremos <math>\frac{i \Delta t}{4 \Delta x^2}</math> de b. Nesse caso, temos:
<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - b\Psi_{j+1}^{n+1} - b\Psi_{j-1}^{n+1} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + b\Psi_{j+1}^{n} +b\Psi_{j-1}^{n}</math>
Definimos os Modos de Fourier:

<math>\Psi_{j}^{n} \approx A^{n}e^{iqjh}</math>

Assim, obtemos:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}e^{iqjh}- bA^{n+1}e^{iq(j-1)h} - bA^{n+1}e^{iq(j+1)h} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n}e^{iqjh} + bA^{n}e^{iq(j-1)h} +bA^{n}e^{iq(j+1)h}</math>

Dividindo tudo por <math>e^{iqjh}</math>:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}- bA^{n+1}e^{-iqh} - bA^{n+1}e^{iqh} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n} + bA^{n}e^{-iqh} +bA^{n}e^{iqh}</math>

Escrito de outra forma:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(e^{-iqh}-2+e^{iqh})\right] = A^{n}\left[1+\frac{i\Delta t}{2} V_j +b(e^{-iqh}-2+e^{iqh})\right] </math>

Utilizando relações trigonométricas:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(-sin^2(\frac{qh}2))\right] = A^{n}\left[1+\frac{i\Delta t}{2} V_j +b(-sin^2(\frac{qh}2))\right]</math>

Assim, obtemos o fator de amplificação:

<math>\frac{A^{n+1}}{A^n} = \frac{1+\frac{i\Delta t}{2} V_j -b(sin^2(\frac{qh}2))}{1+\frac{i\Delta t}{2} V_j +b(sin^2(\frac{qh}2))} \leq 1</math>

Portanto, o método é estavel.

== Exemplos de Potenciais ==

=== Oscilador Harmônico ===

O oscilador harmônico quântico consiste de uma partícula em uma zona onde o potencial possui comportamento harmônico, ou seja, depende de <math> x^2 </math>. Em nosso exemplo, inserimos um pacote de formato gaussiano no interior de um potencial harmônico e observamos sua evolução temporal.

[[Arquivo:Harmonic_oscilator.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_Oh(x)
return 0.005(x-L/2)^2
end

b = dt*im/(4*dx^2)
A = [if i==j; 1 - 2b - V_Oh(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B = @. conj(A)
IB = inv(B)

len = length(x)
ψ = [(sqrt(1/(2*pi))*exp(-(val - L/3)^2/(2))*exp(-im*75*(val-L/3))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves = @. V_Oh(x)

t=0
@gif for t in 0:dt:65
ψ[2:end-1] = IB[2:end-1,2:end-1]*A[2:end-1,2:end-1]*ψ[2:end-1]

plot(x,abs2.(ψ),c="red")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves,c="blue")
plot!(ylim=[0,0.5],title="t=$t")
end every 3

</source> 

=== Barreira (Tunelamento) ===

Na mecânica quântica partículas podem existir, mesmo que com baixa probabilidade, em regiões que classicamente seriam proibidas. Este é o caso da barreira, onde uma partícula atravessa uma região onde sua energia total é negativa. Chamamos este fenômeno de tunelamento. A seguir, há um exemplo onde inserimos um pacote de formato gaussiano se deslocando em direção a uma barreira de potencial, de forma que parte da onda é refletida de volta e outra tunela pela barreira.

[[Arquivo:Barrier3.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_b(x)
if x==L/2
return 1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_b = [if i==j; 1 - 2b - V_b(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_b = @. conj(A_b)
IB_b = inv(B_b)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/3)^2/(8))*exp(-im*75*(val-L/3))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_b = @. V_b(x)

t=0
@gif for t in 0:dt:35
ψ[2:end-1] = IB_b[2:end-1,2:end-1]*A_b[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_b,c="blue",label="V(x)")
plot!(ylim=[0,0.2],title="t=$t")
end

</source> 

=== Poço Finito ===

O contrário da barreira é o chamado de poço de potencial. Ao invés da partícula atravessá-lo por completo, uma parte do pacote de onda é refletido, uma parte é transmitido e outra ainda fica preso no meio do poço. Veja o exemplo.

[[Arquivo:Pit.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_p(x)
if x<L/3
return 0.1
elseif x>2L/3
return 0.1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_p = [if i==j; 1 - 2b - V_p(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_p = @. conj(A_p)
IB_p = inv(B_p)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/4)^2/(8))*exp(-im*75*(val-L/4))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_p = @. V_p(x)

t=0
@gif for t in 0:dt:50
ψ[2:end-1] = IB_p[2:end-1,2:end-1]*A_p[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_p,c="blue",label="V(x)")
plot!(ylim=[0,0.15],title="t=$t")
end

</source>

Equação de Schrödinger Unidimensional

2024-04-23T17:28:23Z

Gustavobpasquali: /* Estabilidade do Método de Crank-Nicolson */

Conhecer o estado de uma partícula, na Mecânica Clássica, é saber, em um determinado momento, sua velocidade e sua posição. A partir destes valores, é possível determinar os seus futuros estados utilizando as relações fornecidas pelas Leis de Newton. Na Mecânica Quântica, porém, o estado físico não é mais caracterizado por uma quantidade discreta de valores numéricos, e sim, por uma função. A chamada Função de Onda, <math> \Psi </math>, nos informa a evolução temporal do estado quântico de uma partícula, seguindo a Equação de Schrödinger. Em uma dimensão, a Equação de Schrödinger é dada por:

<math> -\frac{\hbar}{2m}\frac{\partial^2\Psi}{\partial x^2} + V(x)\Psi = i\hbar\frac{\partial\Psi}{\partial t} </math>

== O Método De Crank-Nicolson==

A fim de resolvê-la numericamente para diferentes tipos de potenciais, usou-se o Método de Crank-Nicolson de resolução de EDPs. Este método baseia-se em combinar os métodos explícito e implícito (com igual contribuição) para aumentar sua estabilidade.

O método explícito FTCS (Forward Time Central Space) consiste em tomar derivadas temporais "para frente" e manter a derivada espacial centrada.

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t} </math>

<math> \frac{\partial\Psi(x,t)}{\partial x} = \frac{\Psi(x+\Delta x,t) - \Psi(x,t)}{\Delta x} </math>; <math> \frac{\partial^2\Psi(x,t)}{\partial^2 x} = \frac{(\Psi(x + \Delta x,t) - \Psi(x+\Delta x -\Delta x,t)) - ( \Psi(x,t) - \Psi(x - \Delta x,t) )}{\Delta x \Delta x} = \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} </math>

Tomando <math> \hbar=m=1 </math> (unidades atômicas), ao inserirmos as expressões acima na Equação de Schrödinger, temos que:

<math> \frac{i}{2} \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} -i V(x)\Psi(x,t) = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t}</math>

Simplificando a notação para <math> \Psi(x,t)=\Psi_j^n </math>, onde <math> j </math> representa a posição e <math> n </math> o tempo, e reorganizando os termos, ficamos com:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

O método implícito é semelhante, porém a derivada temporal é "para trás".

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t) - \Psi(x,t-\Delta t)}{\Delta t} </math>

O que nos leva a ter:

<math> \Psi_j^{n} = \Psi_{j}^{n-1} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

que, com uma substituição de variável (que é muda), nos dá, enfim:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - i\Delta t V_j \Psi_{j}^{n+1}</math>

Como dito, o Método de Crank-Nicolson é uma combinação dos métodos implícito e explícito, cada um tendo o mesmo peso de 1/2. Ao somarmos as duas equações (cada uma multiplicada por 0.5), ficaremos com a seguinte expressão:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n} + \Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - \frac{i\Delta t}{2} V_j (\Psi_{j}^{n} + \Psi_{j}^{n+1})</math>

que, ao reorganizarmos para isolar os <math> \Psi_j^{n+1}</math>, resulta em:

<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>

Para resolver a equação em simultaneamente para ambos os tempos e em todos os pontos, escrevemos a equação em forma matricial:
:<math>
\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{j}^{n+1}\end{pmatrix}
=
\begin{pmatrix}
1 - 2b - \frac{i\Delta t V_j}{2} & +b & 0 & \cdots & 0 & 0 & 0\\ +b & 1 - 2b - \frac{i\Delta t V_j}{2} & +b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & +b & 1 - 2b - \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{j}^{n} \end{pmatrix} ;
</math>
com
<math> b=\frac{i \Delta t}{4 \Delta x^2}</math>.

Podemos notar que a diferença entre as matrizes que multiplicam <math> \Psi_j^{n+1}</math> e <math> \Psi_j^{n}</math> é somente o sinal das componente imaginárias. Assim, a equação final pode ser escrita como:

<math> A\Psi_j^{n+1} = A^*\Psi_j^{n}</math>, onde A é a matriz <math>\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix} </math>

==Estabilidade do Método de Crank-Nicolson==
Para verificar a estabilidade, devemos primeiro olhar para a equação que define o método:
<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>
Para simplificar, chamaremos <math>\frac{i \Delta t}{4 \Delta x^2}</math> de b. Nesse caso, temos:
<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - b\Psi_{j+1}^{n+1} - b\Psi_{j-1}^{n+1} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + b\Psi_{j+1}^{n} +b\Psi_{j-1}^{n}</math>
Definimos os Modos de Fourier:

<math>\Psi_{j}^{n} \approx A^{n}e^{iqjh}</math>

Assim, obtemos:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}e^{iqjh}- bA^{n+1}e^{iq(j-1)h} - bA^{n+1}e^{iq(j+1)h} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n}e^{iqjh} + bA^{n}e^{iq(j-1)h} +bA^{n}e^{iq(j+1)h}</math>

Dividindo tudo por <math>e^{iqjh}</math>:

<math> \left(1+2b+\frac{i\Delta t}{2} V_j\right)A^{n+1}- bA^{n+1}e^{-iqh} - bA^{n+1}e^{iqh} = \left(1-2b-\frac{i\Delta t}{2} V_j\right)A^{n} + bA^{n}e^{-iqh} +bA^{n}e^{iqh}</math>

Escrito de outra forma:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(e^{-iqh}-2+e^{iqh})\right] = A^{n}\left[1+\frac{i\Delta t}{2} V_j +b(e^{-iqh}-2+e^{iqh})\right] </math>

Utilizando relações trigonométricas:

<math> A^{n+1}\left[1+\frac{i\Delta t}{2} V_j -b(-sin^2(\frac{qh}2))\right] = A^{n}\left[1+\frac{i\Delta t}{2} V_j +b(-sin^2(\frac{qh}2))\right]</math>

Assim, obtemos o fator de amplificação:

<math>\frac{A^{n+1}}{A^n} = \frac{1+\frac{i\Delta t}{2} V_j -b(sin^2(\frac{qh}2))}{1+\frac{i\Delta t}{2} V_j +b(sin^2(\frac{qh}2))} \leq 1</math>

Portanto, o método é estavel

== Exemplos de Potenciais ==

=== Oscilador Harmônico ===

O oscilador harmônico quântico consiste de uma partícula em uma zona onde o potencial possui comportamento harmônico, ou seja, depende de <math> x^2 </math>. Em nosso exemplo, inserimos um pacote de formato gaussiano no interior de um potencial harmônico e observamos sua evolução temporal.

[[Arquivo:Harmonic_oscilator.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_Oh(x)
return 0.005(x-L/2)^2
end

b = dt*im/(4*dx^2)
A = [if i==j; 1 - 2b - V_Oh(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B = @. conj(A)
IB = inv(B)

len = length(x)
ψ = [(sqrt(1/(2*pi))*exp(-(val - L/3)^2/(2))*exp(-im*75*(val-L/3))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves = @. V_Oh(x)

t=0
@gif for t in 0:dt:65
ψ[2:end-1] = IB[2:end-1,2:end-1]*A[2:end-1,2:end-1]*ψ[2:end-1]

plot(x,abs2.(ψ),c="red")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves,c="blue")
plot!(ylim=[0,0.5],title="t=$t")
end every 3

</source> 

=== Barreira (Tunelamento) ===

Na mecânica quântica partículas podem existir, mesmo que com baixa probabilidade, em regiões que classicamente seriam proibidas. Este é o caso da barreira, onde uma partícula atravessa uma região onde sua energia total é negativa. Chamamos este fenômeno de tunelamento. A seguir, há um exemplo onde inserimos um pacote de formato gaussiano se deslocando em direção a uma barreira de potencial, de forma que parte da onda é refletida de volta e outra tunela pela barreira.

[[Arquivo:Barrier3.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_b(x)
if x==L/2
return 1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_b = [if i==j; 1 - 2b - V_b(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_b = @. conj(A_b)
IB_b = inv(B_b)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/3)^2/(8))*exp(-im*75*(val-L/3))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_b = @. V_b(x)

t=0
@gif for t in 0:dt:35
ψ[2:end-1] = IB_b[2:end-1,2:end-1]*A_b[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_b,c="blue",label="V(x)")
plot!(ylim=[0,0.2],title="t=$t")
end

</source> 

=== Poço Finito ===

O contrário da barreira é o chamado de poço de potencial. Ao invés da partícula atravessá-lo por completo, uma parte do pacote de onda é refletido, uma parte é transmitido e outra ainda fica preso no meio do poço. Veja o exemplo.

[[Arquivo:Pit.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_p(x)
if x<L/3
return 0.1
elseif x>2L/3
return 0.1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_p = [if i==j; 1 - 2b - V_p(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_p = @. conj(A_p)
IB_p = inv(B_p)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/4)^2/(8))*exp(-im*75*(val-L/4))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_p = @. V_p(x)

t=0
@gif for t in 0:dt:50
ψ[2:end-1] = IB_p[2:end-1,2:end-1]*A_p[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_p,c="blue",label="V(x)")
plot!(ylim=[0,0.15],title="t=$t")
end

</source>

Equação de Schrödinger Unidimensional

2024-04-23T16:54:27Z

Gustavobpasquali: /* O Método De Crank-Nicolson */

Conhecer o estado de uma partícula, na Mecânica Clássica, é saber, em um determinado momento, sua velocidade e sua posição. A partir destes valores, é possível determinar os seus futuros estados utilizando as relações fornecidas pelas Leis de Newton. Na Mecânica Quântica, porém, o estado físico não é mais caracterizado por uma quantidade discreta de valores numéricos, e sim, por uma função. A chamada Função de Onda, <math> \Psi </math>, nos informa a evolução temporal do estado quântico de uma partícula, seguindo a Equação de Schrödinger. Em uma dimensão, a Equação de Schrödinger é dada por:

<math> -\frac{\hbar}{2m}\frac{\partial^2\Psi}{\partial x^2} + V(x)\Psi = i\hbar\frac{\partial\Psi}{\partial t} </math>

== O Método De Crank-Nicolson==

A fim de resolvê-la numericamente para diferentes tipos de potenciais, usou-se o Método de Crank-Nicolson de resolução de EDPs. Este método baseia-se em combinar os métodos explícito e implícito (com igual contribuição) para aumentar sua estabilidade.

O método explícito FTCS (Forward Time Central Space) consiste em tomar derivadas temporais "para frente" e manter a derivada espacial centrada.

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t} </math>

<math> \frac{\partial\Psi(x,t)}{\partial x} = \frac{\Psi(x+\Delta x,t) - \Psi(x,t)}{\Delta x} </math>; <math> \frac{\partial^2\Psi(x,t)}{\partial^2 x} = \frac{(\Psi(x + \Delta x,t) - \Psi(x+\Delta x -\Delta x,t)) - ( \Psi(x,t) - \Psi(x - \Delta x,t) )}{\Delta x \Delta x} = \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} </math>

Tomando <math> \hbar=m=1 </math> (unidades atômicas), ao inserirmos as expressões acima na Equação de Schrödinger, temos que:

<math> \frac{i}{2} \frac{\Psi(x+\Delta x,t) + \Psi(x-\Delta x,t) - 2 \Psi(x,t)}{\Delta x ^2} -i V(x)\Psi(x,t) = \frac{\Psi(x,t+\Delta t) - \Psi(x,t)}{\Delta t}</math>

Simplificando a notação para <math> \Psi(x,t)=\Psi_j^n </math>, onde <math> j </math> representa a posição e <math> n </math> o tempo, e reorganizando os termos, ficamos com:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

O método implícito é semelhante, porém a derivada temporal é "para trás".

<math> \frac{\partial\Psi(x,t)}{\partial t} = \frac{\Psi(x,t) - \Psi(x,t-\Delta t)}{\Delta t} </math>

O que nos leva a ter:

<math> \Psi_j^{n} = \Psi_{j}^{n-1} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n}) - i\Delta t V_j \Psi_{j}^{n}</math>

que, com uma substituição de variável (que é muda), nos dá, enfim:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{2 \Delta x^2} (\Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - i\Delta t V_j \Psi_{j}^{n+1}</math>

Como dito, o Método de Crank-Nicolson é uma combinação dos métodos implícito e explícito, cada um tendo o mesmo peso de 1/2. Ao somarmos as duas equações (cada uma multiplicada por 0.5), ficaremos com a seguinte expressão:

<math> \Psi_j^{n+1} = \Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2} (\Psi_{j+1}^{n} + \Psi_{j-1}^{n} - 2 \Psi_{j}^{n} + \Psi_{j+1}^{n+1} + \Psi_{j-1}^{n+1} - 2 \Psi_{j}^{n+1}) - \frac{i\Delta t}{2} V_j (\Psi_{j}^{n} + \Psi_{j}^{n+1})</math>

que, ao reorganizarmos para isolar os <math> \Psi_j^{n+1}</math>, resulta em:

<math> \left(1+\frac{i \Delta t}{2 \Delta x^2}+\frac{i\Delta t}{2} V_j\right)\Psi_j^{n+1} - \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n+1} - \frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n+1} = \left(1-\frac{i \Delta t}{2 \Delta x^2}-\frac{i\Delta t}{2} V_j\right)\Psi_{j}^{n} + \frac{i \Delta t}{4 \Delta x^2}\Psi_{j+1}^{n} +\frac{i \Delta t}{4 \Delta x^2} \Psi_{j-1}^{n}</math>

Para resolver a equação em simultaneamente para ambos os tempos e em todos os pontos, escrevemos a equação em forma matricial:
:<math>
\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n+1} \\ \Psi_2^{n+1} \\ \vdots \\ \Psi_{j}^{n+1}\end{pmatrix}
=
\begin{pmatrix}
1 - 2b - \frac{i\Delta t V_j}{2} & +b & 0 & \cdots & 0 & 0 & 0\\ +b & 1 - 2b - \frac{i\Delta t V_j}{2} & +b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & +b & 1 - 2b - \frac{i\Delta t V_j}{2}
\end{pmatrix}
\begin{pmatrix} \Psi_1^{n} \\ \Psi_2^{n} \\ \vdots \\ \Psi_{j}^{n} \end{pmatrix} ;
</math>
com
<math> b=\frac{i \Delta t}{4 \Delta x^2}</math>.

Podemos notar que a diferença entre as matrizes que multiplicam <math> \Psi_j^{n+1}</math> e <math> \Psi_j^{n}</math> é somente o sinal das componente imaginárias. Assim, a equação final pode ser escrita como:

<math> A\Psi_j^{n+1} = A^*\Psi_j^{n}</math>, onde A é a matriz <math>\begin{pmatrix}
1 + 2b + \frac{i\Delta t V_j}{2} & -b & 0 & \cdots & 0 & 0 & 0\\ -b & 1 + 2b + \frac{i\Delta t V_j}{2} & -b &0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots& \vdots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & -b & 1 + 2b + \frac{i\Delta t V_j}{2}
\end{pmatrix} </math>

==Estabilidade do Método de Crank-Nicolson==

== Exemplos de Potenciais ==

=== Oscilador Harmônico ===

O oscilador harmônico quântico consiste de uma partícula em uma zona onde o potencial possui comportamento harmônico, ou seja, depende de <math> x^2 </math>. Em nosso exemplo, inserimos um pacote de formato gaussiano no interior de um potencial harmônico e observamos sua evolução temporal.

[[Arquivo:Harmonic_oscilator.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_Oh(x)
return 0.005(x-L/2)^2
end

b = dt*im/(4*dx^2)
A = [if i==j; 1 - 2b - V_Oh(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B = @. conj(A)
IB = inv(B)

len = length(x)
ψ = [(sqrt(1/(2*pi))*exp(-(val - L/3)^2/(2))*exp(-im*75*(val-L/3))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves = @. V_Oh(x)

t=0
@gif for t in 0:dt:65
ψ[2:end-1] = IB[2:end-1,2:end-1]*A[2:end-1,2:end-1]*ψ[2:end-1]

plot(x,abs2.(ψ),c="red")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves,c="blue")
plot!(ylim=[0,0.5],title="t=$t")
end every 3

</source> 

=== Barreira (Tunelamento) ===

Na mecânica quântica partículas podem existir, mesmo que com baixa probabilidade, em regiões que classicamente seriam proibidas. Este é o caso da barreira, onde uma partícula atravessa uma região onde sua energia total é negativa. Chamamos este fenômeno de tunelamento. A seguir, há um exemplo onde inserimos um pacote de formato gaussiano se deslocando em direção a uma barreira de potencial, de forma que parte da onda é refletida de volta e outra tunela pela barreira.

[[Arquivo:Barrier3.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_b(x)
if x==L/2
return 1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_b = [if i==j; 1 - 2b - V_b(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_b = @. conj(A_b)
IB_b = inv(B_b)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/3)^2/(8))*exp(-im*75*(val-L/3))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_b = @. V_b(x)

t=0
@gif for t in 0:dt:35
ψ[2:end-1] = IB_b[2:end-1,2:end-1]*A_b[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_b,c="blue",label="V(x)")
plot!(ylim=[0,0.2],title="t=$t")
end

</source> 

=== Poço Finito ===

O contrário da barreira é o chamado de poço de potencial. Ao invés da partícula atravessá-lo por completo, uma parte do pacote de onda é refletido, uma parte é transmitido e outra ainda fica preso no meio do poço. Veja o exemplo.

[[Arquivo:Pit.gif]]

Abaixo encontra-se o código usado em Júlia.

 

<source lang="julia">

using Plots

dt = 0.1
dx = 0.25
L = 50
x = collect(0:dx:L)
size(x)

function V_p(x)
if x<L/3
return 0.1
elseif x>2L/3
return 0.1
else
return 0
end
end

b = dt*im/(4*dx^2)
A_p = [if i==j; 1 - 2b - V_p(i)*dt*im/2 elseif i==dx+j || i==j-dx; b else 0 end for i=0:dx:L, j=0:dx:L]
B_p = @. conj(A_p)
IB_p = inv(B_p)

len = length(x)
ψ = [(sqrt(1/(8*pi))*exp(-(val - L/4)^2/(8))*exp(-im*75*(val-L/4))) for val in x]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end

ves_p = @. V_p(x)

t=0
@gif for t in 0:dt:50
ψ[2:end-1] = IB_p[2:end-1,2:end-1]*A_p[2:end-1,2:end-1]*ψ[2:end-1]
for (i,num) in enumerate(ψ)
if abs(num)<1e-10;
ψ[i]=0
end
end
plot(x,abs2.(ψ),c="red",label="|ψ|^2")
#plot!([x for x=0:dx:L],imag(ψ),c="blue")
plot!([x for x=0:dx:L],ves_p,c="blue",label="V(x)")
plot!(ylim=[0,0.15],title="t=$t")
end

</source>