我正在尝试在没有势场的盒子模拟中创建粒子。花了我一些时间来发现简单的显式和隐式方法破坏了单一时间的演化,因此我求助于应该是单一的crank-nicolson。但是当我尝试它时,我发现事实并非如此。我不确定我缺少什么。我使用的公式是这样的:
其中T是二阶导数wrt x的三对角Toeplitz矩阵,
系统简化为
A和B矩阵是:
我只是解决了使用稀疏模块的线性系统。数学是有道理的,我在一些论文中发现了相同的数字格式,因此使我相信我的代码就是问题所在。
到目前为止,这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import toeplitz
from scipy.sparse.linalg import spsolve
from scipy import sparse
# Spatial discretisation
N = 100
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
# Time discretisation
K = 10000
t = np.linspace(0, 10, K)
dt = t[1] - t[0]
alpha = (1j * dt) / (2 * (dx ** 2))
A = sparse.csc_matrix(toeplitz([1 + 2 * alpha, -alpha, *np.zeros(N-4)]), dtype=np.cfloat) # 2 less for both boundaries
B = sparse.csc_matrix(toeplitz([1 - 2 * alpha, alpha, *np.zeros(N-4)]), dtype=np.cfloat)
# Initial and boundary conditions (localized gaussian)
psi = np.exp((1j * 50 * x) - (200 * (x - .5) ** 2))
b = B.dot(psi[1:-1])
psi[0], psi[-1] = 0, 0
for index, step in enumerate(t):
# Within the domain
psi[1:-1] = spsolve(A, b)
# Enforce boundaries
# psi[0], psi[N - 1] = 0, 0
b = B.dot(psi[1:-1])
# Square integration to show if it's unitary
print(np.trapz(np.abs(psi) ** 2, dx))
最佳答案
您依靠Toeplitz构造函数生成对称矩阵,因此对角线以下的条目与对角线上方的条目相同。但是,scipy.linalg.toeplitz(c, r=None)
的文档不是说“转置”,而是
*“如果未给出r,则假定r == conjugate(c)。”
因此所得矩阵是自伴随的。在这种情况下,这意味着对角线上方的条目的符号已切换。
首先构造一个密集矩阵然后提取一个稀疏表示是没有意义的。使用scipy.sparse.diags
从一开始就将其构造为稀疏三对角矩阵
A = sparse.diags([ (N-3)*[-alpha], (N-2)*[1+2*alpha], (N-3)*[-alpha]], [-1,0,1], format="csc");
B = sparse.diags([ (N-3)*[ alpha], (N-2)*[1-2*alpha], (N-3)*[ alpha]], [-1,0,1], format="csc");