Shooting method e Método de Crank-Nicolson

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.

Equação de Schrödinger

A equação de Schrödinger unidimensional pode ser escrita da seguinte maneira:

$i\hslash \frac{\partial \mathrm{\Psi }}{\partial t}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2\mathrm{\Psi }}{\partial x^2}+V\mathrm{\Psi }.$

Para resolvê-la é necessário efetuar uma separação de variáveis:

$\mathrm{\Psi }(x,t)=\psi (x)\phi (t).$

Aplicando na primeira equação e separando os termos espaciais dos termos temporais, chega-se a uma equação com o seguinte formato:

$i\hslash \frac{1}{\phi }\frac{d\phi }{dt}=-\frac{{\hslash }^2}{2m}\frac{1}{\psi }\frac{d^2\psi }{d{x}^2}+V.$

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:

$\frac{d\phi }{dt}=-\frac{iE}{\hslash }\phi .$

A solução geral possui o seguinte formato

$\phi (t)=C{e}^{-\frac{-iEt}{\hslash }},$

cuja constante C pode, neste caso, ser absorvida, de modo que

$\phi (t)=e^{-\frac{-iEt}{\hslash }}.$

Parte espacial

Quanto à parte espacial, utilizando o mesmo raciocínio empregado anteriormente, a equação pode ser escrita como

$-\frac{{\hslash }^2}{2m}\frac{d^2\psi }{d{x}^2}+V\psi =E\psi$

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:

Poço de potencial infinito

O potencial pode ser descrito como:

$V(x)=\{\begin{array}{ll}0,& \text{se }0\le x\le L,\\ \mathrm{\infty },& \text{de outra forma.}\end{array}$

Dentro do poço, onde $V=0$, o problema pode ser modelado da seguinte maneira

$-\frac{{\hslash }^2}{2m}\frac{d^2\psi }{d{x}^2}=E\psi ,$

ou

$\frac{d^2\psi }{d{x}^2}=-k^2\psi ,$

onde

$k\equiv \frac{\sqrt{2mE}}{\hslash }.$

A solução é dada por

$\psi (x)=Asen(kx)+Bcos(kx).$

Aplicando as condições de contorno $\psi (0)=\psi (L)=0$ e efetuando a normalização da função de onda, obtém-se a solução geral

${\psi }_n(x)=\sqrt{\frac{2}{L}}sen\left(\frac{n\pi }{L}x\right),$

cujas energias discretizadas são

$E_n=\frac{{\hslash }^2{k}_n^2}{2m}=\frac{n^2{\pi }^2{\hslash }^2}{2m{L}^2}.$

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 $E_1=0,376$ eV, $E_2=1,504$ eV e $E_3=3,384$ 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:

1. Supor um valor para a condição de contorno desconhecida (e.g. $y(0)$ ou $y^{\prime }(0)$);
2. Integrar o problema através de um método conhecido até a próxima condição de contorno (e.g., $y(L)$);
3. Se o chute inicial não fez com que o sistema evoluísse até $y(L)$, 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 $\frac{d^2\psi }{d{x}^2}=-k\psi E$, onde $k=\frac{2m}{{\hslash }^2}$.

Escrevendo com outra notação: $\ddot{\psi }=-k\psi E$.

Dividindo o problema em $\mathrm{\Delta }x$'s pequenos, pode-se reescrever a equação acima da seguinte forma:

$\ddot{\psi }=\frac{\mathrm{\Delta }\dot{\psi }}{\mathrm{\Delta }x}=\frac{\dot{{\psi }_2}-\dot{{\psi }_1}}{\mathrm{\Delta }x}⟹\dot{{\psi }_2}=\ddot{\psi }\mathrm{\Delta }x+\dot{{\psi }_1}$

.

Também:

$\dot{\psi }=\frac{\mathrm{\Delta }\psi }{\mathrm{\Delta }x}=\frac{{\psi }_2-{\psi }_1}{\mathrm{\Delta }x}⟹{\psi }_2=\dot{\psi }\mathrm{\Delta }x+{\psi }_1$

.

Além disso:

$\mathrm{\Delta }x=x_2-x_1⟹x_2=x_1+\mathrm{\Delta }x$

.

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 é, $x=L$.

Com a discretização acima, foi possível implementar o algoritmo. Das condições de contorno do problema, sabe-se que $\psi (0)=0$, de modo que ${\psi }_1=0$. No entanto, o valor da derivada $\dot{{\psi }_1}$ não é conhecido, de modo que supõe-se que seja uma constante, a saber, $\dot{{\psi }_1}=1$. Chutando que $E=0$, utilizando a massa do elétron e $L=1$, obtém-se a primeira solução estacionária:

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:

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:

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

$\frac{\partial f}{\partial t}=L_1\text{r}f(\text{r},t)$

,

onde $L_{\text{r}}$ é um operador diferencial linear em r.

Em forma discretizada no tempo, pode-se escrever

$f^{n+1}(\text{r})-f^n(\text{r})=L_{\text{r}}{f}^n(\text{r})dt$

.

Por simetria, pode-se escrever a equação acima utilizando um f à direita:

$f^{n+1}(\text{r})-f^n(\text{r})=L_{\text{r}}{f}^{n+1}(\text{r})dt.$

A equação acima é dita "explícita" pois, para o cálculo de $f^{n+1}$, só é utilizado o valor já explicitamente calculado $f^n$. Já a equação anterior é chamada implícita pois $f^{n+1}$ 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:

$f^{n+1}(\text{r})-f^n(\text{r})=\frac{dt}{2}(L_{\text{r}}{f}^{n+1}(\text{r})+L_{\text{r}}{f}^n(\text{r})).$

Após alguma álgebra:

$f^{n+1}(\text{r})={(1-\frac{dt}{2}{L}_{\text{r}})}^{-1}(1+\frac{dt}{2}{L}_{\text{r}})f^n(\text{r})$

.

Chamando $M=I+\frac{dt}{2}{L}_{\text{r}}$ e $E=I-\frac{dt}{2}{L}_{\text{r}}$, onde I indica a matriz identidade, pode-se reescrever a equação acima na seguinte maneira:

$f^{n+1}=E^{-1}M{f}^n$

.

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

$i\hslash \frac{\partial \mathrm{\Psi }}{\partial t}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2\mathrm{\Psi }}{\partial x^2}+V\mathrm{\Psi }$

.

Efetuando a discretização das variáveis através do Método de Crank-Nicolson, obtém-se:

$\frac{\partial \mathrm{\Psi }}{\partial t}=\frac{{\mathrm{\Psi }}_j^{n+1}-{\mathrm{\Psi }}_j^n}{\mathrm{\Delta }t};$ $\frac{{\hslash }^2}{2m}\frac{{\partial }^2\mathrm{\Psi }}{\partial x^2}=\frac{{\hslash }^2}{2m}\left[\frac{({\mathrm{\Psi }}_{j+1}^{n+1}-2{\mathrm{\Psi }}_j^{n+1}+{\mathrm{\Psi }}_{j-1}^{n+1})+({\mathrm{\Psi }}_{j+1}^n-2{\mathrm{\Psi }}_j^n+{\mathrm{\Psi }}_{j-1}^n)}{2\mathrm{\Delta }{x}^2}\right];$ $V\mathrm{\Psi }=\frac{1}{2}[V_j^{n+1}{\mathrm{\Psi }}_j^{n+1}+V_j^n{\mathrm{\Psi }}_j^n].$

Substituindo as discretizações na eq. de Schrödinger:

$i\left(\frac{{\mathrm{\Psi }}_j^{n+1}-{\mathrm{\Psi }}_j^n}{\mathrm{\Delta }t}\right)=-\frac{{\hslash }^2}{4m(\mathrm{\Delta }x{)}^2}[({\mathrm{\Psi }}_{j+1}^{n+1}-2{\mathrm{\Psi }}_j^{n+1}+{\mathrm{\Psi }}_{j-1}^{n+1})+({\mathrm{\Psi }}_{j+1}^n-2{\mathrm{\Psi }}_j^n+{\mathrm{\Psi }}_{j-1}^n)]+\frac{1}{2}[V_j^{n+1}{\mathrm{\Psi }}_j^{n+1}+V_j^n{\mathrm{\Psi }}_j^n].$