But still I am quite happy today, because I managed to get a correctly working plots of this harmonic oscillator, but without units, unfortunately. One small problem is that I couldn't get to higher energies because double precision has limited factorials. I would need a working factorial of 100! which currently seems unlikely. And my currently maximum is only 25!
Using Stirling's approximation will help here; see below code. However, you will likely run into other problems with the Hermite polynomials at large n.
So, for your pleasure I am presenting you with a harmonic oscillator code :)
Unfortunately without units :) But if you will try to convert this to units, I might learn even more about using this library. Especially because dimension of wavefunction in 1 dimensional space is 1/sqrt(meter), and I have no idea how to declare this derived_dimension, this didn't work:
Here's a fully-unitized code for the 1D harmonic oscillator: #include <complex> #include <iostream> #include <limits> using namespace std; #include <boost/mpl/list.hpp> // use math constants #include <boost/math/constants/constants.hpp> // use math factorial #include <boost/math/special_functions/factorials.hpp> #include <boost/math/special_functions/hermite.hpp> #include <boost/typeof/std/complex.hpp> #include <boost/units/cmath.hpp> #include <boost/units/io.hpp> #include <boost/units/pow.hpp> #include <boost/units/systems/si.hpp> #include <boost/units/systems/si/prefixes.hpp> #include <boost/units/systems/si/codata/electron_constants.hpp> #include <boost/units/systems/si/codata/physico-chemical_constants.hpp> #include <boost/units/systems/si/codata/universal_constants.hpp> using namespace boost::units; using namespace boost::units::si; // MCS - need codata namespace to find hbar using namespace boost::units::si::constants::codata; // MCS - derived_dimension only works with integer powers; use make_dimension_list here instead typedef make_dimension_list<boost::mpl::list<dim<length_base_dimension,static_rational<-1,2> > > >::type quantum_wavefunction_1D_dimension; typedef unit<quantum_wavefunction_1D_dimension,si::system> quantum_wavefunction_1D; double factorial_to_minus_one_half(int n) { if (n <= 10) { return 1./sqrt(boost::math::factorial<double>(n)); } else { // use Stirling's approximation for n>10 const double pi(boost::math::constants::pi<double>()), logfac(log(2*pi*n)/2. + n*(log(n)-1)); return exp(-0.5*logfac); } } // this function returns 1/sqrt(meter) quantity<quantum_wavefunction_1D> quantum_harmonic_oscillator_wavefunction_1d(int n, const quantity<mass>& m, const quantity<frequency>& omega, const quantity<length>& x) { const double pi(boost::math::constants::pi<double>()), prefac(pow(2.,-n/2.)*factorial_to_minus_one_half(n)); const quantity<energy_time> h_(hbar); const double Hn(boost::math::hermite(n,root<2>(m*omega/h_)*x)), gaussian(exp(-m*omega*pow<2>(x)/(2.*h_))); return prefac*root<4>(m*omega/(pi*h_))*gaussian*Hn; } int main() { // MCS - change -1 -> +1 here otherwise loop doesn't run quantity<length> s(15*nano*meter), ds(0.1*nano*meter); for (quantity<length> x=-s ; x<=s ; x+=ds) { const quantity<temperature> T(273.15*kelvin); const quantity<frequency> omega(k_B*T/hbar.value()); std::cout << x.value() << "\t"; for (int n=0;n<25;++n) { const quantity<quantum_wavefunction_1D> psi(quantum_harmonic_oscillator_wavefunction_1d(n,m_e,omega,x)); std::cout << psi.value() << "\t"; } std::cout << std::endl; } return 0; } and PDFs for n=0..25 : Matthias