用Numerov方法和二分法求解氢原子

Question

给出的代码没有收敛到正确的特征值。请指出可能的问题并进行修正。我怀疑在二分法函数中存在问题。

import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import simps



xlower = 10**-15

xupper = 60

dr = 0.001

x = np.arange(xlower,xupper+ dr,dr)

G = x.shape[0]

delta_x = dr

def potential(x, l=0):

    V = l*(l+1)/2*x**2 - 1/x

    return V

def intervals():

    delta_E = 1

    E = 0

    energy_intervals = []

    while  E>-1:

        energy_intervals.append((E-delta_E, E))

        E = E - delta_E

    return energy_intervals


def Numerov(E, x,v, G):

    psis = [0] * G

    psis[1] = 10**-6

    for j in range(1,len(x)-1):

        m = 2*(1-5/12* delta_x**2 * 2*(E - v[j]))*psis[j]

        n = (1+1/12*delta_x**2*(E - v[j-1])*2)*psis[j-1]

        o = 1 + 1/12*delta_x**2 *(E - v[j+1])*2

        psis[j+1] = (m - n)/o



    return psis



def binary_search(function, left, right,x, v, G):

    print('left', left)

    print('right', right)

    precision = 10**-14

    iter = 0

    while abs(right - left) > precision and iter <= 10:

        mid = (left + right) / 2

        print('mid', mid)

        psis = function(mid, x, v, G)

        left_psis = function(left, x,v, G)

        print('psis', psis[-1])

        if np.isclose(psis[-1], 10**-6):

            return mid, psis

        elif left_psis[-1]*psis[-1] < 0:

            print('psis < 0')

            right = mid

            print('new right', left)

        else:

            print('psis > 0')

            left = mid

            print('new left', right)



        iter += 1

    if np.isclose(psis[-1], 10**-6):

        mid = (left + right) / 2

        return mid, psis



    return None, None



def main():

    n_samples = 1

    potentials = []

    eigenvectors = {}

    eigenvalues = {}

    for i in range(n_samples):

        eigenvectors[i] = []

        eigenvalues[i] = []

        V = potential(x)

        energy_intervals = intervals()





        for interval in energy_intervals:

            print(interval)

            eigenvalue, psis = binary_search(Numerov, interval[0], interval[1], x, V, G)



            if eigenvalue is not None:

                print("eigenvalue found")

                eigenvalues[i].append(eigenvalue)

                eigenvectors[i].append(psis)



        print(f"{i}th sample is done")





    return eigenvalues, eigenvectors



eigenvalues, eigenvectors= main()

我先尝试在实现二分法之前计算解决方案中的节点数量，但那也没有成功。代码没有收敛到正确的特征值。