Boost logo

Boost :

Subject: Re: [boost] [units] - learning to use, continued :)
From: Matthias Schabel (boost_at_[hidden])
Date: 2011-11-30 12:18:30

> 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));
                // 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)
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>()),
        const quantity<energy_time> h_(hbar);
        const double Hn(boost::math::hermite(n,root<2>(m*omega/h_)*x)),
        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),
        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 :


Boost list run by bdawes at, gregod at, cpdaniel at, john at