对量子谐波建模 oscillator/SHM
Modelling a quantum harmonic oscillator/SHM
我需要帮助来弄清楚为什么 b) 的基态图看起来不对，这是完整的问题：（我认为完整地发布它会为我尝试使用的方法提供上下文）
(a) Consider a square potential well with ()=0 in between two infinitely high walls separated by a distance equal to the Bohr radius, i.e. for all x in the interval [0,] .
Write a function solve(energy, func) which takes the parameter energy and a python function </code> . This function should solve the Schrödinger ODE for the case described above and return only the final value of () at the boundary .</p></li>
<li><p>Write a script using the function <code>solve(energy, func) to calculate the ground state energy of an electron in this potential well in units of eV (i.e. divide result by elementary charge value). For the initial condition, see technical hint below. For the number of values to solve for () the value 1000 is recommended.
The result of your calculation should be a number between 134 eV and 135 eV. One of the tests will assess your solve(energy, func) function for a distorted potential well.
(b) Consider the harmonic potential
()=02/2, where 0 and =10−11 m are constants. Limit the infinite range of the variable to the interval [−10,10] with 0=50 eV.
The harmonic oscillator is known to have equidistant energy eigenvalues. Check that this is true, to the precision of your calculation, by calculating the ground state and the first 2 excited states. (Hint: the ground state has energy in the range 100 to 200 eV.)
In order to test your result, write a function result() which must return the difference of calculated energy eigenvalues in eV as a single number. Note that the test with the expected number is hidden and will test your result to a precision of ±1 eV.
Provide a plot title and appropriate axis labels, plot all three wave functions onto a single canvas using color-coded lines and also provide suitable axis limits in x and y to render all curves well visible.
Technical Hint: This is not an initial value problem for the Schrödinger ODE but a boundary value problem! This requires additional effort, a different method to the previous ODE exercises and is hence an 'unseen' problem.
Take a simple initial condition for both problems at 0=0 or 0=−10 , respectively: (0)=0 and (0)/=1 . Use that as a starting point to solve the ODE and vary the energy, , until a solution converges. The task is to evaluate the variation of the energy variable until the second boundary condition (for instance at L for exercise (a)) is satisfied since the first boundary condition is already satisfied due to the choice of initial condition.
This requires an initial guess for the energy interval in which the solution might be and a computational method for root finding. Search scipy for root finding methods and select one which does not require knowledge of the derivative. Root finding is appropriate here since the boundary condition to satisfy is ()=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 eV 到 200 eV 之间的所有能量值，图形形状看起来像这样。我不明白我哪里出错了。我已尝试尽可能多地测试我的代码。
您的函数代码或任务文本没有理由 boundary。如 a) 取边界条件右区间端的值，理想情况下这将是无穷大的值。
您的代码中还有另外两个问题，它们并不是您的错：
由于某些原因，可能不会出现在最新版本的 scipy 中，而且我无法在可用的源代码中找到，根查找器的调用 toms748，您使用的似乎缓存了该函数，而不是用较新的函数替换它。这意味着在 b) 部分中，根查找器调用仍然找到了一个根，但它是 a) 部分中 eigen 的根。只需检查函数值即可确认这一点。我建议使用带有初始括号间隔的更通用的界面，默认情况下使用布伦特方法的一个版本。
第二个问题是能量的尺度是极端的，因此在求根方法的设计参数之外。一种解决方案是操纵绝对和相对公差以适合您的问题域和范围尺度。另一个解决方案是将问题转化为具有更合理缩放的域，以便误差控制 method/heuristic 在其设计范围和默认容差内工作。
sol = optimize.root_scalar(lambda s: eigen2(s*e_el), bracket=(100,200))
groundh = sol.root
groundh_ = sol.root*e_el
完美工作（可能使用布伦特方法作为包围法的标准），找到 138.023972 eV 的基态能量并得到波形图
继续搜索，先是符号变化，然后是根，在50*m*e_el到50*(m+1)*e_el区间内，找到下一个状态
我需要帮助来弄清楚为什么 b) 的基态图看起来不对，这是完整的问题：（我认为完整地发布它会为我尝试使用的方法提供上下文）
(a) Consider a square potential well with ()=0 in between two infinitely high walls separated by a distance equal to the Bohr radius, i.e. for all x in the interval [0,] .
Write a function
solve(energy, func)which takes the parameter energy and a python function</code> . This function should solve the Schrödinger ODE for the case described above and return only the final value of () at the boundary .</p></li> <li><p>Write a script using the function <code>solve(energy, func)to calculate the ground state energy of an electron in this potential well in units of eV (i.e. divide result by elementary charge value). For the initial condition, see technical hint below. For the number of values to solve for () the value 1000 is recommended.The result of your calculation should be a number between 134 eV and 135 eV. One of the tests will assess your solve(energy, func) function for a distorted potential well.
(b) Consider the harmonic potential ()=02/2, where 0 and =10−11 m are constants. Limit the infinite range of the variable to the interval [−10,10] with 0=50 eV.
The harmonic oscillator is known to have equidistant energy eigenvalues. Check that this is true, to the precision of your calculation, by calculating the ground state and the first 2 excited states. (Hint: the ground state has energy in the range 100 to 200 eV.)
In order to test your result, write a function
result()which must return the difference of calculated energy eigenvalues in eV as a single number. Note that the test with the expected number is hidden and will test your result to a precision of ±1 eV.Provide a plot title and appropriate axis labels, plot all three wave functions onto a single canvas using color-coded lines and also provide suitable axis limits in x and y to render all curves well visible.
Technical Hint: This is not an initial value problem for the Schrödinger ODE but a boundary value problem! This requires additional effort, a different method to the previous ODE exercises and is hence an 'unseen' problem.
Take a simple initial condition for both problems at 0=0 or 0=−10 , respectively: (0)=0 and (0)/=1 . Use that as a starting point to solve the ODE and vary the energy, , until a solution converges. The task is to evaluate the variation of the energy variable until the second boundary condition (for instance at L for exercise (a)) is satisfied since the first boundary condition is already satisfied due to the choice of initial condition.
This requires an initial guess for the energy interval in which the solution might be and a computational method for root finding. Search scipy for root finding methods and select one which does not require knowledge of the derivative. Root finding is appropriate here since the boundary condition to satisfy is ()=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 eV 到 200 eV 之间的所有能量值，图形形状看起来像这样。我不明白我哪里出错了。我已尝试尽可能多地测试我的代码。
您的函数代码或任务文本没有理由 boundary。如 a) 取边界条件右区间端的值，理想情况下这将是无穷大的值。
您的代码中还有另外两个问题，它们并不是您的错：
由于某些原因，可能不会出现在最新版本的
scipy中，而且我无法在可用的源代码中找到，根查找器的调用toms748，您使用的似乎缓存了该函数，而不是用较新的函数替换它。这意味着在 b) 部分中，根查找器调用仍然找到了一个根，但它是 a) 部分中eigen的根。只需检查函数值即可确认这一点。我建议使用带有初始括号间隔的更通用的界面，默认情况下使用布伦特方法的一个版本。第二个问题是能量的尺度是极端的，因此在求根方法的设计参数之外。一种解决方案是操纵绝对和相对公差以适合您的问题域和范围尺度。另一个解决方案是将问题转化为具有更合理缩放的域，以便误差控制 method/heuristic 在其设计范围和默认容差内工作。
sol = optimize.root_scalar(lambda s: eigen2(s*e_el), bracket=(100,200))
groundh = sol.root
groundh_ = sol.root*e_el
完美工作（可能使用布伦特方法作为包围法的标准），找到 138.023972 eV 的基态能量并得到波形图
继续搜索，先是符号变化，然后是根，在50*m*e_el到50*(m+1)*e_el区间内，找到下一个状态