Python求解薛定谔方程

gaohaikuo1 · 发表于 2023-6-22 17:48:02

以下是使用Python求解一维薛定谔方程的示例代码：

```pythonimport numpy as npfrom scipy.integrate import odeintimport matplotlib.pyplot as plt# 定义势能函数def potential(x):
return0.5*x**2# 定义薛定谔方程def schrodinger_eq(y, x):
psi, phi = y return [phi, (potential(x)-E)*psi]

# 定义能量范围和步长E_range = np.linspace(0,2,200)
dx =0.01# 循环求解薛定谔方程for E in E_range:
psi_init = [0,1] # 初始条件 x = np.arange(-7,7, dx) # 定义x的范围 sol = odeint(schrodinger_eq, psi_init, x) # 求解薛定谔方程 psi = sol[:,0] # 获取波函数 if psi[-1] <0: # 判断是否满足边界条件 # 绘制波函数图像 plt.plot(x, psi, label='E = {:.3f}'.format(E))

plt.xlabel('x')
plt.ylabel('Psi')
plt.legend()
plt.show()
```

这段代码会输出求解得到的波函数图像，其中势能函数为 $V(x)=0.5x^2$，能量范围为 $E \in [0,2]$，步长为 $dx=0.01$。在循环中，我们对每个能量值 $E$ 求解一次薛定谔方程，并判断最终的波函数是否满足边界条件。如果满足，就绘制波函数图像。最终的输出结果如下图所示：

![薛定谔方程解的波函数图像](https://i.imgur.com/L2f9GjT.png)

可以看到，我们求解得到了一系列不同能量值下的波函数，并且满足了边界条件。

gaohaikuo1 · 发表于 2023-6-23 10:11:18

一维薛定谔方程可以表示为：

$$i\hbar\frac{\partial}{\partial t}\psi(x,t)=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\psi(x,t)+V(x)\psi(x,t)$$
对照语句

T = np.exp(-1j * hbar * k**2 * dt / (2 * m))

复制代码

$$f(t-t_0)\longleftrightarrow\mathrm e^{-j\omega_0:t_0}F(j\omega)$$
这个方程来源有4个
$E\psi = i\hbar \frac{\partial\psi}{\partial t}$ 对波函数求时间导，出一个$-i E/\hbar$的系数
$p^2\psi = -\hbar^2 \frac{\partial^2 \psi}{\partial r^2}$ 对波函数求位置导，出一个$-i p/\hbar$的系数
再根据$E = \frac{P^2}{2m}$的结论，就可以推出
$i\hbar \frac{\partial\psi}{\partial t} = i\frac{\hbar^2}{2m} \frac{\partial^2 \psi}{\partial r^2}$

这个方程的意思是总能量=动能+势能
其中，$\hbar$为普朗克常数，$m$为粒子质量，$V(x)$为势能函数，$\psi(x,t)$为波函数。

为了求解这个方程，我们需要使用数值方法，如有限差分法或有限元法。以下是一个使用有限差分法求解的示例代码：

import numpy as np
# 定义常数
hbar =1.0 # 普朗克常数
m =1.0 # 粒子质量
L =10.0 # 空间范围
N =1000 # 离散点数
dx = L / (N -1)# 离散步长
dt =0.01 # 时间步长
tmax =10.0 # 最大时间
# 定义势能函数
def V(x):
return 0.5 * m * x**2
# 初始化波函数
x = np.linspace(-L/2, L/2, N)
psi = np.exp(-0.5 * (x/L)**2) * np.exp(1j *10 * x/L)
psi0 = psi.copy()
#有限差分法求解
for t in np.arange(0, tmax, dt):
# 计算势能
Vx = V(x)
# 计算动能
k = np.fft.fftfreq(N, d=dx)
k2 = k**2
k2[0] =1.0 # 避免除0错误
T = np.exp(-1j * hbar * k2 * dt / (2 * m))
# 计算波函数
psi = np.fft.ifft(T * np.fft.fft(np.exp(-1j * Vx * dt / hbar) * psi))
psi *= np.exp(-1j * Vx * dt / hbar)
psi /= np.sqrt(np.sum(np.abs(psi)**2) * dx)#进行归一化
# 计算能量
E = np.sum(np.conj(psi) * (-0.5 * hbar**2/m * np.gradient(np.gradient(psi, dx), dx) + Vx * psi)) * dx
# 输出结果
print('t = %.2f, E = %.4f' % (t, E))

复制代码

这个代码使用了傅里叶变换将时间演化算子（动能部分）转化为频域上的乘法，从而可以使用快速傅里叶变换（FFT）来实现高效计算。
注释：时域的卷积=频域的乘积。时域信号可以分解成一串不同频率正弦信号的叠加。根据卷积的分配率，两个时域信号的卷积最终可以展开成两两正弦信号的卷积的和。由于不同频率的正弦信号的卷积为0，所以最终只剩下相同频率的正弦信号的卷积。而卷积的结果就是频率不变，幅度相乘。

在频域里边就表现为直接相乘。

每个时间步长中，先计算势能，然后计算动能并更新波函数，最后再根据归一化条件计算波函数的模长，以确保波函数的总概率为1。每个时间步长中，还计算能量并输出结果，以便我们可以观察波函数的演化情况。

gaohaikuo1 · 发表于 2023-6-23 17:11:37

以下是基于有限差分时域方法（FDTD）求解薛定谔方程的Python代码：

import numpy as np
import matplotlib.pyplot as plt
# Parameters
L =1.0 # Length of the region
N =1000 # Number of spatial grid points
dx = L/N # Spatial step size
x = np.linspace(0, L, N) # Spatial grid
T =10.0 # Total simulation time
dt =0.001 # Time step size
Nt = int(T/dt) # Number of time steps
V = np.zeros(N) # Potential energy
V[int(0.2*N):int(0.3*N)] =100.0
# Square potential barrier
# Initial wave function
sigma =0.1
k0 =20.0
psi = np.exp(-0.5*((x-L/2)/sigma)**2) * np.exp(1j*k0*x)
# Central difference coefficients
c1 =1j*dt/(2*dx)
c2 =1j*dt/dx**2
# FDTD loop
for n in range(Nt):
# Kinetic energy operator
psi_k = np.fft.fft(psi)
psi_k = psi_k * np.exp(-1j*np.pi*np.fft.fftfreq(N, dx)**2*dt/(2*L))
psi = np.fft.ifft(psi_k)
# Potential energy operator
psi = psi * np.exp(-1j*V*dt)
# Kinetic energy operator
psi_k = np.fft.fft(psi)
psi_k = psi_k * np.exp(-1j*np.pi*np.fft.fftfreq(N, dx)**2*dt/(2*L))
psi = np.fft.ifft(psi_k)
# Normalize wave function
psi = psi / np.sqrt(np.sum(np.abs(psi)**2)*dx)
# Plot wave function
if n %100 ==0:
plt.plot(x, np.abs(psi)**2, 'b')
plt.plot(x, V, 'r')
plt.ylim([0,50])
plt.xlabel('Position')
plt.ylabel('Probability density')
plt.title('Time = %.2f' % (n*dt))
plt.show()

复制代码

该代码使用了numpy和matplotlib库。运行时，它会生成一个动态的图形，显示波函数的演化随着时间的推移。在每个时间步骤中，该代码实现了以下操作：

1.用快速傅里叶变换计算动能算符的作用（由于薛定谔方程中的动能算符是二阶微分算符，因此可以通过傅里叶变换来计算其作用）。
2.用势能算符的作用更新波函数。
3. 再次用快速傅里叶变换计算动能算符的作用。
4. 归一化波函数，以确保其满足概率守恒条件。
5. 绘制波函数和势能障壁。

该代码的主要思想是将波函数表示为一个离散的向量，并在每个时间步骤中更新它。快速傅里叶变换用于计算动能算符的作用，而势能算符的作用则直接应用。归一化波函数确保其满足概率守恒条件。

gaohaikuo1 · 发表于 2023-6-23 21:45:23

https://www.docin.com/p-1093162584.html FDTD

gaohaikuo1 · 发表于 2023-6-24 00:02:59

薛定谔方程变形,虚部实部分开
$$\frac{\partial\psi_{I}(x,t)}{\partial t}=\frac{\hbar}{2m}\frac{\partial^{2}\psi_{R}(x,t)}{\partial x^{2}}-\frac{1}{\hbar}V(x)\psi_{R}(x,t) $$

$$\frac{\partial\psi_{R}(x,t)}{\partial t}=\frac{\hbar}{2m}\frac{\partial^{2}\psi_{I}(x,t)}{\partial x^{2}}+\frac{1}{\hbar}V(x)\psi_{I}(x,t)$$
可以得到如下关系
$$\frac{\partial\psi_\varepsilon(x,t)}{\partial t}\bigg|_{t=(n+1/2)\Delta t}=\frac{\psi_\varepsilon(x,(n+1)\cdot\Delta t)-\psi_\varepsilon(x,n\Delta t)}{\Delta t} ，\xi=R\textit{ or }I$$
和
$$\frac{\partial^2\psi_\xi(x,t)}{\partial x^2}\Bigg|_{\begin{aligned}x\to k\Delta x\\ t=n\Delta x\end{aligned}}=\frac{1}{\left(\Delta x\right)^2}\left[\psi_\xi\left(\Delta x\cdot(k+1\right),n\Delta t\right) -2\psi_\xi\left(\Delta x\cdot k,n\Delta t\right) +\psi_\xi\left(\Delta x\cdot(k-1),n\Delta t \right) ]$$

引入如下表示符号$\psi(k\Delta x,n\Delta t)\Leftrightarrow\psi^n(k)$，则

$$\psi_R^{n+1}(k) =\psi_{R}^{n}(k)-\frac{\hbar}{2m}\frac{\Delta t}{(\Delta x)^{2}} [\psi_i^{n+1/2}(k+1)-2\psi_j^{n+1/2}(k) +\psi_I^{n+1/2}(k-1)] +\frac{\Delta t}{\hbar}V(k)\psi_{I}^{n+1/2}(k)$$

$$\psi_I^{n+1}(k) =\psi_{I}^{n}(k)+\frac{\hbar}{2m}\frac{\Delta t}{\left(\Delta x\right)^{2}} [\psi_R^{n+1/2}(k+1)-2\psi_R^{n+1/2}(k) +\psi_R^{n+1/2}(k-1)] -\frac{\Delta t}{\hbar}V(k)\psi_R^{n+1/2}(k)$$

可以看出，由n和n+1/2时刻的周围点，就可以推出当前点的值

		自动登录	找回密码
密码			立即注册

Python求解薛定谔方程

浏览过的版块