我'm writing a code in Python that predicts the energy levels of Hydrogen which I will use as a template for research into quarkonium energy levels. I' m使用 scipy.integrate.odeint()
函数来解决Shroedinger方程,它可以在低至1818992_的较低能量水平下正常工作 . 我不需要超越它,但odeint会返回 Excess work done on this call (perhaps wrong Dfun type).
,这只会鼓励我扩展我能预测的内容 .
我正在使用的Shroedinger方程替换是:
u'' - (l*(l+1)/r**2 - 2mu_e(E-V_emag(r))) * u = 0
=>
u' = v
v' = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*u
然后我在它上面使用 scipy.integrate.odeint()
并迭代能量并使用我定义的其他函数来评估结果中的转折点和节点 . 我找到能量水平的方法是找到尽可能低的E值,其中转折点和节点的数量与它应该的相匹配;然后将 L
递增1并找到新的地面能量,例如如果 L=0
我会找到 n=1
能量,如果 L=3
,我会找到 n=2
能量 .
一旦代码增加到 L=7
,它就不会返回任何有用的内容 . r
的范围已经延长,但我已经阅读了关于雅各比人的内容 . 我已经确定了它们是什么,或者我还需要它们 . 有任何想法吗?
def v_emag(r):
v = -alpha/r
return v
def s_e(y,r,l,E): #Shroedinger equation for electromagntism
x = numpy.zeros_like(y)
x[0] = y[1]
x[1] = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*y[0]
return x
def i_s_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
r = numpy.arange(start,stop,step)
y = odeint(s_e,y0,r,args=(l,E))
return y
def inormalise_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
r = numpy.arange(start,stop,step)
f = i_s_e(l,E,start,stop,step)[:,0]
f2 = f**2
area = numpy.trapz(f2,x=r)
return f/(numpy.sqrt(area))
def inodes_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
x = i_s_e(l,E,start,stop,step)
r = numpy.arange(start,stop,step)
k=0
for i in range(len(r)-1):
if x[i,0]*x[i+1,0] < 0: #If u value times next u value <0,
k+=1 #crossing of u=0 has occured, therefore count node
return k
def iturns_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
x = i_s_e(l,E,start,stop,step)
r = numpy.arange(start,stop,step)
k = 0
for i in range(len(r)-1):
if x[i,1]*x[i+1,1] < 0: #If du/dr value times next du/dr value <0,
k=k+1 #crossing of du/dr=0, therefore a maximum/minimum
return k
l = 0
while l < 10: #The ground state for a system with a non-zero angular momentum will
E1 = -1.5e-08 #be the energy level of principle quantum number l-1, therefore
E3 = 0 #by changing l, we can find n by searching for the ground state
E2 = 0.5*(E1+E3)
i = 0
while i < 40:
N1 = inodes_e(l,E1)
N2 = inodes_e(l,E2)
N3 = inodes_e(l,E3)
T1 = iturns_e(l,E1)
T2 = iturns_e(l,E2)
T3 = iturns_e(l,E3)
if N1 != N2:# and T1 != T2: #Looks in lower half first, therefore will tend to ground state
E3 = E2
E2 = 0.5*(E1+E3)
elif N2 != N3:# and T2 != T3:
E1 = E2
E2 = 0.5*(E1+E3)
else:
print "Can't find satisfactory E in range"
break
i += 1
x = inormalise_e(l,E2)
if x[((l+1)**2)/0.005] > (x[2*((l+1)**2)/0.005]) and iturns_e(l,E2+1e-20)==1:
print 'Energy of state: n =',(l+1),'is: ',(E2*(10**9)),'eV'
l += 1
else:
E1 = E2+10e-20
1 回答
我不完全确定你的
while i<40:
循环是做什么的,所以如果我错了,也许你可以更正以下内容 .如果您希望此系统的波函数为某个
n, l
,则可以将能量计算为E = RH / n ^ 2,其中RH是里德堡常数,因此您无需计算节点数 . 如果确实需要对节点进行计数,那么对应于(n,l)
的数字是n-l-1
,因此您可以改变E并观察固定l的节点数量变化 .主要的问题,在我看来,你的r范围不够大,无法包含所有节点(对于大n~l),并且
odeint
不知道远离其他(非物理)渐近解决方案,psi~exp(cr),所以在某些条件下将psi发送到大于r的±无穷大 .如果它有用的话,这就是我想出来找到SE方程的数值解:你需要根据
n,l
改变r
-range但是要避免上述问题(例如,如果你要求n, l = 10, 9
,看看会发生什么) .