1 /*
2 Copyright 2011-2013 Karsten Ahnert
3 Copyright 2011-2013 Mario Mulansky
4
5 Distributed under the Boost Software License, Version 1.0.
6 (See accompanying file LICENSE_1_0.txt or
7 copy at http://www.boost.org/LICENSE_1_0.txt)
8 */
9
10
11 #include <iostream>
12 #include <vector>
13
14
15 /* WARNING: Compilation in debug mode might consume enormous memory
16 (e.g. ~2GB on gcc 4.4 )
17 */
18
19 // first increase fusion macro variables (now done by odeint itself)
20 //#define BOOST_FUSION_INVOKE_MAX_ARITY 15
21 //#define BOOST_RESULT_OF_NUM_ARGS 15
22
23 //[ units_define_basic_quantities
24
25 #include <boost/numeric/odeint.hpp>
26 #include <boost/numeric/odeint/algebra/fusion_algebra.hpp>
27 #include <boost/numeric/odeint/algebra/fusion_algebra_dispatcher.hpp>
28
29 #include <boost/units/systems/si/length.hpp>
30 #include <boost/units/systems/si/time.hpp>
31 #include <boost/units/systems/si/velocity.hpp>
32 #include <boost/units/systems/si/acceleration.hpp>
33 #include <boost/units/systems/si/io.hpp>
34
35 #include <boost/fusion/container.hpp>
36
37 using namespace std;
38 using namespace boost::numeric::odeint;
39 namespace fusion = boost::fusion;
40 namespace units = boost::units;
41 namespace si = boost::units::si;
42
43 typedef units::quantity< si::time , double > time_type;
44 typedef units::quantity< si::length , double > length_type;
45 typedef units::quantity< si::velocity , double > velocity_type;
46 typedef units::quantity< si::acceleration , double > acceleration_type;
47 typedef units::quantity< si::frequency , double > frequency_type;
48
49 typedef fusion::vector< length_type , velocity_type > state_type;
50 typedef fusion::vector< velocity_type , acceleration_type > deriv_type;
51 //]
52
53
54
55 //[ units_define_ode
56 struct oscillator
57 {
58 frequency_type m_omega;
59
oscillatoroscillator60 oscillator( const frequency_type &omega = 1.0 * si::hertz ) : m_omega( omega ) { }
61
operator ()oscillator62 void operator()( const state_type &x , deriv_type &dxdt , time_type t ) const
63 {
64 fusion::at_c< 0 >( dxdt ) = fusion::at_c< 1 >( x );
65 fusion::at_c< 1 >( dxdt ) = - m_omega * m_omega * fusion::at_c< 0 >( x );
66 }
67 };
68 //]
69
70
71 //[ units_observer
72 struct streaming_observer
73 {
74 std::ostream& m_out;
75
streaming_observerstreaming_observer76 streaming_observer( std::ostream &out ) : m_out( out ) { }
77
78 struct write_element
79 {
80 std::ostream &m_out;
write_elementstreaming_observer::write_element81 write_element( std::ostream &out ) : m_out( out ) { };
82
83 template< class T >
operator ()streaming_observer::write_element84 void operator()( const T &t ) const
85 {
86 m_out << "\t" << t;
87 }
88 };
89
90 template< class State , class Time >
operator ()streaming_observer91 void operator()( const State &x , const Time &t ) const
92 {
93 m_out << t;
94 fusion::for_each( x , write_element( m_out ) );
95 m_out << "\n";
96 }
97 };
98 //]
99
100
main(int argc,char ** argv)101 int main( int argc , char**argv )
102 {
103 // typedef dense_output_runge_kutta
104 // <
105 // controlled_runge_kutta
106 // <
107 // runge_kutta_dopri5< state_type , double , deriv_type , time_type , fusion_algebra >
108 // >
109 // > stepper_type;
110
111 //[ units_define_stepper
112 typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type > stepper_type;
113
114 state_type x( 1.0 * si::meter , 0.0 * si::meter_per_second );
115
116 integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , stepper_type() ) , oscillator( 2.0 * si::hertz ) ,
117 x , 0.0 * si::second , 100.0 * si::second , 0.1 * si::second , streaming_observer( cout ) );
118 //]
119
120 return 0;
121 }
122