如何将python中的crank-nicolson方法应用于像薛定谔这样的波动方程

Question

我正在尝试在没有势场的盒子模拟中做一个粒子。我花了一些时间才发现简单的显式和隐式方法会破坏酉时间演化，所以我求助于应该是酉的曲柄尼科尔森。但是当我尝试它时，我发现它仍然不是这样。我不确定我错过了什么。我使用的公式是这样的：

其中 T 是 x 和 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"，而是

*"如果没有给出 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");