标签：10 Phi 宇称 python 打靶 谐振子 hbar root sol1

1. 理论公式

一维谐振子薛定谔方程：

\[-\frac{\hbar^2}{2m} \frac{ d^2 }{d x^2} \psi(x) + \frac{1}{2}m \omega^2 x^2 \psi(x) = E \psi(x). \]

即

\[\frac{d^2}{dx^2} \psi(x) = [\frac{ m^2 \omega^2 x^2}{\hbar^2} - \frac{2mE}{\hbar^2} ] \psi(x). \]

可以设计一个长度量 \(x_0 = \sqrt{\hbar/(m\omega)}\)，用来让 \(x\) 失去量纲，然后整理方程，估计可以用 \(E/(\hbar \omega/2)\)　这个没有量纲的量，这样来求解方程，比较干净。

取 \(t = x/x_0\)，\(E'=E/(\hbar \omega/2)\), \(\Phi(t) = \psi(t x_0)\)，则可以把上面的方程整理为

\[\frac{ d^2 }{dt^2} \Phi(t) = (t^2 - E')\Phi(t), \]

这样就舒服了。我们知道，偶宇称解为 \(E' = 1, 5, 9, \cdots\)；奇宇称解为 \(E' = 3, 7, 11, \cdots\)。

2. 代码

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

# 一维谐振子： Phi''(t) = (t^2 - E' ) Phi(t)
def deriv(y,t,EE): # EE 即 E' = E/(hbar omega /2 )
    return [y[1], (t*t - EE) * y[0]]
    # y[0] 是 Phi, y[1] 是 Phi', 所以 [Phi, Phi'] 的导数是：
    # [Phi', Phi''] = [Phi', (t*t-E') Phi ]

tmax = 10; n = 1000
def shoot( EE): # 从 tmax 往　t=0 打靶
    t1 = np.linspace(tmax, 0, n); # 生成横坐标，从大到小
    sol1 = odeint(deriv, [1E-10, 0], t1, args=(EE,)) # ODE 从右往左推
    return sol1[-1,0] # 想要偶宇称，就打 sol1[-1,1]; 想要奇宇称，就打 sol1[-1,0]

root = fsolve(shoot, [-0.5]) # 得到本征能量
print("root=",root," shoot(root)=",shoot(root))

# 画个图，显示 [0, tmax] 上的波函数
t1 = np.linspace(tmax, 0, n); 
sol1 = odeint(deriv, [1E-10, 0], t1, args=(root,))
plt.plot(t1, sol1[:,0], label = "sol1"); plt.show();

代码思路：

• 选定 \(t=10\) 的点，取值 \(\Phi(10) = 10^{-10}, \Phi'(10) = 0\)。
• 调用 odeint 函数从 \(t=10\) 往 \(t=0\) 推；如果想要偶宇称解，就打靶 \(\Phi'(0)\)；如果想要奇宇称解，就打靶 \(\Phi(0)\)。
• 打靶即调节 \(E'\)，使得 \(\Phi'(0) \approx 0\)(偶宇称)，或者 \(\Phi(0) \approx 0\) （奇宇称）。

运行结果（以奇宇称第一激发态为例）：

root= [3.00000002]  shoot(root)= 0.0002180337905883789

标签：10,Phi,宇称,python,打靶,谐振子,hbar,root,sol1
From： https://www.cnblogs.com/luyi07/p/16748140.html