时间戳

首页 > 解决方案 > 求解谐波势阱中粒子的薛定谔方程

python - 求解谐波势阱中粒子的薛定谔方程

问题描述

你好（这是我第一次在堆栈溢出中发帖），我正在尝试使用射手方法计算谐波势中粒子的前 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 靠得太近了，但我不太确定如何解决这个问题

标签： pythonnumpy

解决方案

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

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

推荐阅读