Solución numérica de la ecuación de Schrödinger - valores propios

Me alegro de que finalmente pueda responder a esto.

El método de Numerov como se describe en Wikipedia no es cómo desea proceder aquí. Para darle una idea de cómo proceder, comencemos con una versión simplificada del método. Lo que voy a hacer es discretizar ingenuamente el operador diferencial, así:

La última ecuación es solo una definición: estoy tratando el espacio como si fuera una red con puntos$d$aparte y estoy llamando$\psi_n$el valor de la función de onda en el$n$th point.. Ahora la ecuación de Schrödinger se lee en esta notación:

¡ Pero esta es una ecuación matricial ! Permítanme ser más explícito:

Por supuesto, tuve que elegir algún número entero.$N$eso solo corresponde a colocar el sistema en una caja que es "lo suficientemente grande" para tener un sistema finito.

Ahora está claro que lo que tiene es un problema de valor propio de matriz de la forma

y puede proceder a diagonalizarlo de la forma que elija. Tenga en cuenta que llamamos$A$una matriz tridiagonal, por razones obvias. Es posible que desee ocuparse primero de las condiciones de contorno; lo hace configurando$\psi_1 = \psi_n = 0$antes de construir la matriz, lo que corresponde a establecer la primera y la última columna en cero. Esto le dará algunas funciones propias falsas con valor propio cero que simplemente puede descartar. Si tiene diferentes condiciones de contorno, no tiene suerte; no conozco una forma que lo haga funcionar en general.

Ahora solo tienes que rehacer lo anterior con el método completo de Numerov, que será un poco más complicado, y ya está todo listo.

Esta parece ser la forma que está buscando pero, por supuesto, no es la única forma de hacerlo. Griffiths describe uno llamado "movimiento del perro" que consiste en lo siguiente: adivinas un valor inicial para una energía y calculas la función de onda como quieras (RK, por ejemplo). Lo más probable es que no sea un valor propio, por lo que la función de onda explotará en el infinito. podría ir a$+\infty$para grande$x$o podría ir a$-\infty$. Ahora va variando lentamente la energía hasta que el comportamiento en el infinito "cambie", es decir, hasta que la cola "se mueva". Eso le permitirá restringir el valor de un solo valor propio de energía y le dará la forma de la función de onda a un tamaño máximo que aumenta a medida que se elige.$E$se aproxima al valor propio exacto.

Puede buscar valores propios utilizando el método de bisección.

Priliminares: para obtener los valores propios del método de Numerov, necesitará conocer la función de onda en los límites. En general, esto significaría que debe establecer el potencial en infinito en los límites, por lo tanto, poner la función de onda en cero en esos puntos. Para su potencial, modifíquelo de la siguiente manera:

V = infinito si -50fm>r>50fm V = 0 si |r|>8.5fm V = potencial de Woods en caso contrario

El trato real: ahora escriba un programa para calcular la función de onda en el potencial anterior utilizando cualquier energía arbitraria en el dominio -50fm

Aquí hay un código de Python que encuentra un valor propio en el intervalo dado (E1, E2) si existe.

import numpy as np

psiarray = []
r = np.linspace(rmin, rmax, N)
delta = r[1]-r[0]

#Assume psi(rmin) = psibcmin and psi(rmin+delta) = 1 and then
#calculate phi0 and phi1
psiarray.append(psibcmin)
psiarray.append(1)
phi0 = psitophi(psiarray[0], l, E, r[0], delta) #r[0]=rmin
phi1 = psitophi(psiarray[1], l, E, r[1], delta)

#Now populate the psiarray for each value of r
for i in range(2, len(r)):
    phi2 = 2*phi1 - phi0 + (delta**2)*f(l, E, r[i-1])*psiarray[i-1]
    psi = phitopsi(phi2, l, E, r[i], delta)
    psiarray.append(psi)
    phi0 = phi1
    phi1 = phi2

return r, psiarray

for i in range(len(psi)):
    psi[i] = psi[i]/area

return psi

    if psi1*psi2 < 0:
        emid = e1 + (e2 - e1) * 0.5
        r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, emid)
        psimid = psi.pop()

        if psimid*psi1 < 0:
            e2 = emid

        elif psimid*psi2 < 0:
            e1 = emid
    elif psi1*psi2 > 0:
        print "There are either no eigenvalues or too many of them in"+\
        " the given energy interval. For e1={} psi={} and e2={} psi=".format(e1, psi1, e2)+\
        "{}".format(psi2)
        return e1,e2
return e1,e2

e1 = float(sys.argv[1])
e2 = float(sys.argv[2])
l = int(sys.argv[3])
if l == 0:
    rmin = -1
else:
    rmin = 0
tol =1e-5
rmax = 1
psibcmin = 0
psibcmax = 0
N = 100

e1, e2 = locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax,\
        l, e1, e2, tol)

fig = plt.figure()
ax = fig.add_subplot(111)
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
ax.plot(r,psi, 'o')
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l ,e2)
ax.plot(r, psi, 'g')
ax.set_title("l = {} and E_blue = {}, E_green={}".format(l,e1,e2))
plt.show()

if name=="main": main()