为什么我的量子谐振子解会爆炸？
为什么我的量子谐振子解会爆炸？
提问于 2022-08-15 19:20:58
我用odeint来解QHO的能级(Griffith问题2.55)。
我正在从x=0到x=3进行积分。当我绘制结果时，我希望看到一半的高斯，尾巴会向正无穷大或负无穷远爆炸，这取决于我是将能量参数设置在有效能级之上还是低于有效能级。
相反，我的解决方案会立即扩展到正无穷大，不会显示出任何其他行为。
下面是我的代码，包括我在注释中导出的ODEs系统：
#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 为什么输出不符合我的预期？
编辑:这是代码的ascii版本，以防有人对unicode有问题：
#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;
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(
[](const State& q, State& dqdx, const double x) {
auto& dPsi_dx = dqdx[0];
auto& d2Psi_dx2 = dqdx[1];
const auto& psi = q[0];
dPsi_dx = q[1];
d2Psi_dx2 = (std::pow(mass*frequency*x/hbar, 2) - energy) * psi;
},
state_init,
x_init, x_final, x_step_init,
[](const auto& q, auto x) {
std::cout << x << ", " << q[0] << "," << q[1] << std::endl;
});
}回答 1
Stack Overflow用户
回答已采纳
发布于 2022-08-18 17:40:15
这不是编码错误，而是在转换方程时的错误。在注释中的最后一个等式之后，您需要另一个步骤。
⇒ ∂²Ψ/∂x² = (EΨ - ½mω²x²Ψ)/(-ℏ²/2m)
⇒ ∂²Ψ/∂x² = ((mωx)²-2mE/ℏ²)Ψ请注意，因子中的常量与代码中使用的常量不同。
从数学上讲，y''=C·(x²-a²)处于|x|<a的振荡状态和|x|>a的指数状态。振荡频率取决于系数C·(a²-x²)的大小，因此，当a较小时，频率较小，波长可能大于2a，因此不能保证根的存在。对于足够大的a，x=0周围的频率将变得足够大，从而保证接近x=0的零点。中间的某个位置是本征态的最低能量。
通过忽略这个因子，所使用的能量被移动到最低的本征态以下，因此只有指数行为是可见的。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/73365550
复制相关文章
相似问题
腾讯云开发者
