求解谐波势阱中粒子的 Shrodinger 方程

Question

你好（这是我第一次在堆栈溢出中发帖），我正在尝试使用射手法计算谐波势中粒子的前 3 个能级

代码改编自 Mark Newman 在 Computational Physics 中的一个脚本，该脚本计算了盒子中粒子的基态。

这里是链接http://www-personal.umich.edu/~mejn/cp/programs/squarewell.py

这是我改编的代码

import numpy as np

from scipy.constants import m_e,hbar,elementary_charge

#define constants

vo = 50
a = 1e-11

#define limits of integration for adaptive runge-kutta method

xi = -10*a
xf = 10*a
N = 1000
dx = (xf-xi)/N

def V(x):#define harmonic potential well
    return vo*((x**2)/(a**2))

def f(r,x,E): #schrodinger's equation
    psi = r[0]
    phi = r[1]
    fpsi = psi
    fphi = (2*m_e/(hbar**2))*(V(x)-E)*psi
    return np.array([fpsi,fphi],float)

def solve(E): #calculates wave function for an energy E
    psi = 0.0
    phi = 1.0
    r = np.array([psi,phi],float)
    for x in np.arange(xi,xf,dx): #adaptive runge-kutta method
        k1 = dx*f(r,x,E)
        k2 = dx*f(r+0.5*k1,x+0.5*dx,E)
        k3 = dx*f(r+0.5*k2,x+0.5*dx,E)
        k4 = dx*f(r+k3,x+dx,E)
        r += (k1+2*k2+2*k3+k4)/6
    return r[0]

#finds the energy using secant method

E1 = 0.0
E2 = elementary_charge
psi2 = solve(E1)
target = elementary_charge/1000

while abs(E1-E2)>target:
    psi1,psi2 = psi2,solve(E2)
    E1,E2 = E2,E2-psi2*(E2-E1)/(psi2-psi1)
    print (E2/elementary_charge)

运行时出现此错误

RuntimeWarning: invalid value encountered in double_scalars

E1,E2 = E2,E2-psi2*(E2-E1)/(psi2-psi1)

我认为这意味着 psi2 和 psi1 靠得太近了，但我不太确定如何解决这个问题

Answer 1

是的！你说得对，价值观太接近了。您的代码返回 nan。这是因为除以零。

我建议使用校正因子。像max(delta, (psi2-psi1)) 这样的分母，​​delta 仍然可以是一个非常小的值，但它会阻止division by zero。