我以两种方式在数值上求解了薛定谔方程的以下两级时间相关哈密顿量：

import numpy as np

def H(t):
    return np.array([[t,0.5],[0.5,-t]],dtype="complex64")

首先，将状态视为波向量：

from scipy.integrate import solve_ivp
def schrodinger(t, X):
    dXdt = -1j * H(t).dot(X)
    return dXdt

T=80.
t_list = np.arange(-T, T, 2*T / 1000.0)
init_state= np.array([1.,0.],dtype="complex64")
solution = solve_ivp(fun=lambda t, X: schrodinger(t, X), t_span=[-T,T], y0=init_state, t_eval=t_list, method="RK45", vectorized=True)

其次，状态被视为（矢量化）密度矩阵：

def schrodinger_rho(t, X):
    unit = np.eye(2, 2, dtype="complex64")
    dXdt = -1j * (np.kron(unit, H(t)) - np.kron(H(t).T, unit)).dot(X)
    return dXdt
T=80.
t_list = np.arange(-T, T, 2*T / 1000.0)
init_state= np.array([1.,0.,0.,0.],dtype="complex64")
solution_rho = solve_ivp(fun=lambda t, X: schrodinger_rho(t, X), t_span=[-T,T], y0=init_state, t_eval=t_list, method="RK45", vectorized=True)

数值结果如下：

这是一个问题：为什么国家不能保持其统一性，比如说$P_0+P_1=1$，当状态被视为波向量？的值越大T，这种违反倾向越强。此外，当我使用method="BDF".

由于数值成本，我想将状态视为波向量。有什么办法可以改善这种现象吗？