建模量子谐振子/SHM

Question

我需要帮助弄清楚为什么b)的基态图看起来不正确，以下是完整问题：（我认为将其完整发布可以为我尝试使用的方法提供背景）

（a）考虑一个方形势阱，两个无限高的壁之间的距离等于玻尔半径，即对于所有x在区间[0，]中，()= 0。编写一个函数solve（energy，func），该函数采用参数能量和python函数func。此函数应解决上述情况下的薛定谔ODE，并仅返回边界处()的最终值。 编写一个使用函数solve（energy，func）的脚本，以eV为单位计算电子在该势阱中的基态能量（即将结果除以元电荷值）。对于初始条件，请参见以下技术提示。建议解决()的值1000。 您的计算结果应该是134 eV到135 eV之间的数字。其中一个测试将评估您的solve（energy，func）函数是否为扭曲势阱。 （b）考虑谐振势 ()=02/2，其中0和=10-11 m是常数。将变量的无限范围限制在[-10,10]区间内，₀=50 eV。 已知谐振子具有等差能级。通过计算基态和前2个激发态来检查这一点，以您的计算精度为准。 （提示：基态的能量范围在100到200 eV之间。） 为了测试您的结果，编写一个函数result（），它必须返回以eV为单位的计算能量本征值之间的差异作为单个数字。请注意，测试所需的数字已隐藏，并将以±1 eV的精度测试您的结果。 提供绘图标题和适当的轴标签，使用彩色线将所有三个波函数绘制到单个画布上，并提供适当的x和y轴限制，以使所有曲线都可见。 技术提示：这不是薛定谔ODE的初始值问题，而是边界值问题！这需要额外的工作，不同于以前的ODE练习的方法，因此是一个“未见过”的问题。 对于两个问题，采用简单的初始条件，即₀= 0或₀=-10：(₀)=0并且(₀)/= 1。使用此作为解决ODE的起点，并变化能量，直到解决方案收敛。任务是评估能量变化，直到满足第二个边界条件（例如在(a)练习中的L），因为由于初始条件的选择，已经满足了第一个边界条件。 这需要对解决方案可能存在的能量间隔进行初始猜测，并使用根查找的计算方法。在scipy中搜索根查找方法，并选择不需要知道导数的方法。在这里，根查找是适当的，因为要满足的边界条件是()= 0。

量子物理学背景下，两个练习的边界条件是在每个电势边界处()=0。这些条件仅在特定的离散能量值处实现，称为能量本征值，其中最小的被称为基态能量。

m_el   = 9.1094e-31      # mass of electron in [kg]
hbar   = 1.0546e-34      # Planck's constant over 2 pi [Js]
e_el   = 1.6022e-19      # electron charge in [C]
L_bohr = 5.2918e-11      # Bohr radius [m]
V0 = 50*e_el
a = 10**(-11)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
from scipy.integrate import odeint
from scipy import optimize

def eqn(y, x, energy):                       #part a)
    y0 = y[1]
    y1 = -2*m_el*energy*y[0]/hbar**2

    return np.array([y0,y1])

x = np.linspace(0,L_bohr,1000)
def solve(energy, func):
    p0 = 0
    dp0 = 1
    init = np.array([p0,dp0])
    ysolve = odeint(func, init, x, args=(energy,))
    return ysolve[-1,0]

def eigen(energy):

    return solve(energy, eqn)

root_ = optimize.toms748(eigen,134*e_el,135*e_el)
root = root_/e_el

print('Ground state infinite square well',root,'eV')

intervalb = np.linspace(-10*a,10*a,1000)                    #part b)

def heqn(y, x2, energy):
    f0 = y[1]
    f1 = (2.0 * m_el / hbar**2) * (V0 * (x2**2/a**2) - energy) * y[0]

    return np.array([f0,f1])

def solveh(energy, func):
    ph0 = 0
    dph = 1
    init = np.array([ph0,dph])
    ysolve = odeint(func, init, intervalb, args=(energy,))
    return ysolve

def boundary(energy):        #finding the boundary V=E to apply the b.c
    f = a*np.sqrt(energy/V0)
    index = np.argmin(np.abs(intervalb-f))

    return index

def eigen2(energy):

    return solveh(energy,heqn)[boundary(energy),0]

groundh_ = optimize.toms748(eigen2,100*e_el,200*e_el)
groundh = groundh_/e_el


print('Ground state of Harmonic Potential:', groundh, 'eV')
plt.suptitle('Harmonic Potential Well')
plt.xlabel('x (a [pm])')
plt.ylabel('Psi(x)')

groundsol = solveh(groundh_,heqn)[:,0]
plt.plot(intervalb/a, groundsol)

对于能量在100电子伏到200电子伏之间的所有数值，图形的形状都是这样的。我不明白我做错了什么。我已经尽可能地测试了我的代码。

- user5056300

---

你正在以米为单位从-10到10进行绘图。要求你将x的范围限制在-100到100皮米之间（[-10a，10a]）。你不应该针对intervalb而是针对intervalb/a进行绘图，因为intervalb已经变化了从-10*a到10*a。除以a值可以得到10米。 - William Miller

我不知道为什么这个程序无法运行，但你可以轻松地使用VMC/DMC解决这个问题（特别是对于harm.oscillator）。在你的情况下，这两种算法都非常简单，可能需要更少的代码量（尤其是使用pandas）。我记得很多年前用Fortran编写过它... https://en.wikipedia.org/wiki/Diffusion_Monte_Carlo - Oleg O

我认为我做不到，因为问题似乎需要这种方法。 - user5056300

我已经编辑了这个问题，主要是为了让未来的读者更清楚哪些部分是你自己的声音，哪些部分是你的教育者的声音。我无法确定以“量子物理背景”开头的段落是你还是问题陈述 - 如果这不是你的创作，请编辑并引用此段。谢谢。 - halfer

Answer 1

您的代码或任务文本中没有关于函数boundary的原因。与a）相同，对于边界条件，请取右区间端点处的值，理想情况下应该是无穷远处的值。

您的代码中还有两个问题，并不是您的错：

• 由于某种原因，可能在最新版本的scipy中不存在，并且我无法在可用的源代码中找到，您使用的根查找器toms748似乎会缓存函数并将其替换为更新的函数。这意味着在b）部分中，根查找器调用仍然找到了根，但它是来自a）部分的eigen的根。只需检查函数值即可确认此情况。我建议使用初始括号间隔的更通用接口，默认情况下使用Brent方法的某个版本。

• 第二个问题是能量尺度极大，因此超出了根查找方法的设计参数范围。一种解决方案是操纵绝对和相对容差以适合您的问题域和范围比例。另一个解决方案是将问题转换为具有更合理比例的域，以便误差控制方法/启发式方法在其设计范围和默认容差内工作。

sol = optimize.root_scalar(lambda s: eigen2(s*e_el), bracket=(100,200)) 
groundh = sol.root
groundh_ = sol.root*e_el

使用 Brent 方法作为括号方法的标准，完美运行，可在 138.023972 电子伏特处找到基态能量，并生成以下波形图：

继续搜索，在区间 50*m*e_el 到 50*(m+1)*e_el 中先寻找符号变化再寻找根，可找到下一状态：

- Lutz Lehmann