Решение уравнения Шредингера

Стационарное уравнение Шредингера для одномерной потенциальной ямы решается как краевая задача. Решение должно удовлетворять набору условий:

• Непрерывно и непрерывно дифференцируемо
• Ограничено
• Интегрируемо с квадратом

В свзи с этим решение УШ выполняют по следующему алгоритму:

• Задают значение энергии Е
• Выбирают интервал, на котором будет решаться УШ (точки положительной и отрицательной бесконечностей)
• Задают условия на левом и правом концах интервала (обычно значения функции и производной вблизи машинного 0)
• Решают задачу Коши для УШ слева и справа, так, чтобы абсциссы конечных точек совпадали.
• Масштабируем (умножаем на число) одну из функций, чтобы значения обеих функций в точке соприкосновения совпадали (непрерывность).
• Сравниваем производные слева и справа в точке сшивки. Если равны, то полученная функция - волновая функция, а энергия - энергетический уровень такой системы.
• Если производные не совпадают, то меняем значение энергии (например с некоторым шагом).

После получения решения нужно проверить, что шаг интегрирования и выбор «бесконечностей» не влияют на полученный результат.

Фрагменты кода.

Задаем силы:

def force(psi, E):
    """ force function [dx/dx, dpsi/dx, d2psi/dx2] for psi vector [x, psi, psi']
    """
    return np.array([1.0, psi[2], (V(psi[0])-E)*psi[1]])

Делаем один шаг:

def step(f, x, h, E):
    """ one step by Runge-Kutta method
    """
    kt1=f(x, E)*h
    kt2=f(x+kt1/2., E)*h
    kt3=f(x+kt2/2., E)*h
    kt4=f(x+kt3, E)*h
    xn=x+(kt1+2.*kt2+2.*kt3+kt4)/6.
    return xn

Решаем задачу Коши, результат - список точек (координата, функция, производная):

def solve(f, psi, xs, N, E): 
    """f - function, psi - initial conditions, xs - end point, N - number of points
    """
    h=(xs-psi[0])/N
    pm=[psi]
    for i in range(N):
        psi = step(f, psi, h, E )
        pm.append(psi)  
    return pm

Считаем функция справа и слева и находим вронскиан в точке сшивки:

def wronsk(E):
    """ calculate Wronskian for two parts of solution
    """
    Xm, Xp = -20.0, 20.0  # effective infinities 
    Pm, Pp = 1.0e-15, 1.0e-15 # value of function on effective infinities 

    N = 500
    xs = 0.0 #match point

    psi = np.array([-20., 1.0e-15, 1.0e-15*(V(-20)-E)])
    fleft = solve(force, psi, xs, N, E)

    psi = np.array([20., 1.0e-15, -1.0e-15*(V(20)-E)])
    fright = solve(force, psi, xs, N, E)

    return fleft[-1][2]*fright[-1][1]-fleft[-1][1]*fright[-1][2]

Поиск нуля функции методом деления отрезка пополам:

def zero(wronsk, El, d, delta):
    """ find zero of wronsk function in Er...Er+d range, with delta precision
    """
    wl=wronsk(El)
    Er=El+d
    wr=wronsk(Er)
    while (abs(Er-El)>delta):
      En=(El+Er)/2.
      wn=wronsk(En)
      if (wn == 0): 
         return En
      if (wl*wn<0): 
         Er = En
         wr = wn
      if (wn*wr<0):
         El = En
         wl = wn
    return (El+Er)/2.