如何用C/C++数值求解二维玻色-爱因斯坦凝聚基态微分方程？

核心逻辑澄清 #

你已经把问题转化为线性方程组 A v* = b 的形式，这里的关键是非线性项 |vₙ|² 基于当前步 vₙ 计算，属于已知量——半隐式方法的核心就是对非线性项做显式处理，线性项做隐式处理，因此矩阵A在每一步迭代中是完全确定的，不存在无法分离 vₙ₊₁ 的问题。

具体实现步骤 #

• 空间离散与矩阵构建
  1. 对二维(x,y)域做网格离散，比如采用N×N的均匀网格，将波函数v展平为N²维向量（按行/列优先顺序）。
  2. 构建Laplacian矩阵：用二维5点有限差分格式离散，得到N²×N²的稀疏矩阵（每个网格点仅与相邻4个点作用，稀疏性极高）。
  3. 构建Lz矩阵：二维角动量算子 Lz = -i(x∂/∂y - y∂/∂x)，通过有限差分离散后转化为稀疏矩阵。
  4. 势场V直接构建为对角矩阵，对角元对应每个网格点的势场取值。
• 单步迭代流程
  1. 计算当前步 vₙ 的模平方 |vₙ|²，转化为对角矩阵（对角元为对应网格点的模平方值）。
  2. 组装矩阵A：A = Identity/dTau - 0.5*Laplacian + V + β*diag(|vₙ|²) - Ω*Lz，所有矩阵均为N²×N²规模。
  3. 构造右端项b：b = vₙ / dTau，为N²维复向量（波函数是复值）。
  4. 求解线性方程组 A v* = b：由于A是稀疏矩阵，必须用稀疏求解器，比如C/C++中的Eigen库SparseLU、SuperLU，或PETSc库，避免使用稠密矩阵造成内存溢出。
  5. 归一化得到 vₙ₊₁ = v* / ||v*||：BEC波函数需满足归一化条件，这一步是迭代的必要环节。
• 收敛判定
  重复迭代直到 ||vₙ₊₁ - vₙ|| 小于预设阈值（如1e-6），此时的v即为基态波函数。

C/C++实现关键注意事项 #

• 稀疏矩阵存储：若N取200，N²为40000，稠密矩阵会占用1.6e9个元素，完全不可行，必须用COO/CSR等稀疏格式存储矩阵。
• 复数类型支持：波函数为复值，所有矩阵、向量需使用std::complex<double>类型。
• 初始化策略：推荐用归一化的高斯函数作为初始波函数：v₀(x,y) = exp(-(x²+y²)/(2σ²))，σ可根据势场范围调整。
• 步长选择：dTau不宜过大，否则迭代易发散，建议从0.01开始尝试，逐步调整到合适值。

常见误区修正 #

你觉得无法分离 vₙ₊₁，是混淆了隐式与半隐式的逻辑：半隐式中非线性项|vₙ|²用的是当前迭代步的vₙ，而非下一步的vₙ₊₁，因此矩阵A是线性且已知的，直接解线性方程组即可推进迭代。

内容的提问来源于stack exchange，提问作者Pascal Gicquel