libs/numeric/odeint/examples/harmonic_oscillator.cpp

   1 /*
   2  Copyright 2010-2012 Karsten Ahnert
   3  Copyright 2011-2013 Mario Mulansky
   4  Copyright 2013 Pascal Germroth
   5
   6  Distributed under the Boost Software License, Version 1.0.
   7  (See accompanying file LICENSE_1_0.txt or
   8  copy at http://www.boost.org/LICENSE_1_0.txt)
   9  */
  10
  11
  12 #include <iostream>
  13 #include <vector>
  14
  15 #include <boost/numeric/odeint.hpp>
  16
  17
  18
  19 //[ rhs_function
  20 /* The type of container used to hold the state vector */
  21 typedef std::vector< double > state_type;
  22
  23 const double gam = 0.15;
  24
  25 /* The rhs of x' = f(x) */
  26 void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
  27 {
  28     dxdt[0] = x[1];
  29     dxdt[1] = -x[0] - gam*x[1];
  30 }
  31 //]
  32
  33
  34
  35
  36
  37 //[ rhs_class
  38 /* The rhs of x' = f(x) defined as a class */
  39 class harm_osc {
  40
  41     double m_gam;
  42
  43 public:
  44     harm_osc( double gam ) : m_gam(gam) { }
  45
  46     void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
  47     {
  48         dxdt[0] = x[1];
  49         dxdt[1] = -x[0] - m_gam*x[1];
  50     }
  51 };
  52 //]
  53
  54
  55
  56
  57
  58 //[ integrate_observer
  59 struct push_back_state_and_time
  60 {
  61     std::vector< state_type >& m_states;
  62     std::vector< double >& m_times;
  63
  64     push_back_state_and_time( std::vector< state_type > &states , std::vector< double > &times )
  65     : m_states( states ) , m_times( times ) { }
  66
  67     void operator()( const state_type &x , double t )
  68     {
  69         m_states.push_back( x );
  70         m_times.push_back( t );
  71     }
  72 };
  73 //]
  74
  75 struct write_state
  76 {
  77     void operator()( const state_type &x ) const
  78     {
  79         std::cout << x[0] << "\t" << x[1] << "\n";
  80     }
  81 };
  82
  83
  84 int main(int /* argc */ , char** /* argv */ )
  85 {
  86     using namespace std;
  87     using namespace boost::numeric::odeint;
  88
  89
  90     //[ state_initialization
  91     state_type x(2);
  92     x[0] = 1.0; // start at x=1.0, p=0.0
  93     x[1] = 0.0;
  94     //]
  95
  96
  97
  98     //[ integration
  99     size_t steps = integrate( harmonic_oscillator ,
 100             x , 0.0 , 10.0 , 0.1 );
 101     //]
 102
 103
 104
 105     //[ integration_class
 106     harm_osc ho(0.15);
 107     steps = integrate( ho ,
 108             x , 0.0 , 10.0 , 0.1 );
 109     //]
 110
 111
 112
 113
 114
 115     //[ integrate_observ
 116     vector<state_type> x_vec;
 117     vector<double> times;
 118
 119     steps = integrate( harmonic_oscillator ,
 120             x , 0.0 , 10.0 , 0.1 ,
 121             push_back_state_and_time( x_vec , times ) );
 122
 123     /* output */
 124     for( size_t i=0; i<=steps; i++ )
 125     {
 126         cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
 127     }
 128     //]
 129
 130
 131
 132
 133
 134
 135
 136     //[ define_const_stepper
 137     runge_kutta4< state_type > stepper;
 138     integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 139     //]
 140
 141
 142
 143
 144     //[ integrate_const_loop
 145     const double dt = 0.01;
 146     for( double t=0.0 ; t<10.0 ; t+= dt )
 147         stepper.do_step( harmonic_oscillator , x , t , dt );
 148     //]
 149
 150
 151
 152
 153     //[ define_adapt_stepper
 154     typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
 155     //]
 156
 157
 158
 159     //[ integrate_adapt
 160     typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
 161     controlled_stepper_type controlled_stepper;
 162     integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 163     //]
 164
 165     {
 166     //[integrate_adapt_full
 167     double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
 168     controlled_stepper_type controlled_stepper(
 169         default_error_checker< double , range_algebra , default_operations >( abs_err , rel_err , a_x , a_dxdt ) );
 170     integrate_adaptive( controlled_stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 171     //]
 172     }
 173
 174
 175     //[integrate_adapt_make_controlled
 176     integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) ,
 177                         harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 178     //]
 179
 180
 181
 182
 183     //[integrate_adapt_make_controlled_alternative
 184     integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) ,
 185                         harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 186     //]
 187
 188     #ifdef BOOST_NUMERIC_ODEINT_CXX11
 189     //[ define_const_stepper_cpp11
 190     {
 191     runge_kutta4< state_type > stepper;
 192     integrate_const( stepper , []( const state_type &x , state_type &dxdt , double t ) {
 193             dxdt[0] = x[1]; dxdt[1] = -x[0] - gam*x[1]; }
 194         , x , 0.0 , 10.0 , 0.01 );
 195     }
 196     //]
 197
 198
 199
 200     //[ harm_iterator_const_step]
 201     std::for_each( make_const_step_time_iterator_begin( stepper , harmonic_oscillator, x , 0.0 , 0.1 , 10.0 ) ,
 202                    make_const_step_time_iterator_end( stepper , harmonic_oscillator, x ) ,
 203                    []( std::pair< const state_type & , const double & > x ) {
 204                        cout << x.second << " " << x.first[0] << " " << x.first[1] << "\n"; } );
 205     //]
 206     #endif
 207
 208
 209
 210
 211 }