Python - неверные значения из численного решения дифференциального уравнения и нелинейного уравнения

Question

Я пытался создать программу для решения уравнений Шредингера для одномерного гармонического осциллятора. Одна из вещей, которые мне нужно было сделать, - это создать решение Рунге-Кутты (РК) для дифференциального уравнения, что я и сделал. Мне также нужно было найти допустимые собственные значения для энергии, которая в конечном итоге оказалась нелинейным уравнением, для решения которого я решил использовать метод секущих. Но проблема в том, что я получаю неверные значения из решения уравнения Шредингера с РК и для энергий с помощью метода секанса. Волновая функция должна стремиться к нулю по мере приближения к плюсу или минусу бесконечности (которое я аппроксимирую как "большое" расстояние по сравнению с колодцем, и мое собственное значение для энергии (первое решение или основное состояние) должно быть между 100 эВ. и 200 эВ.Но когда я рисую его и получаю собственные значения, я получаю что-то, что не похоже на волну, и собственное значение около 10000

Я пробовал менять единицы измерения, которые использовал (я работаю в электронвольтах и ​​нанометрах), но все равно получаю смешные значения, которые не близки к тому, что должно быть. Я попробовал свою программу на более простой задаче (бесконечный квадратный колодец) и получил правильные значения. Также, когда я расширяю диапазон оценки моей волновой функции, собственные значения энергии будут увеличиваться, хотя они не должны значительно увеличиваться, чем больше вы увеличиваете диапазон оценки. Я также попробовал метод градиентного спуска и метод релаксации для нахождения собственных значений, но они все еще были далеки от точных аналитических результатов (по крайней мере на порядок).

import numpy as np

#constants in ev-nm units
hbarc = 197.326    #hbar times speed of light
mc = 5.112e5        #mass energy of electron
a = 0.01           #scaling factor of oscillator
N = 10000           #number of steps
h = a/N            #step size
V0 = 50            #potential energy scaling factor

#potential energy function
def v(x):
    return V0*(x**2)/(a**2)
#wave function differential equation
def f(r,x,E):
    psi = r[0]
    phi = r[1]
    fpsi = phi
    fphi = 2*mc/(hbarc**2)*(v(x) - E)*psi
    return np.array([fpsi, fphi], float)

#Runge-Kutta function
def RK(E):
    psi = 0             #state of the wave function
    phi = 1
    r = np.array([psi, phi],float)
    for x in np.linspace(-100*a,100*a,N):              #Evaluation Range of the wavefunction
        k1 = h*f(r, x, E)
        k2 = h*f(r+ 0.5*k1, x + 0.5*h, E)
        k3 = h*f(r + 0.5*k2, x + 0.5*h, E)
        k4 = h*f(r + k3, x + h, E)
        r += (k1 + 2*k2 + 2*k3 + k4)/6
    return r[0]

#Secant method for finding the eigenenergies of the system
while abs(E1-E2) > error:
    psi1 = psi2
    psi2 = RK(E2)
    E2 = E2 - psi2*(E2 - E1)/(psi2 - psi1)
    E1 = E2
print("E = ", E2, "ev")

Я также попробовал метод градиентного спуска и метод релаксации для нахождения собственных значений, но они все еще были далеки от точных аналитических результатов (по крайней мере на порядок).