I am using odeint to solve for the energy levels of the QHO (Griffiths problem 2.55).
I am integrating from x=0 to x=3. When I plot the results, I expect to see half of a gaussian with a tail that explodes towards positive or negative infinity, depending on whether I set the energy parameter to be above or below a valid energy level.
Instead, my solution blows up to positive infinity right away, and will not show any other behavior.
Here is my code, including my derivation of the system of ODEs in a comment:
#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;
});
}
Here is some example output:
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
Why does the output not match my expectations?
edit: here is an ascii version of the code in case anyone has issues with 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;
});
}
This is not a coding error but an error in transforming the equations. After the last equation you have in the comments, you need another step
⇒ ∂²Ψ/∂x² = (EΨ - ½mω²x²Ψ)/(-ℏ²/2m)
⇒ ∂²Ψ/∂x² = ((mωx)²-2mE/ℏ²)Ψ
Note that the constant in the factor is different from the one you use in code.
Mathematically, y''=C·(x²-a²) is in an oscillating regime for |x|<a and in an exponential regime for |x|>a. The oscillating frequency depends on the size of the factor C·(a²-x²), so if a is small the frequency is small and the wave length can be larger than 2a so that there is no guarantee for a root. For a large enough the frequency around x=0 will become large enough so that zeros close to x=0 are guaranteed. Somewhere in-between is the lowest energy of an eigenstate.
By omitting the factor the energy used was shifted below the lowest eigenstate, so that only the exponential behavior was visible.