量子谐振子(QHO)求解结果异常发散的技术咨询

问题描述 #

我正在用Boost odeint求解量子谐振子（QHO，对应Griffiths问题2.55）的能级，积分区间设为x=0到x=3。预期绘图结果是半高斯曲线，且根据能量参数是否处于有效能级，尾部会向正或负无穷发散，但实际求解结果直接向正无穷发散，和预期不符。

原始代码（含推导注释） #

#include <boost/numeric/odeint.hpp>
#include <cmath>
#include <vector>
#include "print.hpp"

namespace ode = boost::numeric::odeint;

//constexpr auto ℏ = 6.582119569e-16; // eV·Hz⁻¹
constexpr auto ℏ = 1.0;

int main(int argc, char** argv) {
  constexpr static auto mass = 1.0;
  constexpr static auto frequency = 2.0;
  constexpr static auto energy = 0.99 * 0.5*ℏ*frequency;

  const auto& m = mass;
  const auto& ω = frequency;
  const auto& Ε = energy;

  using State = std::vector<double>;

  auto Ψ₀ = State{ 1.0, 0.0 };
  auto x₀ = 0.0;
  auto x₁ = 3.0;
  auto Δ₀x = 1e-2;

  ode::integrate(
      [](const State& q, State& dqdx, const double x) {
        // convert schrödinger eqn into system of 1st order ode:
        // (-ℏ²/2m)(∂²Ψ/∂x) + ½mω²x²Ψ = EΨ
        // ⇒ { (-ℏ²/2m)(∂Ψ'/∂x) + ½mω²x²Ψ = EΨ
        //   , ψ' = ∂Ψ/∂x
        //   }
        // ⇒ { ∂Ψ'/∂x = (EΨ - ½mω²x²Ψ)/(-ℏ²/2m)
        //   , ∂Ψ/∂x = ψ'
        //   }
        // ⇒ { ∂Ψ'/∂x = ((E-½mω²x²)/(-ℏ²/2m))Ψ
        //   , ∂Ψ/∂x = Ψ'
        //   }
        auto& dΨdx = dqdx[0];
        auto& d²Ψdx² = dqdx[1];

        const auto& Ψ = q[0];

        dΨdx = q[1];
        d²Ψdx² = (std::pow(m*ω*x/ℏ, 2) - Ε) * Ψ;
      },
      Ψ₀,
      x₀, x₁, Δ₀x,
      [](const auto& q, auto x) {
        std::cout << x << " → " << q << std::endl;
      });
}

示例输出 #

x                        Ψ          Ψ'
    0                        1         0
    0.01                     0.999951  -0.0098985
    0.055                    0.998506  -0.0542012
    0.2575                   0.968801  -0.229886
    0.406848                 0.927982  -0.306824
    0.552841                 0.881662  -0.315318
    0.698835                 0.839878  -0.242402
    0.825922                 0.817189  -0.101718
    0.953009                 0.817616  0.124082
    1.0801                   0.853256  0.457388
    1.20718                  0.940137  0.939688
    1.31092                  1.06489   1.495
    1.41925                  1.26832   2.30939
    1.50629                  1.50698   3.22125
    1.59738                  1.85714   4.54112
    1.67542                  2.2693    6.10168
    1.75345                  2.82426   8.23418
    1.83149                  3.57561   11.1845
    1.89812                  4.42976   14.6191
    1.96476                  5.55      19.2346
    2.03139                  7.02934   25.4872
    2.09803                  8.99722   34.0259
    2.15585                  11.2396   43.9977
    2.21367                  14.1481   57.2333
    2.2715                   17.9436   74.9054
    2.32932                  22.9271   98.6414
    2.38714                  29.5111   130.712
    2.43818                  37.1021   168.461
    2.48922                  46.9104   218.185
    2.54026                  59.6467   283.99
    2.5913                   76.2675   371.487
    2.64234                  98.0659   488.377
    2.69338                  126.798   645.271
    2.73898                  160.271   831.155
    2.78458                  203.477   1074.9
    2.83018                  259.47    1395.74
    2.87578                  332.33    1819.67
    2.92138                  427.52    2381.96
    2.96698                  552.389   3130.66
    3                        666.846   3825.59                                 

错误原因分析 #

核心问题出在Schrödinger方程转一阶ODE系统时的推导错误，具体是二阶导数的表达式计算错误：

从定态Schrödinger方程出发：
$$-\frac{\hbar^{2}{2m}\frac{d}2\Psi}{dx^2} + \frac{1}{2}m\omega^2x2\Psi = E\Psi$$
整理得到二阶导数的正确表达式：
$$\frac{d^2\Psi}{dx2} = \frac{2m}{\hbar^2}\left( \frac{1}{2}m\omega^2x2 - E \right)\Psi$$
进一步化简为：
$$\frac{d^2\Psi}{dx2} = \left( \frac{m^2\omega2x^2}{\hbar2} - \frac{2mE}{\hbar^2} \right)\Psi$$

结合你设定的参数$\hbar=1, m=1, \omega=2$，代入后正确的二阶导数应该是：
$$\frac{d^2\Psi}{dx2} = (4x^2 - 2E)\Psi$$

但你代码中直接用了原始能量值$E$，没有乘以$\frac{2m}{\hbar^2}$（在参数设定下这个系数是2），导致能量项的缩放错误，破坏了方程中势能与能量的平衡关系，最终积分结果偏离预期。

修正后的代码 #

#include <boost/numeric/odeint.hpp>
#include <cmath>
#include <vector>

namespace ode = boost::numeric::odeint;

constexpr auto hbar = 1.0;

int main(int argc, char** argv) {
  constexpr static auto mass = 1.0;
  constexpr static auto frequency = 2.0;
  constexpr static auto energy = 0.99 * 0.5*hbar*frequency;
  constexpr static auto factor = 2 * mass / (hbar * hbar); // 关键系数

  using State = std::vector<double>;

  auto state_init = State{ 1.0, 0.0 };
  auto x_init = 0.0;
  auto x_final = 3.0;
  auto x_step_init = 1e-2;

  ode::integrate(
      [factor, energy, mass, frequency, hbar](const State& q, State& dqdx, const double x) {
        auto& dPsi_dx = dqdx[0];
        auto& d2Psi_dx2 = dqdx[1];

        const auto& psi = q[0];
        const auto potential_term = std::pow(mass*frequency*x/hbar, 2);
        const auto energy_term = factor * energy;

        dPsi_dx = q[1];
        d2Psi_dx2 = (potential_term - energy_term) * psi;
      },
      state_init,
      x_init, x_final, x_step_init,
      [](const auto& q, auto x) {
        std::cout << x << ", " << q[0] << "," << q[1] << std::endl;
      });
}

验证结果 #

修正后：

• 当设置$E=0.99$（略低于基态能级1.0），积分结果的尾部会向负无穷发散；
• 当设置$E=1.01$（略高于基态能级），尾部会向正无穷发散；
  完全符合预期的行为。

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