Численное решение уравнения Шредингера - собственные значения

Я рад, что наконец могу ответить на этот вопрос.

Метод Нумерова, описанный в Википедии, здесь не подходит . Чтобы дать вам представление о том, как действовать, давайте начнем с упрощенной версии метода. Что я собираюсь сделать, так это просто наивно дискретизировать дифференциальный оператор, например так:

Последнее уравнение — это просто определение — я рассматриваю пространство как решетку с точками.$d$врозь и я звоню$\psi_n$значение волновой функции на$n$й точке.. Теперь уравнение Шредингера читается в этих обозначениях:

Но это матричное уравнение! Позвольте мне быть более явным:

Конечно, я должен был выбрать какое-то целое число$N$это просто соответствует помещению системы в ящик, который «достаточно велик», чтобы иметь конечную систему.

Теперь ясно, что у вас есть задача на собственные значения матрицы вида

и вы можете продолжать диагонализовать его любым способом, который вы выберете. Обратите внимание, что мы называем$A$трехдиагональная матрица по понятным причинам. Возможно, вы захотите сначала позаботиться о граничных условиях — вы сделаете это, установив$\psi_1 = \psi_n = 0$прежде чем построить матрицу, которая соответствует установке нуля в первом и последнем столбцах. Это даст вам несколько ложных собственных функций с нулевым собственным значением, которые вы можете просто отбросить. Если у вас разные граничные условия, вам не повезло — я не знаю, как заставить это работать вообще.

Теперь вам просто нужно повторить вышеописанное с полным методом Нумерова, который будет немного сложнее, и все готово.

Кажется, это то, что вы ищете, но, конечно, это не единственный способ сделать это. Гриффитс описывает один из них, называемый «вилять собакой», который состоит из следующего: вы угадываете начальное значение энергии и вычисляете волновую функцию, как хотите (например, RK). Скорее всего, это не будет собственным значением, поэтому волновая функция взорвется на бесконечности. Это может пойти на$+\infty$для больших$x$или это может пойти на$-\infty$. Теперь вы медленно меняете энергию до тех пор, пока поведение на бесконечности не «перевернется», то есть пока хвост не «виляет». Это позволит вам ограничить значение одного собственного значения энергии и даст вам форму волновой функции до некоторого максимального размера, который увеличивается с выбранным значением.$E$приближается к точному собственному значению.

Вы можете искать собственные значения, используя метод деления пополам.

Предварительные сведения: чтобы получить собственные значения методом Нумерова, вам необходимо знать волновую функцию на границах. Обычно это означает, что вам нужно установить потенциал на бесконечность на границах, следовательно, обнулить волновую функцию в этих точках. Для вашего потенциала измените его следующим образом:

V = бесконечность, если -50 фм>r>50 фм V = 0, если |r|>8,5 фм V = потенциал Вуда в противном случае

Суть дела: теперь напишите программу для расчета волновой функции в приведенном выше потенциале, используя любую произвольную энергию в области -50 фм.

Вот код Python, который находит собственное значение в заданном интервале (E1, E2), если оно существует.

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()