Trabalhos 2022/2: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Desfeita a edição 9611 de Dzorrer (Discussão))
 
(10 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 1: Linha 1:
=== [[ Equações de Laplace e Poisson - Eletrostática ]] ===
=== [[ Equações de Laplace e Poisson - Eletrostática ]] ===
=== [[Shooting method e Método de Crank-Nicolson aplicados à Equação de Schrödinger]] ===
=== [[Shooting method e Método de Crank-Nicolson]] ===
=== [[ Lennard-Jones e propriedades termodinâmicas ]] ===
O objetivo deste trabalho é aplicar o Shooting method (método do chute) para encontrar as primeiras funções de onda espaciais da Equação de Schrödinger para o caso do poço de potencial infinito. Após, será realizada a evolução temporal através do Método de Crank-Nicolson.
=== [[ Modelo de Blume-Capel bidimensional ]] ===
 
==Equação de Schrödinger==
 
A equação de Schrödinger unidimensional pode ser escrita da seguinte maneira:
<center><math> i\hbar\frac{\partial\Psi}{\partial t}=-\frac{\hbar^2}{2m}\frac{\partial^2 \Psi}{\partial x^2}+V\Psi.</math></center>
Para resolvê-la é necessário efetuar uma separação de variáveis:
<center><math>\Psi(x,t)=\psi(x)\varphi(t).</math></center> 
 
Aplicando na primeira equação e separando os termos espaciais dos termos temporais, chega-se a uma equação com o seguinte formato:
 
<center><math>i\hbar\frac{1}{\varphi}\frac{d\varphi}{dt}=-\frac{\hbar^2}{2m}\frac{1}{\psi}\frac{d^2\psi}{dx^2}+V. </math></center>
 
Pelo fato da parte da esquerda ser dependente de ''t'' e a parte da direita ser dependente de ''x'' e de ambas estarem relacionadas por uma igualdade, é necessário que ambos os lados sejam constantes: em outras palavras, não é possível modificar um lado sem necessariamente alterar o outro. Através de um raciocínio perspicaz, a constante em questão será denominada ''E''.
 
===Parte temporal===
A parte que diz respeito à evolução temporal:
<center><math> \frac{d\varphi}{dt}=-\frac{iE}{\hbar}\varphi.</math></center>
A solução geral possui o seguinte formato
<center><math> \varphi(t)=Ce^{-\frac{-iEt}{\hbar}} ,</math></center>
cuja constante C pode, neste caso, ser absorvida, de modo que
<center><math> \varphi(t)=e^{-\frac{-iEt}{\hbar}} . </math></center>
 
===Parte espacial===
Quanto à parte espacial, utilizando o mesmo raciocínio empregado anteriormente, a equação pode ser escrita como
<center><math>-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2}+V\psi=E\psi </math></center>
 
Para este caso, no entanto, não há uma única solução, pois esta depende do potencial ''V'' escolhido. Para o presente trabalho optou-se por trabalhar com o caso do poço infinito de potencial pelo fato das soluções analíticas já serem conhecidas, de modo a tornar possível avaliar os resultados numéricos obtidos à luz da solução analítica.
 
==Poço de potencial infinito==
Esquematicamente, tem-se:
[[Arquivo:Poço.png|200px|thumb|center|Poço de potencial infinito]]
O potencial pode ser descrito como:
<center><math>
  V(x) =
      \begin{cases}
            0, & \mbox{se } 0\leq x\leq L, \\
            \infty, & \mbox{de outra forma.}
      \end{cases}
</math></center>
 
Dentro do poço, onde $V=0$, o problema pode ser modelado da seguinte maneira
<center><math>-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2}=E\psi , </math></center>
ou
<center><math>\frac{d^2\psi}{dx^2}=-k^2\psi,</math></center> onde <center><math>k\equiv \frac{\sqrt{2mE}}{\hbar}.</math></center>
 
A solução é dada por
<center><math>\psi(x)=Asen(kx)+Bcos(kx).</math></center>
 
Aplicando as condições de contorno <math>\psi(0)=\psi(L)=0 </math> e efetuando a normalização da função de onda, obtém-se a solução geral
<center><math>\psi_n(x)=\sqrt{\frac{2}{L}}sen\left(\frac{n\pi}{L}x\right), </math></center>
 
cujas energias discretizadas são
 
<center><math>E_n=\frac{\hbar^2k_n^2}{2m}=\frac{n^2\pi^2\hbar^2}{2mL^2}.</math></center>
 
Utilizando a equação acima, pode-se calcular os valores da energia de cada estado estacionário. Para o caso de um elétron, as energias referentes aos três estados estacionários são <math>E_1=0,376</math> eV, <math>E_2=1,504</math> eV e <math>E_3=3,384</math> eV.
 
Na próxima seção será feita uma estimativa dos valores acima expostos através do "Shooting method".
 
==Shooting Method==
Muitos métodos numéricos (e.g. Runge-Kutta, Forward Euler) requerem os valores da função e de sua derivada no ponto inicial. Acontece que podem haver problemas em que estes valores não estarão disponíveis, principalmente o valor da derivada em questão. Uma alternativa seria conjecturar o valor da condição inicial e integrar, através de um método apropriado, em direção à outra condição de contorno: um "chute" apropriado faria com que a integração evoluísse e retornasse um valor muito próximo, a depender da acurácia necessária, ao da condição de contorno. A ideia seria executar os seguintes passos:
 
# Supor um valor para a condição de contorno desconhecida (e.g. <math>y(0)</math> ou <math>y'(0)</math>);
# Integrar o problema através de um método conhecido até a próxima condição de contorno (e.g., <math>y(L)</math>);
# Se o chute inicial não fez com que o sistema evoluísse até <math>y(L)</math>, então deve-se supor outro valor para a condição inicial e repetir o procedimento.
 
O método descrito acima de forma simplificada recebe o nome, em inglês, de ''Shooting method'', o que em português seria algo como "Método do tiro" ou "Método do chute". Na próxima seção esse método será aplicado para o caso do poço infinito de potencial.
 
===Poço de potencial infinito===
Seja a equação <math>\frac{d^2\psi}{dx^2}=-k\psi E</math>, onde <math>k=\frac{2m}{\hbar^2}</math>.
 
Escrevendo com outra notação: <math>\ddot{\psi}=-k\psi E</math>.
 
Dividindo o problema em <math>\Delta x </math>'s pequenos, pode-se reescrever a equação acima da seguinte forma:
 
<center><math> \ddot{\psi}=\frac{\Delta \dot{\psi}}{\Delta x}=\frac{\dot{\psi_2}-\dot{\psi_1}}{\Delta x} \implies \dot{\psi_2}=\ddot{\psi} \Delta x + \dot{\psi_1}</math></center>.
 
Também:
<center><math>\dot{\psi}=\frac{\Delta \psi}{\Delta x}=\frac{\psi_2-\psi_1}{\Delta x} \implies \psi_2 = \dot{\psi}\Delta x + \psi_1</math></center>.
 
Além disso:
<center><math>\Delta x = x_2-x_1 \implies x_2 = x_1 +\Delta{x}</math></center>.
 
A integração, então, é realizada utilizando as relações 8, 9, 10 e 11, até que se atinja a borda do poço, isto é, <math>x=L</math>.
 
Com a discretização acima, foi possível implementar o algoritmo. Das condições de contorno do problema, sabe-se que <math>\psi(0)=0</math>, de modo que <math>\psi_1=0</math>. No entanto, o valor da derivada <math>\dot{\psi_1}</math> não é conhecido, de modo que supõe-se que seja uma constante, a saber, <math>\dot{\psi_1}=1</math>. Chutando que <math>E=0</math>, utilizando a massa do elétron e <math>L=1</math>, obtém-se a primeira solução estacionária:
[[Arquivo:n=1.png|300px|thumb|center|Solução estacionária (n=1)]]
Pode-se observar que o valor de energia obtido numericamente é cerca de 4% menor do que aquele obtido analiticamente.
 
Para o caso n=2:
[[Arquivo:n=2.png|300px|thumb|center|Solução estacionária (n=2)]]
Aqui, o valor obtido numericamente é aproximadamente 5% maior do que o valor obtido analiticamente.
 
Para o caso n=3:
[[Arquivo:n=3.png|300px|thumb|center|Solução estacionária (n=3)]]
Para este caso, o valor numérico é cerca de 1% menor do que o valor analítico.
 
==Método de Crank-Nicolson==
Seja a equação diferencial
<center><math>\frac{\partial f}{\partial t}=L_1{\textbf{r}}f(\textbf{r},t)</math></center>,
onde <math>L_{\textbf{r}}</math> é um operador diferencial linear em '''r'''.
 
Em forma discretizada no tempo, pode-se escrever
<center><math>
f^{n+1}(\textbf{r})-f^{n}(\textbf{r})= L_{\textbf{r}}f^{n}(\textbf{r})dt
</math></center>.
 
Por simetria, pode-se escrever a equação acima utilizando um f à direita:
<center><math>
f^{n+1}(\textbf{r})-f^{n}(\textbf{r})= L_{\textbf{r}}f^{n+1}(\textbf{r})dt .
</math></center>
 
A equação acima é dita "explícita" pois, para o cálculo de <math>f^{n+1}</math>, só é utilizado o valor já explicitamente calculado <math>f^{n}</math>. Já a equação anterior é chamada implícita pois <math>f^{n+1}</math> está presente explicitamente. Em termos numéricos, um método peca pelo excesso enquanto o outro o faz pela falta, de modo que um resultado mais satisfatório pode ser obtido ao tomar-se a média dos dois:
<center><math>
f^{n+1}(\textbf{r})-f^{n}(\textbf{r})=\frac{dt}{2}(L_{\textbf{r}}f^{n+1}(\textbf{r})+L_{\textbf{r}}f^n(\textbf{r})).
</math></center>
 
Após alguma álgebra:
<center><math> f^{n+1}(\textbf{r})=\left(1-\frac{dt}{2}L_{\textbf{r}}\right)^{-1}\left(1+\frac{dt}{2}L_{\textbf{r}}\right)f^{n}(\textbf{r}) </math></center>.
 
Chamando <math> M=I+\frac{dt}{2}L_{\textbf{r}} </math> e <math> E=I-\frac{dt}{2}L_{\textbf{r}} </math>, onde ''I'' indica a matriz identidade, pode-se reescrever a equação acima na seguinte maneira:
 
<center><math> f^{n+1}=E^{-1}Mf^{n} </math></center>.
 
Trata-se do método de Crank-Nicolson, mais estável e preciso do que os métodos implícito e explícito. Caso o problema apresentar condições de contorno, estas serão devidamente implementadas nos elementos das matrizes ''M'' e ''E''.
 
===Equação de Schrödinger===
Seja a equação de Schrödinger unidimensional
 
<center><math> i\hbar\frac{\partial\Psi}{\partial t}=-\frac{\hbar^2}{2m}\frac{\partial^2 \Psi}{\partial x^2}+V\Psi </math></center>.
 
Efetuando a discretização das variáveis através do Método de Crank-Nicolson, obtém-se:
 
<center><math>\frac{\partial \Psi}{\partial t}= \frac{\Psi_{j}^{n+1}-\Psi_{j}^{n}}{\Delta t} ;</math></center>
 
<center><math>\frac{\hbar^2}{2m}\frac{\partial^2 \Psi}{\partial x^2} = \frac{\hbar^2}{2m}\left[\frac{(\Psi_{j+1}^{n+1}-2\Psi_{j}^{n+1}+\Psi_{j-1}^{n+1})+ (\Psi_{j+1}^{n}-2\Psi_{j}^{n}+\Psi_{j-1}^{n})}{2\Delta x^2}\right] ;</math></center>
 
<center><math>V\Psi = \frac{1}{2}[V_{j}^{n+1}\Psi_{j}^{n+1}+V_{j}^{n}\Psi_{j}^{n}] .</math></center>
 
Substituindo as discretizações na eq. de Schrödinger:
 
<center><math>i\left(\frac{\Psi_{j}^{n+1}-\Psi_{j}^{n}}{\Delta t}\right)=-\frac{\hbar^2}{4m(\Delta x)^2} \left[(\Psi_{j+1}^{n+1}-2\Psi_{j}^{n+1}+\Psi_{j-1}^{n+1})+ (\Psi_{j+1}^{n}-2\Psi_{j}^{n}+\Psi_{j-1}^{n})\right]+\frac{1}{2}[V_{j}^{n+1}\Psi_{j}^{n+1}+V_{j}^{n}\Psi_{j}^{n}].</math></center>
 
Supondo <math>\hbar=m=1</math> e separando as partes explícita e implícita, obtém-se, após alguma álgebra:
 
<center><math>\Psi_{j}^{n+1}\left[1+\frac{i\Delta t}{2}\left(\frac{1}{\Delta x^2}+V_{j}^{n+1}\right)\right]+\Psi_{j-1}^{n+1}\left[-\frac{i\Delta t}{4\Delta x^2}\right]+\Psi_{j+1}^{n+1}\left[-\frac{i\Delta t}{4\Delta x^2}\right] = \Psi_{j}^{n}\left[1-\frac{i\Delta t}{2}\left(\frac{1}{\Delta x^2}+V_{j}^{n}\right)\right]+\Psi_{j-1}^{n}\left[-\frac{i\Delta t}{4\Delta x^2}\right]+\Psi_{j+1}^{n}\left[\frac{i\Delta t}{4\Delta x^2}\right].</math></center>
 
Definindo
<center><math>a\equiv\frac{-i\Delta t}{4(\Delta x)^2}</math></center> e <center><math>b_{j}\equiv\left(1+\frac{i\Delta t}{2}\right)\left(\frac{1}{\Delta x^2}+V_j\right),</math></center> obtém-se:
 
<center><math>\Psi_{j}^{n+1}b_{j}+\Psi_{j-1}^{n+1}a+\Psi_{j+1}^{n+1}a = \Psi_{j}^{n}b_{j}^{*}+\Psi_{j-1}^{n}a^{*}+\Psi_{j+1}^{n}a^{*} .</math></center>
 
A equação acima pode ser escrita em forma matricial, de modo que: 
<center><math>\hat{C}\Psi^{n+1}=\hat{D}\Psi^{n},</math></center>
onde
 
<center><math>
\hat{C} =
\begin{bmatrix}
b _{0}           &a     &0 & &...  &0 \\
a & b_{1}   &a  &0 &...     &0\\
0 & a     &b_{2} &a & 0     &0 \\
0 & 0   & \ddots & \ddots &\ddots \\
...           & ...   &... & a            &b_{j-1} &a \\
0           & 0   &... & 0            &a &b_{j}\\
\end{bmatrix}
</math></center>
e
<center><math>
\hat{D} =
\begin{bmatrix}
b_{0}^{*}           &a     &0 & &...  &0 \\
a^{*} & b_{1}^{*}   &a^{*} &0 &...     &0\\
0 & a^{*}     &b_{2}^{*} &a^{*} & 0     &0 \\
0 & 0   & \ddots & \ddots &\ddots \\
...           &...   &... & a^{*}            &b_{j-1}^{*} &a^{*}  \\
0           & 0   &... & 0            &a^{*}  &b_{j}^{*}\\
\end{bmatrix}
</math></center>
 
Para avaliar a evolução temporal do sistema, é necessário encontrar <math>\Psi^{n+1}</math>. Utilizando resultados anteriores, pode-se obter <math>\Psi^{n+1}</math> através da seguinte relação:
<center><math>
\Psi^{n+1}=\hat{C}^{-1}\hat{C^{*}}\Psi^{n}
</math></center>
 
===Poço de potencial infinito===
 
Para o presente caso a ideia é obter a evolução temporal do sistema, impondo condições de contorno iguais a zero, de modo que os operadores <math>\hat{C}</math> e <math>\hat{D}</math> ficam:
 
<center><math>
\hat{C} =
\begin{bmatrix}
1           &0     &0 & &...    &0 \\
a & b   &a  &0 &...     &0\\
0 & a     &b &a & 0     &0 \\
0 & 0   & \ddots & \ddots &\ddots \\
...           &...   &... & a            &b &a \\
0           & 0   &... & 0            &0 &1\\
\end{bmatrix}
</math></center>
e
 
<center><math>
\hat{D} =
\begin{bmatrix}
1           &0     &0 & &...    &0 \\
a^{*}  & b^{*}   &a^{*}  &0 &0...     &0\\
0 & a^{*}     &b^{*} &a^{*} & 0     &0 \\
0 & 0   & \ddots & \ddots &\ddots \\
...           &...   &... & a^{*}            &b^{*} &a^{*}  \\
0           & 0   &... & 0            &0 &1\\
\end{bmatrix}
</math></center>
 
A ideia é que o primeiro e o último termo do tanto do vetor <math>\Psi^{n}</math> quanto do vetor <math>\Psi^{n+1}</math> seja constante, o que satisfaz as condições de contorno do presente caso. Também é interessante notar que os índices <math>b</math> são todos constantes, visto que no presente caso o potencial dentro do poço é nulo.
 
Implementando o algoritmo descrito acima, obteve-se:
[[Arquivo:Animação_n_1.gif|300px|thumb|center|Evolução temporal (n=1)]]
Evolução temporal para o caso n=1. Nesta animação e nas subsequentes, foram sobrepostas as partes real e imaginária da equação de Schrödinger: a linha azul diz respeito à parte real enquanto a amarela, à imaginária.
 
[[Arquivo:Animação_n_2.gif|300px|thumb|center|Evolução temporal (n=2)]]
Na figura acima, tem-se a evolução do caso n=2.
[[Arquivo:Animação_n_3.gif|300px|thumb|center|Evolução temporal (n=3)]]
Por último, o caso n=3.

Edição atual tal como às 15h56min de 18 de abril de 2023

Equações de Laplace e Poisson - Eletrostática

Shooting method e Método de Crank-Nicolson

Lennard-Jones e propriedades termodinâmicas

Modelo de Blume-Capel bidimensional

Disponível em “http://fiscomp.if.ufrgs.br/index.php?title=Trabalhos_2022/2&oldid=9612”