如何将python中的cranknolson方法应用于薛定谔波动方程

Question

我试着在一个没有势场的盒子里模拟粒子。我花了一些时间才发现简单的显式和隐式方法打破了酉时间演化，所以我求助于crank nicolson，它应该是酉的。但当我尝试它时，我发现它仍然不是这样。我不知道我错过了什么。。我使用的公式是：

{1美元^

其中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))

Answer 1

您依赖Toeplitz构造函数来生成一个对称矩阵，因此对角线下方的条目与对角线上方的条目相同。但是，scipy.linalg.toeplitz(c, r=None)的文档并不是说“transpose”，而是

*"If r is not given, r == conjugate(c) is assumed."

这样得到的矩阵是自伴的。在这种情况下，这意味着对角线上方的条目的符号被切换了。在

---

首先构造一个稠密矩阵，然后提取稀疏表示是没有意义的。使用^{}从一开始就将其构造为稀疏三对角矩阵

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");