...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
Home | Libraries | People | FAQ | More |
odeint also works well with Boost.Units
- a library for compile time unit and dimension analysis. It works by decoding
unit information into the types of values. For a one-dimensional unit you
can just use the Boost.Unit types as state type, deriv type and time type
and hand the vector_space_algebra
to the stepper definition and everything works just fine:
typedef units::quantity< si::time , double > time_type; typedef units::quantity< si::length , double > length_type; typedef units::quantity< si::velocity , double > velocity_type; typedef runge_kutta4< length_type , double , velocity_type , time_type , vector_space_algebra > stepper_type;
If you want to solve more-dimensional problems the individual entries typically
have different units. That means that the state_type
is now possibly heterogeneous, meaning that every entry might have a different
type. To solve this problem, compile-time sequences from Boost.Fusion
can be used.
To illustrate how odeint works with Boost.Units we use the harmonic oscillator as primary example. We start with defining all quantities
#include <boost/numeric/odeint.hpp> #include <boost/numeric/odeint/algebra/fusion_algebra.hpp> #include <boost/numeric/odeint/algebra/fusion_algebra_dispatcher.hpp> #include <boost/units/systems/si/length.hpp> #include <boost/units/systems/si/time.hpp> #include <boost/units/systems/si/velocity.hpp> #include <boost/units/systems/si/acceleration.hpp> #include <boost/units/systems/si/io.hpp> #include <boost/fusion/container.hpp> using namespace std; using namespace boost::numeric::odeint; namespace fusion = boost::fusion; namespace units = boost::units; namespace si = boost::units::si; typedef units::quantity< si::time , double > time_type; typedef units::quantity< si::length , double > length_type; typedef units::quantity< si::velocity , double > velocity_type; typedef units::quantity< si::acceleration , double > acceleration_type; typedef units::quantity< si::frequency , double > frequency_type; typedef fusion::vector< length_type , velocity_type > state_type; typedef fusion::vector< velocity_type , acceleration_type > deriv_type;
Note, that the state_type
and the deriv_type
are now
a compile-time fusion sequences. deriv_type
represents x' and is now different from the state type
as it has different unit definitions. Next, we define the ordinary differential
equation which is completely equivalent to the example in Harmonic
Oscillator:
struct oscillator { frequency_type m_omega; oscillator( const frequency_type &omega = 1.0 * si::hertz ) : m_omega( omega ) { } void operator()( const state_type &x , deriv_type &dxdt , time_type t ) const { fusion::at_c< 0 >( dxdt ) = fusion::at_c< 1 >( x ); fusion::at_c< 1 >( dxdt ) = - m_omega * m_omega * fusion::at_c< 0 >( x ); } };
Next, we instantiate an appropriate stepper. We must explicitly parametrize
the stepper with the state_type
,
deriv_type
, time_type
.
typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type > stepper_type; state_type x( 1.0 * si::meter , 0.0 * si::meter_per_second ); integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , oscillator( 2.0 * si::hertz ) , x , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , streaming_observer( cout ) );
Note | |
---|---|
When using compile-time sequences, the iteration over vector elements is
done by the |
It is quite easy but the compilation time might take very long. Furthermore, the observer is defined a bit different
struct streaming_observer { std::ostream& m_out; streaming_observer( std::ostream &out ) : m_out( out ) { } struct write_element { std::ostream &m_out; write_element( std::ostream &out ) : m_out( out ) { }; template< class T > void operator()( const T &t ) const { m_out << "\t" << t; } }; template< class State , class Time > void operator()( const State &x , const Time &t ) const { m_out << t; fusion::for_each( x , write_element( m_out ) ); m_out << "\n"; } };
Caution | |
---|---|
Using Boost.Units works nicely but compilation can be very time and memory consuming. For example the unit test for the usage of Boost.Units in odeint take up to 4 GB of memory at compilation. |
The full cpp file for this example can be found here harmonic_oscillator_units.cpp.