8.4: Ejemplo - Problema de partículas en una caja

Para demostrar el uso de matrices dispersas para resolver ecuaciones de diferencia finita, consideramos el problema de partículas en caja 1D. Esto consiste en la ecuación de onda de Schrödinger independiente del tiempo 1D,

\[-\frac{1}{2} \, \frac{d^2\psi}{dx^2} = E\psi(x), \quad 0 \le x \le L,\]

junto con las condiciones de límite de Dirichlet

\[\psi(0) = \psi(L) = 0.\]

La solución analítica es bien conocida por nosotros; hasta un factor de normalización, los autoestados y los valores propios de energía son

\[\psi_m(x) = \sin\Big(m\pi x / L\Big),\quad E_m = \frac{1}{2}\, \left(\frac{m\pi}{L}\right)^2, \quad m = 1, 2, 3, \dots\]

Escribiremos un programa que busque una solución numérica. Usando la regla de tres puntos para discretizar la segunda derivada, las ecuaciones matriciales de diferencias finitas se convierten en

\[-\frac{1}{2h^2}\begin{bmatrix} -2 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & 1 \\ & & 1 & -2\end{bmatrix} \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix} = E \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix},\]

donde

\[h = \frac{L}{N + 1}, \quad \psi_n = \psi\Big(x = h(n+1)\Big)\]

El siguiente programa construye la ecuación de matriz de diferencia finita y muestra las tres primeras soluciones:

from scipy import *
import scipy.sparse as sp
import scipy.sparse.linalg as spl
import matplotlib.pyplot as plt

## Solve the 1D particle-in-a-box problem for box length L,
## using N discretization points.  The parameter nev is the
## number of eigenvalues/eigenvectors to find.  Return three
## arrays E, psi, and x.  E stores the energy eigenvalues;
## psi stores the (non-normalized) eigenstates; and x stores
## the discretization points.
def particle_in_a_box(L, N, nev=3):
    dx = L/(N+1.0)
    x  = linspace(dx, L-dx, N)
    I  = ones(N)
    ## Set up the finite-difference matrix.
    H  = sp.dia_matrix(([I, -2*I, I], [-1,0,1]), shape=(N,N))
    H *= -0.5/(dx*dx)
    ## Find the lowest eigenvalues and eigenvectors.
    E, psi = spl.eigsh(H, k=nev, sigma=-1.0)
    return E, psi, x

def particle_in_a_box_demo():
    E, psi, x = particle_in_a_box(1.0, 1000)

    ## Print the energy eigenvalues.
    print(E)

    ## Plot |psi(x)|^2 vs x for each solution found.
    fig = plt.figure()
    axs = plt.subplot(1,1,1)
    for n in range(len(E)):
        fig_label = "State #" + str(n)
        plt.plot(x, abs(psi[:,n])**2, label=fig_label, linewidth=2)
    plt.xlabel('x')
    plt.ylabel('|psi(x)|^2')

    ## Shrink the axis by 20%, so the legend can fit.
    box = axs.get_position()
    axs.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ## Print the legend.
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    plt.show()

particle_in_a_box_demo()

La matriz hamiltoniana, que es tridiagonal, se construye en el formato de matriz dispersa DIA. Los valores propios y los vectores propios se encuentran con eigsh, los cuales pueden ser utilizados debido a que se sabe que la matriz hamiltoniana es hermitiana. Observe que llamamos eigsh usando el parámetro sigma, diciéndole al solucionador propio que encuentre los valores propios más cercanos en valor a\(-1.0\):

    E, psi = spl.eigsh(H, k=nev, sigma=-1.0)

Esto encontrará los valores propios de energía más bajos porque, en este caso, todos los valores propios de energía son positivos. (Usamos\(-1.0\) en lugar de 0.0, porque el algoritmo no funciona bien cuando sigma es exactamente cero.) Si hay un potencial negativo presente, tendríamos que encontrar una estimación diferente para el límite inferior de los valores propios de energía, y pasarlo a sigma.

Alternativamente, podríamos haber llamado eigsh con una entrada que = 'sa'. Esto le diría al solucionador propio que encuentre el valor propio con el valor más pequeño. Evitamos hacer esto porque el algoritmo ARPACK eigensolver es relativamente ineficiente para encontrar pequeños autovalores en qué modo (y a veces incluso puede fallar en converger, si k es demasiado pequeño). Normalmente, si eres capaz de deducir un límite inferior para los valores propios deseados, es preferible usar sigma.

Al ejecutar el programa se imprime los valores propios de energía más bajos:

[  4.93479815  19.73914399  44.41289171]

Esto concuerda bien con los resultados analíticos\(E_1 = \pi^2/2 = 4.934802\),\(E_2 = 2\pi^2 = 19.739208\), y\(E_3=9\pi^2/2 = 44.413219\). También produce la trama que se muestra a continuación, la cual es igualmente la esperada.

Figura\(\PageIndex{1}\): Gráfica de\(|\psi(x)|^2\) versus\(x\) para el problema de partículas en caja.