薛定谔方程实战：用线性代数理解量子态演化（含Matlab/Julia代码片段）
薛定谔方程实战：用线性代数理解量子态演化（含Matlab/Julia代码片段）
量子力学的数学框架本质上是一个线性代数问题——这个观点可能会让初学者感到惊讶。当我们剥离那些令人望而生畏的物理概念，剩下的核心其实是如何在无限维的希尔伯特空间中处理向量和矩阵运算。本文将带您从计算数学的角度重新审视薛定谔方程，通过可执行的代码示例揭示量子态演化的矩阵本质。
1. 希尔伯特空间中的矩阵表述
量子态的本质是一个状态向量，存在于抽象的希尔伯特空间中。这个空间与我们熟悉的三维欧几里得空间有着惊人的相似性——只是维度可能无限，且向量的分量可以是复数。让我们用一个简单的两能级系统来说明：
% 两能级系统的量子态表示 psi_up = [1; 0]; % 自旋向上态 psi_down = [0; 1]; % 自旋向下态 psi_superposition = (psi_up + psi_down)/sqrt(2); % 叠加态在坐标表象下，波函数ψ(x)实际上就是这个抽象态向量在位置基矢上的"分量"。这种对应关系可以通过狄拉克符号清晰地表达：
| 概念 | 线性代数类比 | 量子力学表述 |
|---|---|---|
| 向量 | v ∈ ℂⁿ | |
| 基矢 | {e₁,e₂,...,eₙ} | { |
| 内积 | ⟨u,v⟩ = ∑uᵢ*vᵢ | ⟨φ |
| 线性算子 | 矩阵A | 算符Â |
提示：在数值计算中，我们总是需要将无限维的希尔伯特空间截断为有限维，这是量子计算模拟的基础近似。
2. 表象变换的数值实现
坐标表象和动量表象之间的转换，本质上是一个基变换问题。在离散化处理后，这个变换可以通过傅里叶变换矩阵来实现：
using FFTW function position_to_momentum(psi_x, dx) N = length(psi_x) dp = 2π/(N*dx) # 动量空间格点间距 psi_p = fftshift(fft(psi_x)) * dx / sqrt(2π) return psi_p, dp end # 示例：高斯波包变换 x = range(-5, 5, length=512) dx = x[2] - x[1] psi_x = exp.(-x.^2 / 2) .* exp.(1im * 2π * x) # 带有动量偏移的高斯波包 psi_p, dp = position_to_momentum(psi_x, dx)这个简单的代码演示了如何将波函数从位置表象转换到动量表象。值得注意的是，离散傅里叶变换引入的周期性边界条件在实际物理问题中需要特别处理。
3. 定态薛定谔方程的矩阵解法
对于势场V(x)中的粒子，其哈密顿量在位置表象下可以离散化为一个稀疏矩阵。我们采用有限差分法来近似二阶导数：
function H = construct_hamiltonian(x, V) N = length(x); dx = x(2) - x(1); e = ones(N,1); % 动能项（三对角矩阵） K = -1/(2*dx^2) * spdiags([e -2*e e], -1:1, N, N); % 势能项（对角矩阵） U = spdiags(V(:), 0, N, N); H = K + U; % 总哈密顿量 end % 方势阱示例 x = linspace(-5, 5, 501); V = double(abs(x) > 2) * 10; % 高度为10的方势垒 H = construct_hamiltonian(x, V); % 求解本征态 [psi, E] = eigs(H, 10, 'smallestreal');这种方法将微分方程转化为矩阵特征值问题，是现代计算量子力学的基础。下表比较了不同数值方法的特点：
| 方法 | 精度 | 内存效率 | 适用场景 |
|---|---|---|---|
| 有限差分法 | 中等 | 高 | 一般势场 |
| 谱方法 | 高 | 低 | 周期性系统 |
| 有限元法 | 可变 | 中等 | 复杂几何 |
| 离散变量表示法 | 很高 | 低 | 精确对角化 |
4. 含时演化的数值模拟
量子态的动力学演化由含时薛定谔方程控制，可以通过指数算子来实现数值求解：
using LinearAlgebra function time_evolution(psi0, H, tlist; dt=0.01) U = exp(-1im * H * dt) # 时间演化算子 psi_t = [psi0] for t in tlist[2:end] push!(psi_t, U * psi_t[end]) end return psi_t end # 谐振子初始态 x = range(-5, 5, length=512) psi0 = exp.(-(x.+1).^2 / 0.5) # 偏离中心的高斯波包 H = construct_hamiltonian(x, 0.5 * x.^2) # 谐振子势能 t = range(0, 10, length=100) psi_t = time_evolution(psi0, H, t)对于大系统，更高效的算法如Crank-Nicolson方法或分裂算子法可能更合适：
分裂算子法步骤：
- 将哈密顿量分为动能和势能部分：H = T + V
- 近似时间演化算子：exp(-iHΔt) ≈ exp(-iVΔt/2)exp(-iTΔt)exp(-iVΔt/2)
- 在动量空间计算动能部分，在位置空间计算势能部分
误差控制：
- 时间步长Δt需要满足Δt ≪ 1/ΔE，其中ΔE是能级间距
- 空间离散化dx应远小于波函数振荡的特征长度
量子模拟中最令人惊叹的或许是，即使在这些数值近似下，我们仍然能够观察到量子干涉、隧穿等纯粹量子现象的数值再现。比如在双势阱系统中运行上述代码，可以清晰地看到波函数在两个势阱间的周期性振荡——这正是量子叠加原理的直观体现。