无限井问题的数值解

计算科学 有限差分 Python 特征值 量子力学
2021-12-15 18:56:48

我已经使用以下代码来实现它

import numpy as np
from scipy.linalg import eig
import matplotlib.pyplot as plt

W = 10 #width of well
n = 1000 #number of points excluding 0,W
hbar = 1;m=1 #atomic units
x = np.linspace(1,W-1,num=n) #exclusion of 0,W
h = x[1]-x[0]#step size

#Construction of the hamiltonian matrix
H = np.zeros([n,n],dtype=complex)
for i in range(n):
    H[i][i] = -2
    if(i-1>=0):
        H[i][i-1] = 1
    if(i+1<n):
        H[i][i+1] = 1
H = -1*(hbar**2)/(2*m*(h**2))*H 
E,V = eig(H)

plt.plot(x,np.abs(V[:][0]),'.')
plt.plot(x,np.abs(V[:][1]),'.')
plt.show()

结果： 在此处输入图像描述

我想知道为什么这段代码不够优化，无法重现接近分析解决方案的结果，以及错误的来源是什么。
编辑：
N=10,000 的输出

但是，它运行大约需要 45 分钟

1个回答

我认为主要问题可能与您正在使用的求解器有关。在这种情况下，哈密顿量（矩阵）是 Hermitian，它甚至是对称的，因为它是纯实数。您可以使用eigh而不是eig利用这一点。

此外，您不仅要删除第一个点和最后一个点，还要删除每端大小为 1 的间隔。接下来，我将向您展示一个包含此更改的片段。

import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as plt

W = 10 # width of well
n = 1000 # number of points
hbar = 1
m = 1
x = np.linspace(0, W, num=n+2)
h = x[1] - x[0] # step size

# Construction of the hamiltonian matrix
H = np.zeros([n, n])
for i in range(n):
    H[i, i] = -2
    if(i-1>=0):
        H[i, i-1] = 1
    if(i+1<n):
        H[i, i+1] = 1
H = -1*(hbar**2)/(2*m*(h**2))*H

# Solution
E, V = eigh(H, eigvals=(0, 10))

# Visualization
eigfun = np.zeros((n+2))
eigfun[1:-1] = V[:, 0]
plt.plot(x, np.abs(eigfun)**2)
eigfun[1:-1] = V[:, 1]
plt.plot(x, np.abs(eigfun)**2)
plt.show()

此外，您正在可视化|ψ|，但我认为最好是可视化|ψ|2.

在此处输入图像描述

最后，您可以利用矩阵的稀疏性。看到这个答案。

其它你可能感兴趣的问题
上一篇来自 scipy 的 solve_ivp 没有整合整个 tspan 范围 下一篇弹性方程混合形式的稳定有限元